Mercurial > octave-nkf

comparison liboctave/EIG.cc @ 4725:fa612b2cbfe9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2004-01-23 16:42:51 by jwe]
author jwe
date Fri, 23 Jan 2004 16:42:51 +0000
parents 6f3382e08a52
children c322edde72ac
comparison
equal deleted inserted replaced
4724:bdacd0383fbd 4725:fa612b2cbfe9
69 F77_CHAR_ARG_LEN_DECL 69 F77_CHAR_ARG_LEN_DECL
70 F77_CHAR_ARG_LEN_DECL); 70 F77_CHAR_ARG_LEN_DECL);
71 } 71 }
72 72
73 int 73 int
74 EIG::init (const Matrix& a) 74 EIG::init (const Matrix& a, bool calc_ev)
75 { 75 {
76 if (a.is_symmetric ()) 76 if (a.is_symmetric ())
77 return symmetric_init (a); 77 return symmetric_init (a, calc_ev);
78 78
79 int n = a.rows (); 79 int n = a.rows ();
80 80
81 if (n != a.cols ()) 81 if (n != a.cols ())
82 { 82 {
93 double *pwr = wr.fortran_vec (); 93 double *pwr = wr.fortran_vec ();
94 94
95 Array<double> wi (n); 95 Array<double> wi (n);
96 double *pwi = wi.fortran_vec (); 96 double *pwi = wi.fortran_vec ();
97 97
98 Matrix vr (n, n); 98 int nvr = calc_ev ? n : 0;
99 Matrix vr (nvr, nvr);
99 double *pvr = vr.fortran_vec (); 100 double *pvr = vr.fortran_vec ();
100 101
101 // XXX FIXME XXX -- it might be possible to choose a better value of 102 // XXX FIXME XXX -- it might be possible to choose a better value of
102 // lwork that would result in more efficient computations. 103 // lwork that would result in more efficient computations.
103 104
107 108
108 double *dummy = 0; 109 double *dummy = 0;
109 int idummy = 1; 110 int idummy = 1;
110 111
111 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), 112 F77_XFCN (dgeev, DGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
112 F77_CONST_CHAR_ARG2 ("V", 1), 113 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
113 n, tmp_data, n, pwr, pwi, dummy, 114 n, tmp_data, n, pwr, pwi, dummy,
114 idummy, pvr, n, pwork, lwork, info 115 idummy, pvr, n, pwork, lwork, info
115 F77_CHAR_ARG_LEN (1) 116 F77_CHAR_ARG_LEN (1)
116 F77_CHAR_ARG_LEN (1))); 117 F77_CHAR_ARG_LEN (1)));
117 118
122 if (info > 0) 123 if (info > 0)
123 (*current_liboctave_error_handler) ("dgeev failed to converge"); 124 (*current_liboctave_error_handler) ("dgeev failed to converge");
124 else 125 else
125 { 126 {
126 lambda.resize (n); 127 lambda.resize (n);
127 v.resize (n, n); 128 v.resize (nvr, nvr);
128 129
129 for (int j = 0; j < n; j++) 130 for (int j = 0; j < n; j++)
130 { 131 {
131 if (wi.elem (j) == 0.0) 132 if (wi.elem (j) == 0.0)
132 { 133 {
133 lambda.elem (j) = Complex (wr.elem (j)); 134 lambda.elem (j) = Complex (wr.elem (j));
134 for (int i = 0; i < n; i++) 135 for (int i = 0; i < nvr; i++)
135 v.elem (i, j) = vr.elem (i, j); 136 v.elem (i, j) = vr.elem (i, j);
136 } 137 }
137 else 138 else
138 { 139 {
139 if (j+1 >= n) 140 if (j+1 >= n)
144 } 145 }
145 146
146 lambda.elem(j) = Complex (wr.elem(j), wi.elem(j)); 147 lambda.elem(j) = Complex (wr.elem(j), wi.elem(j));
147 lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1)); 148 lambda.elem(j+1) = Complex (wr.elem(j+1), wi.elem(j+1));
148 149
149 for (int i = 0; i < n; i++) 150 for (int i = 0; i < nvr; i++)
150 { 151 {
151 double real_part = vr.elem (i, j); 152 double real_part = vr.elem (i, j);
152 double imag_part = vr.elem (i, j+1); 153 double imag_part = vr.elem (i, j+1);
153 v.elem (i, j) = Complex (real_part, imag_part); 154 v.elem (i, j) = Complex (real_part, imag_part);
154 v.elem (i, j+1) = Complex (real_part, -imag_part); 155 v.elem (i, j+1) = Complex (real_part, -imag_part);
161 162
162 return info; 163 return info;
163 } 164 }
164 165
165 int 166 int
166 EIG::symmetric_init (const Matrix& a) 167 EIG::symmetric_init (const Matrix& a, bool calc_ev)
167 { 168 {
168 int n = a.rows (); 169 int n = a.rows ();
169 170
170 if (n != a.cols ()) 171 if (n != a.cols ())
171 { 172 {
186 187
187 int lwork = 8*n; 188 int lwork = 8*n;
188 Array<double> work (lwork); 189 Array<double> work (lwork);
189 double *pwork = work.fortran_vec (); 190 double *pwork = work.fortran_vec ();
190 191
191 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 ("V", 1), 192 F77_XFCN (dsyev, DSYEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
192 F77_CONST_CHAR_ARG2 ("U", 1), 193 F77_CONST_CHAR_ARG2 ("U", 1),
193 n, tmp_data, n, pwr, pwork, lwork, info 194 n, tmp_data, n, pwr, pwork, lwork, info
194 F77_CHAR_ARG_LEN (1) 195 F77_CHAR_ARG_LEN (1)
195 F77_CHAR_ARG_LEN (1))); 196 F77_CHAR_ARG_LEN (1)));
196 197
199 else if (info > 0) 200 else if (info > 0)
200 (*current_liboctave_error_handler) ("dsyev failed to converge"); 201 (*current_liboctave_error_handler) ("dsyev failed to converge");
201 else 202 else
202 { 203 {
203 lambda = ComplexColumnVector (wr); 204 lambda = ComplexColumnVector (wr);
204 v = ComplexMatrix (atmp); 205 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
205 } 206 }
206 207
207 return info; 208 return info;
208 } 209 }
209 210
210 int 211 int
211 EIG::init (const ComplexMatrix& a) 212 EIG::init (const ComplexMatrix& a, bool calc_ev)
212 { 213 {
213 if (a.is_hermitian ()) 214 if (a.is_hermitian ())
214 return hermitian_init (a); 215 return hermitian_init (a, calc_ev);
215 216
216 int n = a.rows (); 217 int n = a.rows ();
217 218
218 if (n != a.cols ()) 219 if (n != a.cols ())
219 { 220 {
227 Complex *tmp_data = atmp.fortran_vec (); 228 Complex *tmp_data = atmp.fortran_vec ();
228 229
229 ComplexColumnVector w (n); 230 ComplexColumnVector w (n);
230 Complex *pw = w.fortran_vec (); 231 Complex *pw = w.fortran_vec ();
231 232
232 ComplexMatrix vtmp (n, n); 233 int nvr = calc_ev ? n : 0;
234 ComplexMatrix vtmp (nvr, nvr);
233 Complex *pv = vtmp.fortran_vec (); 235 Complex *pv = vtmp.fortran_vec ();
234 236
235 // XXX FIXME XXX -- it might be possible to choose a better value of 237 // XXX FIXME XXX -- it might be possible to choose a better value of
236 // lwork that would result in more efficient computations. 238 // lwork that would result in more efficient computations.
237 239
245 247
246 Complex *dummy = 0; 248 Complex *dummy = 0;
247 int idummy = 1; 249 int idummy = 1;
248 250
249 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1), 251 F77_XFCN (zgeev, ZGEEV, (F77_CONST_CHAR_ARG2 ("N", 1),
250 F77_CONST_CHAR_ARG2 ("V", 1), 252 F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
251 n, tmp_data, n, pw, dummy, idummy, 253 n, tmp_data, n, pw, dummy, idummy,
252 pv, n, pwork, lwork, prwork, info 254 pv, n, pwork, lwork, prwork, info
253 F77_CHAR_ARG_LEN (1) 255 F77_CHAR_ARG_LEN (1)
254 F77_CHAR_ARG_LEN (1))); 256 F77_CHAR_ARG_LEN (1)));
255 257
265 267
266 return info; 268 return info;
267 } 269 }
268 270
269 int 271 int
270 EIG::hermitian_init (const ComplexMatrix& a) 272 EIG::hermitian_init (const ComplexMatrix& a, bool calc_ev)
271 { 273 {
272 int n = a.rows (); 274 int n = a.rows ();
273 275
274 if (n != a.cols ()) 276 if (n != a.cols ())
275 { 277 {
294 296
295 int lrwork = 3*n; 297 int lrwork = 3*n;
296 Array<double> rwork (lrwork); 298 Array<double> rwork (lrwork);
297 double *prwork = rwork.fortran_vec (); 299 double *prwork = rwork.fortran_vec ();
298 300
299 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 ("V", 1), 301 F77_XFCN (zheev, ZHEEV, (F77_CONST_CHAR_ARG2 (calc_ev ? "V" : "N", 1),
300 F77_CONST_CHAR_ARG2 ("U", 1), 302 F77_CONST_CHAR_ARG2 ("U", 1),
301 n, tmp_data, n, pwr, pwork, lwork, prwork, info 303 n, tmp_data, n, pwr, pwork, lwork, prwork, info
302 F77_CHAR_ARG_LEN (1) 304 F77_CHAR_ARG_LEN (1)
303 F77_CHAR_ARG_LEN (1))); 305 F77_CHAR_ARG_LEN (1)));
304 306
307 else if (info > 0) 309 else if (info > 0)
308 (*current_liboctave_error_handler) ("zheev failed to converge"); 310 (*current_liboctave_error_handler) ("zheev failed to converge");
309 else 311 else
310 { 312 {
311 lambda = ComplexColumnVector (wr); 313 lambda = ComplexColumnVector (wr);
312 v = ComplexMatrix (atmp); 314 v = calc_ev ? ComplexMatrix (atmp) : ComplexMatrix ();
313 } 315 }
314 316
315 return info; 317 return info;
316 } 318 }
317 319