Mercurial > hg > octave-lyh

comparison liboctave/dbleQR.cc @ 8597:c86718093c1b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve & fix QR classes
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 27 Jan 2009 12:40:06 +0100
parents a6edd5c23cb5
children e9cb742df9eb
comparison
equal deleted inserted replaced
8596:8833c0b18eb2 8597:c86718093c1b
1 /* 1 /*
2 2
3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007 3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
4 John W. Eaton 4 John W. Eaton
5 Copyright (C) 2008, 2009 Jaroslav Hajek 5 Copyright (C) 2008, 2009 Jaroslav Hajek
6 Copyright (C) 2009 VZLU Prague
6 7
7 This file is part of Octave. 8 This file is part of Octave.
8 9
9 Octave is free software; you can redistribute it and/or modify it 10 Octave is free software; you can redistribute it and/or modify it
10 under the terms of the GNU General Public License as published by the 11 under the terms of the GNU General Public License as published by the
89 QR::init (const Matrix& a, QR::type qr_type) 90 QR::init (const Matrix& a, QR::type qr_type)
90 { 91 {
91 octave_idx_type m = a.rows (); 92 octave_idx_type m = a.rows ();
92 octave_idx_type n = a.cols (); 93 octave_idx_type n = a.cols ();
93 94
94 if (m == 0 || n == 0)
95 {
96 (*current_liboctave_error_handler) ("QR must have non-empty matrix");
97 return;
98 }
99
100 octave_idx_type min_mn = m < n ? m : n; 95 octave_idx_type min_mn = m < n ? m : n;
101 Array<double> tau (min_mn); 96 OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
102 double *ptau = tau.fortran_vec ();
103
104 octave_idx_type lwork = 32*n;
105 Array<double> work (lwork);
106 double *pwork = work.fortran_vec ();
107 97
108 octave_idx_type info = 0; 98 octave_idx_type info = 0;
109 99
110 Matrix A_fact = a; 100 Matrix afact = a;
111 if (m > n && qr_type != QR::economy) 101 if (m > n && qr_type == QR::std)
112 A_fact.resize (m, m, 0.0); 102 afact.resize (m, m);
113 103
114 double *tmp_data = A_fact.fortran_vec (); 104 if (m > 0)
115 105 {
116 F77_XFCN (dgeqrf, DGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info)); 106 // workspace query.
107 double rlwork;
108 F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau, &rlwork, -1, info));
109
110 // allocate buffer and do the job.
111 octave_idx_type lwork = rlwork; lwork = std::max (lwork, 1);
112 OCTAVE_LOCAL_BUFFER (double, work, lwork);
113 F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau, work, lwork, info));
114 }
115
116 form (n, afact, tau, qr_type);
117 }
118
119 void QR::form (octave_idx_type n, Matrix& afact,
120 double *tau, QR::type qr_type)
121 {
122 octave_idx_type m = afact.rows (), min_mn = std::min (m, n);
123 octave_idx_type info;
117 124
118 if (qr_type == QR::raw) 125 if (qr_type == QR::raw)
119 { 126 {
120 for (octave_idx_type j = 0; j < min_mn; j++) 127 for (octave_idx_type j = 0; j < min_mn; j++)
121 { 128 {
122 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; 129 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
123 for (octave_idx_type i = limit + 1; i < m; i++) 130 for (octave_idx_type i = limit + 1; i < m; i++)
124 A_fact.elem (i, j) *= tau.elem (j); 131 afact.elem (i, j) *= tau[j];
125 } 132 }
126 133
127 r = A_fact; 134 r = afact;
128 135 }
129 if (m > n) 136 else
130 r.resize (m, n); 137 {
131 } 138 // Attempt to minimize copying.
132 else 139 if (m >= n)
133 { 140 {
134 octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; 141 // afact will become q.
135 142 q = afact;
136 if (qr_type == QR::economy && m > n) 143 octave_idx_type k = qr_type == QR::economy ? n : m;
137 r.resize (n, n, 0.0); 144 r = Matrix (k, n);
145 for (octave_idx_type j = 0; j < n; j++)
146 {
147 octave_idx_type i = 0;
148 for (; i <= j; i++)
149 r.xelem (i, j) = afact.xelem (i, j);
150 for (;i < k; i++)
151 r.xelem (i, j) = 0;
152 }
153 afact = Matrix (); // optimize memory
154 }
138 else 155 else
139 r.resize (m, n, 0.0); 156 {
140 157 // afact will become r.
141 for (octave_idx_type j = 0; j < n; j++) 158 q = Matrix (m, m);
142 { 159 for (octave_idx_type j = 0; j < m; j++)
143 octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; 160 for (octave_idx_type i = j + 1; i < m; i++)
144 for (octave_idx_type i = 0; i <= limit; i++) 161 {
145 r.elem (i, j) = tmp_data[m*j+i]; 162 q.xelem (i, j) = afact.xelem (i, j);
146 } 163 afact.xelem (i, j) = 0;
147 164 }
148 lwork = 32 * n2; 165 r = afact;
149 work.resize (lwork); 166 }
150 double *pwork2 = work.fortran_vec (); 167
151 168
152 F77_XFCN (dorgqr, DORGQR, (m, n2, min_mn, tmp_data, m, ptau, 169 if (m > 0)
153 pwork2, lwork, info)); 170 {
154 171 octave_idx_type k = q.columns ();
155 q = A_fact; 172 // workspace query.
156 q.resize (m, n2); 173 double rlwork;
174 F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
175 &rlwork, -1, info));
176
177 // allocate buffer and do the job.
178 octave_idx_type lwork = rlwork; lwork = std::max (lwork, 1);
179 OCTAVE_LOCAL_BUFFER (double, work, lwork);
180 F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
181 work, lwork, info));
182 }
157 } 183 }
158 } 184 }
159 185
160 QR::QR (const Matrix& q_arg, const Matrix& r_arg) 186 QR::QR (const Matrix& q_arg, const Matrix& r_arg)
161 { 187 {