Mercurial > hg > octave-lojdl

comparison liboctave/dbleQRP.cc @ 3883:69b6bd271277

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2002-04-02 21:05:10 by jwe]
author jwe
date Tue, 02 Apr 2002 21:05:10 +0000
parents 8b262e771614
children 7da18459c08b
comparison
equal deleted inserted replaced
3882:c8c1ead8474f 3883:69b6bd271277
65 { 65 {
66 (*current_liboctave_error_handler) ("QR must have non-empty matrix"); 66 (*current_liboctave_error_handler) ("QR must have non-empty matrix");
67 return; 67 return;
68 } 68 }
69 69
70 Array<double> tau (m < n ? m : n); 70 int min_mn = m < n ? m : n;
71 Array<double> tau (min_mn);
71 double *ptau = tau.fortran_vec (); 72 double *ptau = tau.fortran_vec ();
72 73
73 int lwork = 3*n > 32*m ? 3*n : 32*m; 74 int lwork = 3*n > 32*m ? 3*n : 32*m;
74 Array<double> work (lwork); 75 Array<double> work (lwork);
75 double *pwork = work.fortran_vec (); 76 double *pwork = work.fortran_vec ();
76 77
77 int info = 0; 78 int info = 0;
78 79
79 Matrix A_fact = a; 80 Matrix A_fact = a;
80 if (m > n) 81 if (m > n && qr_type != QR::economy)
81 A_fact.resize (m, m, 0.0); 82 A_fact.resize (m, m, 0.0);
82 83
83 double *tmp_data = A_fact.fortran_vec (); 84 double *tmp_data = A_fact.fortran_vec ();
84 85
85 Array<int> jpvt (n, 0); 86 Array<int> jpvt (n, 0);
107 p.resize (n, n, 0.0); 108 p.resize (n, n, 0.0);
108 for (int j = 0; j < n; j++) 109 for (int j = 0; j < n; j++)
109 p.elem (jpvt.elem (j) - 1, j) = 1.0; 110 p.elem (jpvt.elem (j) - 1, j) = 1.0;
110 } 111 }
111 112
113 int n2 = (qr_type == QR::economy) ? min_mn : m;
114
112 if (qr_type == QR::economy && m > n) 115 if (qr_type == QR::economy && m > n)
113 r.resize (n, n, 0.0); 116 r.resize (n, n, 0.0);
114 else 117 else
115 r.resize (m, n, 0.0); 118 r.resize (m, n, 0.0);
116
117 int min_mn = m < n ? m : n;
118 119
119 for (int j = 0; j < n; j++) 120 for (int j = 0; j < n; j++)
120 { 121 {
121 int limit = j < min_mn-1 ? j : min_mn-1; 122 int limit = j < min_mn-1 ? j : min_mn-1;
122 for (int i = 0; i <= limit; i++) 123 for (int i = 0; i <= limit; i++)
123 r.elem (i, j) = A_fact.elem (i, j); 124 r.elem (i, j) = A_fact.elem (i, j);
124 } 125 }
125
126 int n2 = m;
127 if (qr_type == QR::economy)
128 n2 = min_mn;
129 126
130 F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau, 127 F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau,
131 pwork, lwork, info)); 128 pwork, lwork, info));
132 129
133 if (f77_exception_encountered) 130 if (f77_exception_encountered)