Mercurial > hg > octave-lyh

comparison liboctave/CmplxQR.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
91 ComplexQR::init (const ComplexMatrix& a, QR::type qr_type) 92 ComplexQR::init (const ComplexMatrix& a, QR::type qr_type)
92 { 93 {
93 octave_idx_type m = a.rows (); 94 octave_idx_type m = a.rows ();
94 octave_idx_type n = a.cols (); 95 octave_idx_type n = a.cols ();
95 96
96 if (m == 0 || n == 0)
97 {
98 (*current_liboctave_error_handler)
99 ("ComplexQR must have non-empty matrix");
100 return;
101 }
102
103 octave_idx_type min_mn = m < n ? m : n; 97 octave_idx_type min_mn = m < n ? m : n;
104 98 OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn);
105 Array<Complex> tau (min_mn);
106 Complex *ptau = tau.fortran_vec ();
107
108 octave_idx_type lwork = 32*n;
109 Array<Complex> work (lwork);
110 Complex *pwork = work.fortran_vec ();
111 99
112 octave_idx_type info = 0; 100 octave_idx_type info = 0;
113 101
114 ComplexMatrix A_fact; 102 ComplexMatrix afact = a;
115 if (m > n && qr_type != QR::economy) 103 if (m > n && qr_type == QR::std)
116 { 104 afact.resize (m, m);
117 A_fact.resize (m, m); 105
118 A_fact.insert (a, 0, 0); 106 if (m > 0)
119 } 107 {
120 else 108 // workspace query.
121 A_fact = a; 109 Complex clwork;
122 110 F77_XFCN (zgeqrf, ZGEQRF, (m, n, afact.fortran_vec (), m, tau, &clwork, -1, info));
123 Complex *tmp_data = A_fact.fortran_vec (); 111
124 112 // allocate buffer and do the job.
125 F77_XFCN (zgeqrf, ZGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info)); 113 octave_idx_type lwork = clwork.real (); lwork = std::max (lwork, 1);
114 OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
115 F77_XFCN (zgeqrf, ZGEQRF, (m, n, afact.fortran_vec (), m, tau, work, lwork, info));
116 }
117
118 form (n, afact, tau, qr_type);
119 }
120
121 void ComplexQR::form (octave_idx_type n, ComplexMatrix& afact,
122 Complex *tau, QR::type qr_type)
123 {
124 octave_idx_type m = afact.rows (), min_mn = std::min (m, n);
125 octave_idx_type info;
126 126
127 if (qr_type == QR::raw) 127 if (qr_type == QR::raw)
128 { 128 {
129 for (octave_idx_type j = 0; j < min_mn; j++) 129 for (octave_idx_type j = 0; j < min_mn; j++)
130 { 130 {
131 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; 131 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
132 for (octave_idx_type i = limit + 1; i < m; i++) 132 for (octave_idx_type i = limit + 1; i < m; i++)
133 A_fact.elem (i, j) *= tau.elem (j); 133 afact.elem (i, j) *= tau[j];
134 } 134 }
135 135
136 r = A_fact; 136 r = afact;
137 137 }
138 if (m > n) 138 else
139 r.resize (m, n); 139 {
140 } 140 // Attempt to minimize copying.
141 else 141 if (m >= n)
142 { 142 {
143 octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; 143 // afact will become q.
144 144 q = afact;
145 if (qr_type == QR::economy && m > n) 145 octave_idx_type k = qr_type == QR::economy ? n : m;
146 r.resize (n, n, 0.0); 146 r = ComplexMatrix (k, n);
147 for (octave_idx_type j = 0; j < n; j++)
148 {
149 octave_idx_type i = 0;
150 for (; i <= j; i++)
151 r.xelem (i, j) = afact.xelem (i, j);
152 for (;i < k; i++)
153 r.xelem (i, j) = 0;
154 }
155 afact = ComplexMatrix (); // optimize memory
156 }
147 else 157 else
148 r.resize (m, n, 0.0); 158 {
149 159 // afact will become r.
150 for (octave_idx_type j = 0; j < n; j++) 160 q = ComplexMatrix (m, m);
151 { 161 for (octave_idx_type j = 0; j < m; j++)
152 octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; 162 for (octave_idx_type i = j + 1; i < m; i++)
153 for (octave_idx_type i = 0; i <= limit; i++) 163 {
154 r.elem (i, j) = A_fact.elem (i, j); 164 q.xelem (i, j) = afact.xelem (i, j);
155 } 165 afact.xelem (i, j) = 0;
156 166 }
157 lwork = 32 * n2; 167 r = afact;
158 work.resize (lwork); 168 }
159 Complex *pwork2 = work.fortran_vec (); 169
160 170
161 F77_XFCN (zungqr, ZUNGQR, (m, n2, min_mn, tmp_data, m, ptau, 171 if (m > 0)
162 pwork2, lwork, info)); 172 {
163 173 octave_idx_type k = q.columns ();
164 q = A_fact; 174 // workspace query.
165 q.resize (m, n2); 175 Complex clwork;
176 F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
177 &clwork, -1, info));
178
179 // allocate buffer and do the job.
180 octave_idx_type lwork = clwork.real (); lwork = std::max (lwork, 1);
181 OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
182 F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
183 work, lwork, info));
184 }
166 } 185 }
167 } 186 }
168 187
169 ComplexQR::ComplexQR (const ComplexMatrix& q_arg, const ComplexMatrix& r_arg) 188 ComplexQR::ComplexQR (const ComplexMatrix& q_arg, const ComplexMatrix& r_arg)
170 { 189 {