25 #if defined (HAVE_CONFIG_H) 58 : q (q_arg), r (r_arg)
66 if (! (q_nc == r_nr && (q_nr == q_nc || (q_nr > q_nc && r_nr == r_nc))))
67 (*current_liboctave_error_handler) (
"QR dimensions mismatch");
76 if (! q.isempty () && q.issquare ())
78 else if (q.rows () > q.cols () && r.issquare ())
106 #if ! defined (HAVE_QRUPDATE) 113 static bool warned =
false;
117 (*current_liboctave_warning_with_id_handler)
118 (
"Octave:missing-dependency",
119 "In this version of Octave, QR & Cholesky updating routines " 120 "simply update the matrix and recalculate factorizations. " 121 "To use fast algorithms, link Octave with the qrupdate library. " 122 "See <http://sourceforge.net/projects/qrupdate>.");
128 template <
typename T>
137 if (
u.numel () != m || v.numel () != n)
140 init (q*r + T (
u) * T (v).hermitian (), get_type ());
143 template <
typename T>
152 if (
u.rows () != m || v.rows () != n ||
u.cols () != v.cols ())
155 init (q*r +
u * v.hermitian (), get_type ());
158 template <
typename T,
typename CV_T>
163 T
retval (
a.rows (),
a.cols () + 1);
172 template <
typename T,
typename RV_T>
177 T
retval (
a.rows () + 1,
a.cols ());
186 template <
typename T>
196 template <
typename T>
206 template <
typename T>
228 template <
typename T>
241 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
243 init (math::insert_col (q*r, j,
u), get_type ());
246 template <
typename T>
260 dups = dups && js(
i) == js(
i+1);
263 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
265 if (
u.numel () != m ||
u.cols () != nj)
268 if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
269 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
275 a = math::insert_col (
a, js(
i),
u.column (
i));
277 init (
a, get_type ());
281 template <
typename T>
289 if (j < 0 || j > n-1)
290 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
292 init (math::delete_col (q*r, j), get_type ());
295 template <
typename T>
308 dups = dups && js(
i) == js(
i+1);
311 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
313 if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
320 a = math::delete_col (
a, js(
i));
322 init (
a, get_type ());
326 template <
typename T>
335 if (! q.issquare () ||
u.numel () != n)
339 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
341 init (math::insert_row (q*r, j,
u), get_type ());
344 template <
typename T>
355 if (j < 0 || j > m-1)
356 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
358 init (math::delete_row (q*r, j), get_type ());
361 template <
typename T>
369 if (i < 0 || i > n-1 || j < 0 || j > n-1)
370 (*current_liboctave_error_handler) (
"qrshift: index out of range");
372 init (math::shift_cols (q*r,
i, j), get_type ());
384 F77_INT n = to_f77_int (n_arg);
391 for (
F77_INT j = 0; j < min_mn; j++)
393 F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
395 afact.
elem (
i, j) *= tau[j];
409 for (
F77_INT j = 0; j < n; j++)
413 r.xelem (
i, j) = afact.
xelem (
i, j);
423 for (
F77_INT j = 0; j < m; j++)
426 q.xelem (
i, j) = afact.
xelem (
i, j);
437 F77_XFCN (dorgqr, DORGQR, (m,
k, min_mn, q.fortran_vec (), m,
438 tau, &rlwork, -1, info));
442 lwork =
std::max (lwork, static_cast<F77_INT> (1));
444 F77_XFCN (dorgqr, DORGQR, (m,
k, min_mn, q.fortran_vec (), m,
445 tau, work, lwork, info));
457 F77_INT min_mn = (m < n ? m : n);
475 lwork =
std::max (lwork, static_cast<F77_INT> (1));
484 #if defined (HAVE_QRUPDATE) 490 F77_INT m = to_f77_int (q.rows ());
491 F77_INT n = to_f77_int (r.cols ());
494 F77_INT u_nel = to_f77_int (
u.numel ());
497 if (u_nel != m || v_nel != n)
498 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
503 F77_XFCN (dqr1up, DQR1UP, (m, n,
k, q.fortran_vec (), m,
504 r.fortran_vec (),
k, utmp.fortran_vec (),
512 F77_INT m = to_f77_int (q.rows ());
513 F77_INT n = to_f77_int (r.cols ());
516 F77_INT u_rows = to_f77_int (
u.rows ());
517 F77_INT u_cols = to_f77_int (
u.cols ());
522 if (u_rows != m || v_rows != n || u_cols != v_cols)
523 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
530 F77_XFCN (dqr1up, DQR1UP, (m, n,
k, q.fortran_vec (), m,
540 F77_INT j = to_f77_int (j_arg);
542 F77_INT m = to_f77_int (q.rows ());
543 F77_INT n = to_f77_int (r.cols ());
546 F77_INT u_nel = to_f77_int (
u.numel ());
549 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
552 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
562 F77_INT ldq = to_f77_int (q.rows ());
563 F77_INT ldr = to_f77_int (r.rows ());
567 F77_XFCN (dqrinc, DQRINC, (m, n,
k, q.fortran_vec (), ldq,
568 r.fortran_vec (), ldr, j + 1,
576 F77_INT m = to_f77_int (q.rows ());
577 F77_INT n = to_f77_int (r.cols ());
585 dups = dups && js(
i) == js(
i+1);
588 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
590 F77_INT u_nel = to_f77_int (
u.numel ());
591 F77_INT u_cols = to_f77_int (
u.cols ());
593 if (u_nel != m || u_cols != nj)
594 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
596 F77_INT js_beg = to_f77_int (js(0));
597 F77_INT js_end = to_f77_int (js(nj-1));
599 if (nj > 0 && (js_beg < 0 || js_end > n))
600 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
608 r.resize (kmax, n + nj);
611 r.resize (
k, n + nj);
613 F77_INT ldq = to_f77_int (q.rows ());
614 F77_INT ldr = to_f77_int (r.rows ());
621 F77_INT js_elt = to_f77_int (js(ii));
624 r.fortran_vec (), ldr, js_elt + 1,
634 F77_INT j = to_f77_int (j_arg);
636 F77_INT m = to_f77_int (q.rows ());
638 F77_INT n = to_f77_int (r.cols ());
640 if (j < 0 || j > n-1)
641 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
643 F77_INT ldq = to_f77_int (q.rows ());
644 F77_INT ldr = to_f77_int (r.rows ());
647 F77_XFCN (dqrdec, DQRDEC, (m, n,
k, q.fortran_vec (), ldq,
648 r.fortran_vec (), ldr, j + 1,
w));
663 F77_INT m = to_f77_int (q.rows ());
664 F77_INT n = to_f77_int (r.cols ());
672 dups = dups && js(
i) == js(
i+1);
675 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
677 F77_INT js_beg = to_f77_int (js(0));
678 F77_INT js_end = to_f77_int (js(nj-1));
680 if (nj > 0 && (js_beg > n-1 || js_end < 0))
681 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
685 F77_INT ldq = to_f77_int (q.rows ());
686 F77_INT ldr = to_f77_int (r.rows ());
692 F77_INT js_elt = to_f77_int (js(ii));
693 F77_XFCN (dqrdec, DQRDEC, (m, n - ii, (
k == m ?
k :
k - ii),
694 q.fortran_vec (), ldq,
695 r.fortran_vec (), ldr,
701 q.resize (m,
k - nj);
702 r.resize (
k - nj, n - nj);
705 r.resize (
k, n - nj);
713 F77_INT j = to_f77_int (j_arg);
715 F77_INT m = to_f77_int (r.rows ());
716 F77_INT n = to_f77_int (r.cols ());
719 F77_INT u_nel = to_f77_int (
u.numel ());
721 if (! q.issquare () || u_nel != n)
725 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
727 q.resize (m + 1, m + 1);
730 F77_INT ldq = to_f77_int (q.rows ());
731 F77_INT ldr = to_f77_int (r.rows ());
735 F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), ldq,
736 r.fortran_vec (), ldr,
737 j + 1, utmp.fortran_vec (),
w));
745 F77_INT j = to_f77_int (j_arg);
747 F77_INT m = to_f77_int (r.rows ());
748 F77_INT n = to_f77_int (r.cols ());
753 if (j < 0 || j > m-1)
754 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
756 F77_INT ldq = to_f77_int (q.rows ());
757 F77_INT ldr = to_f77_int (r.rows ());
760 F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), ldq,
761 r.fortran_vec (), ldr, j + 1,
w));
763 q.resize (m - 1, m - 1);
772 F77_INT j = to_f77_int (j_arg);
774 F77_INT m = to_f77_int (q.rows ());
776 F77_INT n = to_f77_int (r.cols ());
778 if (i < 0 || i > n-1 || j < 0 || j > n-1)
779 (*current_liboctave_error_handler) (
"qrshift: index out of range");
781 F77_INT ldq = to_f77_int (q.rows ());
782 F77_INT ldr = to_f77_int (r.rows ());
786 q.fortran_vec (), ldq,
787 r.fortran_vec (), ldr,
798 F77_INT n = to_f77_int (n_arg);
805 for (
F77_INT j = 0; j < min_mn; j++)
807 F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
809 afact.
elem (
i, j) *= tau[j];
823 for (
F77_INT j = 0; j < n; j++)
827 r.xelem (
i, j) = afact.
xelem (
i, j);
837 for (
F77_INT j = 0; j < m; j++)
840 q.xelem (
i, j) = afact.
xelem (
i, j);
851 F77_XFCN (sorgqr, SORGQR, (m,
k, min_mn, q.fortran_vec (), m,
852 tau, &rlwork, -1, info));
856 lwork =
std::max (lwork, static_cast<F77_INT> (1));
858 F77_XFCN (sorgqr, SORGQR, (m,
k, min_mn, q.fortran_vec (), m,
859 tau, work, lwork, info));
871 F77_INT min_mn = (m < n ? m : n);
889 lwork =
std::max (lwork, static_cast<F77_INT> (1));
898 #if defined (HAVE_QRUPDATE) 904 F77_INT m = to_f77_int (q.rows ());
905 F77_INT n = to_f77_int (r.cols ());
908 F77_INT u_nel = to_f77_int (
u.numel ());
911 if (u_nel != m || v_nel != n)
912 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
917 F77_XFCN (sqr1up, SQR1UP, (m, n,
k, q.fortran_vec (), m,
918 r.fortran_vec (),
k, utmp.fortran_vec (),
926 F77_INT m = to_f77_int (q.rows ());
927 F77_INT n = to_f77_int (r.cols ());
930 F77_INT u_rows = to_f77_int (
u.rows ());
931 F77_INT u_cols = to_f77_int (
u.cols ());
936 if (u_rows != m || v_rows != n || u_cols != v_cols)
937 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
944 F77_XFCN (sqr1up, SQR1UP, (m, n,
k, q.fortran_vec (), m,
955 F77_INT j = to_f77_int (j_arg);
957 F77_INT m = to_f77_int (q.rows ());
958 F77_INT n = to_f77_int (r.cols ());
961 F77_INT u_nel = to_f77_int (
u.numel ());
964 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
967 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
977 F77_INT ldq = to_f77_int (q.rows ());
978 F77_INT ldr = to_f77_int (r.rows ());
982 F77_XFCN (sqrinc, SQRINC, (m, n,
k, q.fortran_vec (), ldq,
983 r.fortran_vec (), ldr, j + 1,
992 F77_INT m = to_f77_int (q.rows ());
993 F77_INT n = to_f77_int (r.cols ());
1001 dups = dups && js(
i) == js(
i+1);
1004 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1006 F77_INT u_nel = to_f77_int (
u.numel ());
1007 F77_INT u_cols = to_f77_int (
u.cols ());
1009 if (u_nel != m || u_cols != nj)
1010 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
1012 F77_INT js_beg = to_f77_int (js(0));
1013 F77_INT js_end = to_f77_int (js(nj-1));
1015 if (nj > 0 && (js_beg < 0 || js_end > n))
1016 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1024 r.resize (kmax, n + nj);
1027 r.resize (
k, n + nj);
1029 F77_INT ldq = to_f77_int (q.rows ());
1030 F77_INT ldr = to_f77_int (r.rows ());
1037 F77_INT js_elt = to_f77_int (js(ii));
1040 r.fortran_vec (), ldr, js_elt + 1,
1050 F77_INT j = to_f77_int (j_arg);
1052 F77_INT m = to_f77_int (q.rows ());
1053 F77_INT k = to_f77_int (r.rows ());
1054 F77_INT n = to_f77_int (r.cols ());
1056 if (j < 0 || j > n-1)
1057 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1059 F77_INT ldq = to_f77_int (q.rows ());
1060 F77_INT ldr = to_f77_int (r.rows ());
1063 F77_XFCN (sqrdec, SQRDEC, (m, n,
k, q.fortran_vec (), ldq,
1064 r.fortran_vec (), ldr, j + 1,
w));
1069 r.resize (
k-1, n-1);
1079 F77_INT m = to_f77_int (q.rows ());
1080 F77_INT n = to_f77_int (r.cols ());
1081 F77_INT k = to_f77_int (q.cols ());
1088 dups = dups && js(
i) == js(
i+1);
1091 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1093 F77_INT js_beg = to_f77_int (js(0));
1094 F77_INT js_end = to_f77_int (js(nj-1));
1096 if (nj > 0 && (js_beg > n-1 || js_end < 0))
1097 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1101 F77_INT ldq = to_f77_int (q.rows ());
1102 F77_INT ldr = to_f77_int (r.rows ());
1108 F77_INT js_elt = to_f77_int (js(ii));
1109 F77_XFCN (sqrdec, SQRDEC, (m, n - ii, (
k == m ?
k :
k - ii),
1110 q.fortran_vec (), ldq,
1111 r.fortran_vec (), ldr,
1117 q.resize (m,
k - nj);
1118 r.resize (
k - nj, n - nj);
1121 r.resize (
k, n - nj);
1130 F77_INT j = to_f77_int (j_arg);
1132 F77_INT m = to_f77_int (r.rows ());
1133 F77_INT n = to_f77_int (r.cols ());
1136 F77_INT u_nel = to_f77_int (
u.numel ());
1138 if (! q.issquare () || u_nel != n)
1142 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1144 q.resize (m + 1, m + 1);
1145 r.resize (m + 1, n);
1147 F77_INT ldq = to_f77_int (q.rows ());
1148 F77_INT ldr = to_f77_int (r.rows ());
1152 F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), ldq,
1153 r.fortran_vec (), ldr,
1154 j + 1, utmp.fortran_vec (),
w));
1162 F77_INT j = to_f77_int (j_arg);
1164 F77_INT m = to_f77_int (r.rows ());
1165 F77_INT n = to_f77_int (r.cols ());
1167 if (! q.issquare ())
1170 if (j < 0 || j > m-1)
1171 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1173 F77_INT ldq = to_f77_int (q.rows ());
1174 F77_INT ldr = to_f77_int (r.rows ());
1177 F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), ldq,
1178 r.fortran_vec (), ldr, j + 1,
1181 q.resize (m - 1, m - 1);
1182 r.resize (m - 1, n);
1190 F77_INT j = to_f77_int (j_arg);
1192 F77_INT m = to_f77_int (q.rows ());
1193 F77_INT k = to_f77_int (r.rows ());
1194 F77_INT n = to_f77_int (r.cols ());
1196 if (i < 0 || i > n-1 || j < 0 || j > n-1)
1197 (*current_liboctave_error_handler) (
"qrshift: index out of range");
1199 F77_INT ldq = to_f77_int (q.rows ());
1200 F77_INT ldr = to_f77_int (r.rows ());
1204 q.fortran_vec (), ldq,
1205 r.fortran_vec (), ldr,
1216 F77_INT n = to_f77_int (n_arg);
1223 for (
F77_INT j = 0; j < min_mn; j++)
1225 F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
1227 afact.
elem (
i, j) *= tau[j];
1241 for (
F77_INT j = 0; j < n; j++)
1245 r.xelem (
i, j) = afact.
xelem (
i, j);
1255 for (
F77_INT j = 0; j < m; j++)
1258 q.xelem (
i, j) = afact.
xelem (
i, j);
1266 F77_INT k = to_f77_int (q.cols ());
1269 F77_XFCN (zungqr, ZUNGQR, (m,
k, min_mn,
1277 lwork =
std::max (lwork, static_cast<F77_INT> (1));
1279 F77_XFCN (zungqr, ZUNGQR, (m,
k, min_mn,
1292 F77_INT m = to_f77_int (
a.rows ());
1293 F77_INT n = to_f77_int (
a.cols ());
1295 F77_INT min_mn = (m < n ? m : n);
1315 lwork =
std::max (lwork, static_cast<F77_INT> (1));
1326 #if defined (HAVE_QRUPDATE) 1333 F77_INT m = to_f77_int (q.rows ());
1334 F77_INT n = to_f77_int (r.cols ());
1335 F77_INT k = to_f77_int (q.cols ());
1337 F77_INT u_nel = to_f77_int (
u.numel ());
1340 if (u_nel != m || v_nel != n)
1341 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
1358 F77_INT m = to_f77_int (q.rows ());
1359 F77_INT n = to_f77_int (r.cols ());
1360 F77_INT k = to_f77_int (q.cols ());
1362 F77_INT u_rows = to_f77_int (
u.rows ());
1363 F77_INT u_cols = to_f77_int (
u.cols ());
1368 if (u_rows != m || v_rows != n || u_cols != v_cols)
1369 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
1373 for (
volatile F77_INT i = 0;
i < u_cols;
i++)
1391 F77_INT j = to_f77_int (j_arg);
1393 F77_INT m = to_f77_int (q.rows ());
1394 F77_INT n = to_f77_int (r.cols ());
1395 F77_INT k = to_f77_int (q.cols ());
1397 F77_INT u_nel = to_f77_int (
u.numel ());
1400 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
1403 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1408 r.resize (
k+1, n+1);
1413 F77_INT ldq = to_f77_int (q.rows ());
1414 F77_INT ldr = to_f77_int (r.rows ());
1429 F77_INT m = to_f77_int (q.rows ());
1430 F77_INT n = to_f77_int (r.cols ());
1431 F77_INT k = to_f77_int (q.cols ());
1438 dups = dups && js(
i) == js(
i+1);
1441 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1443 F77_INT u_nel = to_f77_int (
u.numel ());
1444 F77_INT u_cols = to_f77_int (
u.cols ());
1446 if (u_nel != m || u_cols != nj)
1447 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
1449 F77_INT js_beg = to_f77_int (js(0));
1450 F77_INT js_end = to_f77_int (js(nj-1));
1452 if (nj > 0 && (js_beg < 0 || js_end > n))
1453 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1461 r.resize (kmax, n + nj);
1464 r.resize (
k, n + nj);
1466 F77_INT ldq = to_f77_int (q.rows ());
1467 F77_INT ldr = to_f77_int (r.rows ());
1474 F77_INT js_elt = to_f77_int (js(ii));
1490 F77_INT j = to_f77_int (j_arg);
1492 F77_INT m = to_f77_int (q.rows ());
1493 F77_INT k = to_f77_int (r.rows ());
1494 F77_INT n = to_f77_int (r.cols ());
1496 if (j < 0 || j > n-1)
1497 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1499 F77_INT ldq = to_f77_int (q.rows ());
1500 F77_INT ldr = to_f77_int (r.rows ());
1510 r.resize (
k-1, n-1);
1520 F77_INT m = to_f77_int (q.rows ());
1521 F77_INT n = to_f77_int (r.cols ());
1522 F77_INT k = to_f77_int (q.cols ());
1529 dups = dups && js(
i) == js(
i+1);
1532 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1534 F77_INT js_beg = to_f77_int (js(0));
1535 F77_INT js_end = to_f77_int (js(nj-1));
1537 if (nj > 0 && (js_beg > n-1 || js_end < 0))
1538 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1542 F77_INT ldq = to_f77_int (q.rows ());
1543 F77_INT ldr = to_f77_int (r.rows ());
1549 F77_INT js_elt = to_f77_int (js(ii));
1550 F77_XFCN (zqrdec, ZQRDEC, (m, n - ii, (
k == m ?
k :
k - ii),
1554 ldr, js_elt + 1, rw));
1559 q.resize (m,
k - nj);
1560 r.resize (
k - nj, n - nj);
1563 r.resize (
k, n - nj);
1572 F77_INT j = to_f77_int (j_arg);
1574 F77_INT m = to_f77_int (r.rows ());
1575 F77_INT n = to_f77_int (r.cols ());
1578 F77_INT u_nel = to_f77_int (
u.numel ());
1580 if (! q.issquare () || u_nel != n)
1584 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1586 q.resize (m + 1, m + 1);
1587 r.resize (m + 1, n);
1589 F77_INT ldq = to_f77_int (q.rows ());
1590 F77_INT ldr = to_f77_int (r.rows ());
1605 F77_INT j = to_f77_int (j_arg);
1607 F77_INT m = to_f77_int (r.rows ());
1608 F77_INT n = to_f77_int (r.cols ());
1610 if (! q.issquare ())
1613 if (j < 0 || j > m-1)
1614 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1616 F77_INT ldq = to_f77_int (q.rows ());
1617 F77_INT ldr = to_f77_int (r.rows ());
1625 q.resize (m - 1, m - 1);
1626 r.resize (m - 1, n);
1635 F77_INT j = to_f77_int (j_arg);
1637 F77_INT m = to_f77_int (q.rows ());
1638 F77_INT k = to_f77_int (r.rows ());
1639 F77_INT n = to_f77_int (r.cols ());
1641 if (i < 0 || i > n-1 || j < 0 || j > n-1)
1642 (*current_liboctave_error_handler) (
"qrshift: index out of range");
1644 F77_INT ldq = to_f77_int (q.rows ());
1645 F77_INT ldr = to_f77_int (r.rows ());
1662 F77_INT n = to_f77_int (n_arg);
1669 for (
F77_INT j = 0; j < min_mn; j++)
1671 F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
1673 afact.
elem (
i, j) *= tau[j];
1687 for (
F77_INT j = 0; j < n; j++)
1691 r.xelem (
i, j) = afact.
xelem (
i, j);
1701 for (
F77_INT j = 0; j < m; j++)
1704 q.xelem (
i, j) = afact.
xelem (
i, j);
1712 F77_INT k = to_f77_int (q.cols ());
1715 F77_XFCN (cungqr, CUNGQR, (m,
k, min_mn,
1722 lwork =
std::max (lwork, static_cast<F77_INT> (1));
1724 F77_XFCN (cungqr, CUNGQR, (m,
k, min_mn,
1736 F77_INT m = to_f77_int (
a.rows ());
1737 F77_INT n = to_f77_int (
a.cols ());
1739 F77_INT min_mn = (m < n ? m : n);
1758 lwork =
std::max (lwork, static_cast<F77_INT> (1));
1768 #if defined (HAVE_QRUPDATE) 1775 F77_INT m = to_f77_int (q.rows ());
1776 F77_INT n = to_f77_int (r.cols ());
1777 F77_INT k = to_f77_int (q.cols ());
1779 F77_INT u_nel = to_f77_int (
u.numel ());
1782 if (u_nel != m || v_nel != n)
1783 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
1801 F77_INT m = to_f77_int (q.rows ());
1802 F77_INT n = to_f77_int (r.cols ());
1803 F77_INT k = to_f77_int (q.cols ());
1805 F77_INT u_rows = to_f77_int (
u.rows ());
1806 F77_INT u_cols = to_f77_int (
u.cols ());
1811 if (u_rows != m || v_rows != n || u_cols != v_cols)
1812 (*current_liboctave_error_handler) (
"qrupdate: dimensions mismatch");
1816 for (
volatile F77_INT i = 0;
i < u_cols;
i++)
1833 F77_INT j = to_f77_int (j_arg);
1835 F77_INT m = to_f77_int (q.rows ());
1836 F77_INT n = to_f77_int (r.cols ());
1837 F77_INT k = to_f77_int (q.cols ());
1839 F77_INT u_nel = to_f77_int (
u.numel ());
1842 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
1845 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1850 r.resize (
k+1, n+1);
1855 F77_INT ldq = to_f77_int (q.rows ());
1856 F77_INT ldr = to_f77_int (r.rows ());
1870 F77_INT m = to_f77_int (q.rows ());
1871 F77_INT n = to_f77_int (r.cols ());
1872 F77_INT k = to_f77_int (q.cols ());
1879 dups = dups && js(
i) == js(
i+1);
1882 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1884 F77_INT u_nel = to_f77_int (
u.numel ());
1885 F77_INT u_cols = to_f77_int (
u.cols ());
1887 if (u_nel != m || u_cols != nj)
1888 (*current_liboctave_error_handler) (
"qrinsert: dimensions mismatch");
1890 F77_INT js_beg = to_f77_int (js(0));
1891 F77_INT js_end = to_f77_int (js(nj-1));
1893 if (nj > 0 && (js_beg < 0 || js_end > n))
1894 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1902 r.resize (kmax, n + nj);
1905 r.resize (
k, n + nj);
1907 F77_INT ldq = to_f77_int (q.rows ());
1908 F77_INT ldr = to_f77_int (r.rows ());
1914 F77_INT js_elt = to_f77_int (js(ii));
1929 F77_INT j = to_f77_int (j_arg);
1931 F77_INT m = to_f77_int (q.rows ());
1932 F77_INT k = to_f77_int (r.rows ());
1933 F77_INT n = to_f77_int (r.cols ());
1935 if (j < 0 || j > n-1)
1936 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
1938 F77_INT ldq = to_f77_int (q.rows ());
1939 F77_INT ldr = to_f77_int (r.rows ());
1949 r.resize (
k-1, n-1);
1959 F77_INT m = to_f77_int (q.rows ());
1960 F77_INT n = to_f77_int (r.cols ());
1961 F77_INT k = to_f77_int (q.cols ());
1968 dups = dups && js(
i) == js(
i+1);
1971 (*current_liboctave_error_handler) (
"qrinsert: duplicate index detected");
1973 F77_INT js_beg = to_f77_int (js(0));
1974 F77_INT js_end = to_f77_int (js(nj-1));
1976 if (nj > 0 && (js_beg > n-1 || js_end < 0))
1977 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
1981 F77_INT ldq = to_f77_int (q.rows ());
1982 F77_INT ldr = to_f77_int (r.rows ());
1988 F77_INT js_elt = to_f77_int (js(ii));
1989 F77_XFCN (cqrdec, CQRDEC, (m, n - ii, (
k == m ?
k :
k - ii),
1997 q.resize (m,
k - nj);
1998 r.resize (
k - nj, n - nj);
2001 r.resize (
k, n - nj);
2010 F77_INT j = to_f77_int (j_arg);
2012 F77_INT m = to_f77_int (r.rows ());
2013 F77_INT n = to_f77_int (r.cols ());
2016 F77_INT u_nel = to_f77_int (
u.numel ());
2018 if (! q.issquare () || u_nel != n)
2022 (*current_liboctave_error_handler) (
"qrinsert: index out of range");
2024 q.resize (m + 1, m + 1);
2025 r.resize (m + 1, n);
2027 F77_INT ldq = to_f77_int (q.rows ());
2028 F77_INT ldr = to_f77_int (r.rows ());
2043 F77_INT j = to_f77_int (j_arg);
2045 F77_INT m = to_f77_int (r.rows ());
2046 F77_INT n = to_f77_int (r.cols ());
2048 if (! q.issquare ())
2051 if (j < 0 || j > m-1)
2052 (*current_liboctave_error_handler) (
"qrdelete: index out of range");
2054 F77_INT ldq = to_f77_int (q.rows ());
2055 F77_INT ldr = to_f77_int (r.rows ());
2063 q.resize (m - 1, m - 1);
2064 r.resize (m - 1, n);
2073 F77_INT j = to_f77_int (j_arg);
2075 F77_INT m = to_f77_int (q.rows ());
2076 F77_INT k = to_f77_int (r.rows ());
2077 F77_INT n = to_f77_int (r.cols ());
2079 if (i < 0 || i > n-1 || j < 0 || j > n-1)
2080 (*current_liboctave_error_handler) (
"qrshift: index out of range");
2082 F77_INT ldq = to_f77_int (q.rows ());
2083 F77_INT ldr = to_f77_int (r.rows ());
octave_idx_type rows(void) const
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)
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
const T * data(void) const
static const idx_vector colon
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)
const T * fortran_vec(void) const
type get_type(void) const
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 cols(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)
Array< T > sort(int dim=0, sortmode mode=ASCENDING) const
void update(const CV_T &u, const CV_T &v)
void init(const T &a, type qr_type)
octave_value & assign(assign_op op, const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
static octave::math::qr< T >::type qr_type(int nargout, bool economy)
std::complex< double > w(std::complex< double > z, double relerr=0)
octave_idx_type rows(void) const
ColumnVector column(octave_idx_type i) const
void form(octave_idx_type n, T &afact, ELT_T *tau, type qr_type)
void delete_row(octave_idx_type j)
T & xelem(octave_idx_type n)
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_f77_int_type F77_INT
FloatComplexColumnVector column(octave_idx_type i) const
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
#define F77_CONST_CMPLX_ARG(x)
ComplexColumnVector column(octave_idx_type i) const
std::complex< float > FloatComplex
std::complex< double > Complex
FloatColumnVector column(octave_idx_type i) const
octave_idx_type numel(void) const
Number of elements in the array.
Vector representing the dimensions (size) of an Array.
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 const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
charNDArray min(char d, const charNDArray &m)