25 #if defined (HAVE_CONFIG_H)
55 : q (q_arg), r (r_arg)
63 if (! (q_nc == r_nr && (q_nr == q_nc || (q_nr > q_nc && r_nr == r_nc))))
64 (*current_liboctave_error_handler) (
"QR dimensions mismatch");
73 if (! q.is_empty () && q.is_square ())
75 else if (q.rows () > q.columns () && r.is_square ())
103 #if ! defined (HAVE_QRUPDATE)
110 static bool warned =
false;
114 (*current_liboctave_warning_with_id_handler)
115 (
"Octave:missing-dependency",
116 "In this version of Octave, QR & Cholesky updating routines "
117 "simply update the matrix and recalculate factorizations. "
118 "To use fast algorithms, link Octave with the qrupdate library. "
119 "See <http://sourceforge.net/projects/qrupdate>.");
125 template <
typename T>
134 if (u.numel () != m || v.numel () != n)
137 init (q*r + T (u) * T (v).hermitian (), get_type ());
140 template <
typename T>
149 if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
152 init (q*r + u * v.hermitian (), get_type ());
155 template <
typename T,
typename CV_T>
160 T
retval (a.rows (), a.columns () + 1);
169 template <
typename T,
typename RV_T>
174 T
retval (a.rows () + 1, a.columns ());
183 template <
typename T>
193 template <
typename T>
203 template <
typename T>
225 template <
typename T>
237 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
239 init (octave::math::insert_col (q*r, j, u), get_type ());
242 template <
typename T>
256 dups = dups && js(i) == js(i+1);
259 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
260 if (u.numel () != m || u.columns () != nj)
262 if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
263 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
269 a = octave::math::insert_col (a, js(i), u.
column (i));
270 init (a, get_type ());
274 template <
typename T>
282 if (j < 0 || j > n-1)
283 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
285 init (octave::math::delete_col (q*r, j), get_type ());
288 template <
typename T>
301 dups = dups && js(i) == js(i+1);
304 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
305 if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
312 a = octave::math::delete_col (a, js(i));
313 init (a, get_type ());
317 template <
typename T>
326 if (! q.is_square () || u.numel () != n)
329 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
331 init (octave::math::insert_row (q*r, j, u), get_type ());
334 template <
typename T>
342 if (! q.is_square ())
344 if (j < 0 || j > m-1)
345 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
347 init (octave::math::delete_row (q*r, j), get_type ());
350 template <
typename T>
359 (*current_liboctave_error_handler) (
"qrshift: index out of range");
361 init (octave::math::shift_cols (q*r, i, j), get_type ());
382 afact.
elem (i, j) *= tau[j];
400 r.xelem (i, j) = afact.
xelem (i, j);
413 q.xelem (i, j) = afact.
xelem (i, j);
414 afact.
xelem (i, j) = 0;
424 F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (),
m, tau,
429 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
431 F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (),
m, tau,
462 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
468 form (n, afact, tau, qr_type);
471 #if defined (HAVE_QRUPDATE)
487 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
488 m, r.fortran_vec (),
k,
508 F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
509 m, r.fortran_vec (),
k,
526 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
540 F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (),
541 r.fortran_vec (), r.rows (), j + 1,
558 dups = dups && js(i) == js(i+1);
561 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
564 if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
565 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
573 r.resize (kmax, n + nj);
577 r.resize (k, n + nj);
587 r.fortran_vec (), r.rows (), js(ii) + 1,
601 if (j < 0 || j > n-1)
602 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
605 F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
606 r.fortran_vec (), r.rows (), j + 1,
w));
632 dups = dups && js(i) == js(i+1);
635 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
636 if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
645 F77_XFCN (dqrdec, DQRDEC, (m, n - ii, k == m ? k : k - ii,
646 q.fortran_vec (), q.rows (),
647 r.fortran_vec (), r.rows (),
652 q.resize (m, k - nj);
653 r.resize (k - nj, n - nj);
657 r.resize (k, n - nj);
671 if (! q.is_square () || u.
numel () != n)
674 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
676 q.resize (m + 1, m + 1);
680 F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (),
681 r.fortran_vec (), r.rows (),
693 if (! q.is_square ())
695 if (j < 0 || j > m-1)
696 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
699 F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (),
700 r.fortran_vec (), r.rows (), j + 1,
703 q.resize (m - 1, m - 1);
716 (*current_liboctave_error_handler) (
"qrshift: index out of range");
720 q.fortran_vec (), q.rows (),
721 r.fortran_vec (), r.rows (),
742 afact.
elem (i, j) *= tau[j];
760 r.xelem (i, j) = afact.
xelem (i, j);
773 q.xelem (i, j) = afact.
xelem (i, j);
774 afact.
xelem (i, j) = 0;
784 F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (),
m, tau,
789 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
791 F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (),
m, tau,
822 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
828 form (n, afact, tau, qr_type);
831 #if defined (HAVE_QRUPDATE)
847 F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
848 m, r.fortran_vec (),
k,
868 F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
869 m, r.fortran_vec (),
k,
886 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
900 F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), q.rows (),
901 r.fortran_vec (), r.rows (), j + 1,
919 dups = dups && js(i) == js(i+1);
922 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
925 if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
926 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
934 r.resize (kmax, n + nj);
938 r.resize (k, n + nj);
948 r.fortran_vec (), r.rows (), js(ii) + 1,
962 if (j < 0 || j > n-1)
963 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
966 F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
967 r.fortran_vec (), r.rows (), j + 1,
w));
993 dups = dups && js(i) == js(i+1);
996 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
997 if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
1006 F77_XFCN (sqrdec, SQRDEC, (m, n - ii, k == m ? k : k - ii,
1007 q.fortran_vec (), q.rows (),
1008 r.fortran_vec (), r.rows (),
1013 q.resize (m, k - nj);
1014 r.resize (k - nj, n - nj);
1018 r.resize (k, n - nj);
1032 if (! q.is_square () || u.
numel () != n)
1035 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1037 q.resize (m + 1, m + 1);
1038 r.resize (m + 1, n);
1041 F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), q.rows (),
1042 r.fortran_vec (), r.rows (),
1054 if (! q.is_square ())
1056 if (j < 0 || j > m-1)
1057 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1060 F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), q.rows (),
1061 r.fortran_vec (), r.rows (), j + 1,
1064 q.resize (m - 1, m - 1);
1065 r.resize (m - 1, n);
1077 (*current_liboctave_error_handler) (
"qrshift: index out of range");
1080 F77_XFCN (sqrshc, SQRSHC, (m, n, k,
1081 q.fortran_vec (), q.rows (),
1082 r.fortran_vec (), r.rows (),
1103 afact.
elem (i, j) *= tau[j];
1121 r.xelem (i, j) = afact.
xelem (i, j);
1134 q.xelem (i, j) = afact.
xelem (i, j);
1135 afact.
xelem (i, j) = 0;
1151 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
1186 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
1193 form (n, afact, tau, qr_type);
1196 #if defined (HAVE_QRUPDATE)
1257 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1262 r.resize (k+1, n+1);
1291 dups = dups && js(i) == js(i+1);
1294 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1297 if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
1298 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1306 r.resize (kmax, n + nj);
1310 r.resize (k, n + nj);
1334 if (j < 0 || j > n-1)
1335 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1345 r.resize (k-1, n-1);
1366 dups = dups && js(i) == js(i+1);
1369 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1370 if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
1379 F77_XFCN (zqrdec, ZQRDEC, (m, n - ii, k == m ? k : k - ii,
1386 q.resize (m, k - nj);
1387 r.resize (k - nj, n - nj);
1391 r.resize (k, n - nj);
1405 if (! q.is_square () || u.
numel () != n)
1408 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1410 q.resize (m + 1, m + 1);
1411 r.resize (m + 1, n);
1428 if (! q.is_square ())
1430 if (j < 0 || j > m-1)
1431 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1440 q.resize (m - 1, m - 1);
1441 r.resize (m - 1, n);
1453 (*current_liboctave_error_handler) (
"qrshift: index out of range");
1457 F77_XFCN (zqrshc, ZQRSHC, (m, n, k,
1480 afact.
elem (i, j) *= tau[j];
1498 r.xelem (i, j) = afact.
xelem (i, j);
1511 q.xelem (i, j) = afact.
xelem (i, j);
1512 afact.
xelem (i, j) = 0;
1528 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
1563 lwork =
std::max (lwork, static_cast<octave_idx_type> (1));
1570 form (n, afact, tau, qr_type);
1573 #if defined (HAVE_QRUPDATE)
1634 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1639 r.resize (k+1, n+1);
1667 dups = dups && js(i) == js(i+1);
1670 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1673 if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
1674 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1682 r.resize (kmax, n + nj);
1686 r.resize (k, n + nj);
1709 if (j < 0 || j > n-1)
1710 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1719 r.resize (k-1, n-1);
1740 dups = dups && js(i) == js(i+1);
1743 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1744 if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
1753 F77_XFCN (cqrdec, CQRDEC, (m, n - ii, k == m ? k : k - ii,
1760 q.resize (m, k - nj);
1761 r.resize (k - nj, n - nj);
1765 r.resize (k, n - nj);
1780 if (! q.is_square () || u.
numel () != n)
1783 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1785 q.resize (m + 1, m + 1);
1786 r.resize (m + 1, n);
1802 if (! q.is_square ())
1804 if (j < 0 || j > m-1)
1805 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1813 q.resize (m - 1, m - 1);
1814 r.resize (m - 1, n);
1826 (*current_liboctave_error_handler) (
"qrshift: index out of range");
1830 F77_XFCN (cqrshc, CQRSHC, (m, n, k,
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
void resize(octave_idx_type nr, octave_idx_type nc, float rfv=0)
Octave interface to the compression and uncompression libraries.
octave_idx_type rows(void) const
static octave::math::qr< T >::type qr_type(int nargin, int nargout)
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
static const idx_vector colon
octave_idx_type numel(void) const
Number of elements in the array.
OCTAVE_EXPORT octave_value_list while another program execution is suspended until the graphics object the function returns immediately In the second form
void delete_col(octave_idx_type j)
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
#define F77_DBLE_CMPLX_ARG(x)
void shift_cols(octave_idx_type i, octave_idx_type j)
void insert_row(const RV_T &u, octave_idx_type j)
T & elem(octave_idx_type n)
#define F77_XFCN(f, F, args)
void warn_qrupdate_once(void)
octave_idx_type rows(void) const
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
void insert_col(const CV_T &u, octave_idx_type j)
void update(const CV_T &u, const CV_T &v)
void init(const T &a, type qr_type)
octave_idx_type columns(void) const
FloatComplexColumnVector column(octave_idx_type i) const
Array< T > sort(int dim=0, sortmode mode=ASCENDING) const
octave_idx_type index(const T *src, octave_idx_type n, T *dest) const
nd deftypefn *octave_map m
std::complex< double > w(std::complex< double > z, double relerr=0)
const T * data(void) const
octave_value & assign(assign_op op, const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
Array< T > column(octave_idx_type k) const
Extract column: A(:,k+1).
void form(octave_idx_type n, T &afact, ELT_T *tau, type qr_type)
FloatColumnVector column(octave_idx_type i) const
void delete_row(octave_idx_type j)
T & xelem(octave_idx_type n)
type get_type(void) const
charNDArray max(char d, const charNDArray &m)
#define F77_CONST_DBLE_CMPLX_ARG(x)
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
octave_idx_type rows(void) const
=val(i)}if ode{val(i)}occurs in table i
ColumnVector column(octave_idx_type i) const
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
#define F77_CONST_CMPLX_ARG(x)
std::complex< float > FloatComplex
std::complex< double > Complex
const T * fortran_vec(void) const
octave_idx_type cols(void) const
Vector representing the dimensions (size) of an Array.
ComplexColumnVector column(octave_idx_type i) const
octave_idx_type columns(void) const
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
charNDArray min(char d, const charNDArray &m)