Mercurial > hg > octave-lyh

comparison liboctave/fCmplxQR.cc @ 9713:7918eb15040c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
refactor the QR classes onto a templated base
author Jaroslav Hajek <highegg@gmail.com>
date Mon, 12 Oct 2009 10:42:49 +0200
parents b03953732530
children 4c0cdbe0acca
comparison
equal deleted inserted replaced
9712:25e1e368618c 9713:7918eb15040c
32 #include "lo-error.h" 32 #include "lo-error.h"
33 #include "Range.h" 33 #include "Range.h"
34 #include "idx-vector.h" 34 #include "idx-vector.h"
35 #include "oct-locbuf.h" 35 #include "oct-locbuf.h"
36 36
37 #include "base-qr.cc"
38
39 template class base_qr<FloatComplexMatrix>;
40
37 extern "C" 41 extern "C"
38 { 42 {
39 F77_RET_T 43 F77_RET_T
40 F77_FUNC (cgeqrf, CGEQRF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*, 44 F77_FUNC (cgeqrf, CGEQRF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*,
41 const octave_idx_type&, FloatComplex*, FloatComplex*, 45 const octave_idx_type&, FloatComplex*, FloatComplex*,
80 FloatComplex*, float*); 84 FloatComplex*, float*);
81 85
82 #endif 86 #endif
83 } 87 }
84 88
85 FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& a, QR::type qr_type) 89 FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& a, qr_type_t qr_type)
86 : q (), r ()
87 { 90 {
88 init (a, qr_type); 91 init (a, qr_type);
89 } 92 }
90 93
91 void 94 void
92 FloatComplexQR::init (const FloatComplexMatrix& a, QR::type qr_type) 95 FloatComplexQR::init (const FloatComplexMatrix& a, qr_type_t qr_type)
93 { 96 {
94 octave_idx_type m = a.rows (); 97 octave_idx_type m = a.rows ();
95 octave_idx_type n = a.cols (); 98 octave_idx_type n = a.cols ();
96 99
97 octave_idx_type min_mn = m < n ? m : n; 100 octave_idx_type min_mn = m < n ? m : n;
98 OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn); 101 OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn);
99 102
100 octave_idx_type info = 0; 103 octave_idx_type info = 0;
101 104
102 FloatComplexMatrix afact = a; 105 FloatComplexMatrix afact = a;
103 if (m > n && qr_type == QR::std) 106 if (m > n && qr_type == qr_type_std)
104 afact.resize (m, m); 107 afact.resize (m, m);
105 108
106 if (m > 0) 109 if (m > 0)
107 { 110 {
108 // workspace query. 111 // workspace query.
118 121
119 form (n, afact, tau, qr_type); 122 form (n, afact, tau, qr_type);
120 } 123 }
121 124
122 void FloatComplexQR::form (octave_idx_type n, FloatComplexMatrix& afact, 125 void FloatComplexQR::form (octave_idx_type n, FloatComplexMatrix& afact,
123 FloatComplex *tau, QR::type qr_type) 126 FloatComplex *tau, qr_type_t qr_type)
124 { 127 {
125 octave_idx_type m = afact.rows (), min_mn = std::min (m, n); 128 octave_idx_type m = afact.rows (), min_mn = std::min (m, n);
126 octave_idx_type info; 129 octave_idx_type info;
127 130
128 if (qr_type == QR::raw) 131 if (qr_type == qr_type_raw)
129 { 132 {
130 for (octave_idx_type j = 0; j < min_mn; j++) 133 for (octave_idx_type j = 0; j < min_mn; j++)
131 { 134 {
132 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; 135 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
133 for (octave_idx_type i = limit + 1; i < m; i++) 136 for (octave_idx_type i = limit + 1; i < m; i++)
141 // Attempt to minimize copying. 144 // Attempt to minimize copying.
142 if (m >= n) 145 if (m >= n)
143 { 146 {
144 // afact will become q. 147 // afact will become q.
145 q = afact; 148 q = afact;
146 octave_idx_type k = qr_type == QR::economy ? n : m; 149 octave_idx_type k = qr_type == qr_type_economy ? n : m;
147 r = FloatComplexMatrix (k, n); 150 r = FloatComplexMatrix (k, n);
148 for (octave_idx_type j = 0; j < n; j++) 151 for (octave_idx_type j = 0; j < n; j++)
149 { 152 {
150 octave_idx_type i = 0; 153 octave_idx_type i = 0;
151 for (; i <= j; i++) 154 for (; i <= j; i++)
185 work, lwork, info)); 188 work, lwork, info));
186 } 189 }
187 } 190 }
188 } 191 }
189 192
190 FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& q_arg, const FloatComplexMatrix& r_arg)
191 {
192 octave_idx_type qr = q_arg.rows (), qc = q_arg.columns ();
193 octave_idx_type rr = r_arg.rows (), rc = r_arg.columns ();
194 if (qc == rr && (qr == qc || (qr > qc && rr == rc)))
195 {
196 q = q_arg;
197 r = r_arg;
198 }
199 else
200 (*current_liboctave_error_handler) ("QR dimensions mismatch");
201 }
202
203 QR::type
204 FloatComplexQR::get_type (void) const
205 {
206 QR::type retval;
207 if (!q.is_empty () && q.is_square ())
208 retval = QR::std;
209 else if (q.rows () > q.columns () && r.is_square ())
210 retval = QR::economy;
211 else
212 retval = QR::raw;
213 return retval;
214 }
215
216 #ifdef HAVE_QRUPDATE 193 #ifdef HAVE_QRUPDATE
217 194
218 void 195 void
219 FloatComplexQR::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v) 196 FloatComplexQR::update (const FloatComplexColumnVector& u, const FloatComplexColumnVector& v)
220 { 197 {