Mercurial > hg > octave-lyh

comparison src/xpow.cc @ 11682:b7a30502f0c9 release-3-0-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
handle possible error from EIG
author John W. Eaton <jwe@octave.org>
date Tue, 11 Mar 2008 10:49:33 -0400
parents bc4b8d973a3a
children 82be108cc558 72830070a17b
comparison
equal deleted inserted replaced
11681:114dcfc15432 11682:b7a30502f0c9
99 if (nr == 0 || nc == 0 || nr != nc) 99 if (nr == 0 || nc == 0 || nr != nc)
100 error ("for x^A, A must be square"); 100 error ("for x^A, A must be square");
101 else 101 else
102 { 102 {
103 EIG b_eig (b); 103 EIG b_eig (b);
104 ComplexColumnVector lambda (b_eig.eigenvalues ()); 104
105 ComplexMatrix Q (b_eig.eigenvectors ()); 105 if (! error_state)
106 106 {
107 for (octave_idx_type i = 0; i < nr; i++) 107 ComplexColumnVector lambda (b_eig.eigenvalues ());
108 { 108 ComplexMatrix Q (b_eig.eigenvectors ());
109 Complex elt = lambda (i); 109
110 if (std::imag (elt) == 0.0) 110 for (octave_idx_type i = 0; i < nr; i++)
111 lambda (i) = std::pow (a, std::real (elt)); 111 {
112 else 112 Complex elt = lambda(i);
113 lambda (i) = std::pow (a, elt); 113 if (std::imag (elt) == 0.0)
114 } 114 lambda(i) = std::pow (a, std::real (elt));
115 ComplexDiagMatrix D (lambda); 115 else
116 116 lambda(i) = std::pow (a, elt);
117 retval = ComplexMatrix (Q * D * Q.inverse ()); 117 }
118 ComplexDiagMatrix D (lambda);
119
120 retval = ComplexMatrix (Q * D * Q.inverse ());
121 }
122 else
123 error ("xpow: matrix diagonalization failed");
118 } 124 }
119 125
120 return retval; 126 return retval;
121 } 127 }
122 128
142 if (nr == 0 || nc == 0 || nr != nc) 148 if (nr == 0 || nc == 0 || nr != nc)
143 error ("for x^A, A must be square"); 149 error ("for x^A, A must be square");
144 else 150 else
145 { 151 {
146 EIG b_eig (b); 152 EIG b_eig (b);
147 ComplexColumnVector lambda (b_eig.eigenvalues ()); 153
148 ComplexMatrix Q (b_eig.eigenvectors ()); 154 if (! error_state)
149 155 {
150 for (octave_idx_type i = 0; i < nr; i++) 156 ComplexColumnVector lambda (b_eig.eigenvalues ());
151 { 157 ComplexMatrix Q (b_eig.eigenvectors ());
152 Complex elt = lambda (i); 158
153 if (std::imag (elt) == 0.0) 159 for (octave_idx_type i = 0; i < nr; i++)
154 lambda (i) = std::pow (a, std::real (elt)); 160 {
155 else 161 Complex elt = lambda(i);
156 lambda (i) = std::pow (a, elt); 162 if (std::imag (elt) == 0.0)
157 } 163 lambda(i) = std::pow (a, std::real (elt));
158 ComplexDiagMatrix D (lambda); 164 else
159 165 lambda(i) = std::pow (a, elt);
160 retval = ComplexMatrix (Q * D * Q.inverse ()); 166 }
167 ComplexDiagMatrix D (lambda);
168
169 retval = ComplexMatrix (Q * D * Q.inverse ());
170 }
171 else
172 error ("xpow: matrix diagonalization failed");
161 } 173 }
162 174
163 return retval; 175 return retval;
164 } 176 }
165 177
226 } 238 }
227 } 239 }
228 else 240 else
229 { 241 {
230 EIG a_eig (a); 242 EIG a_eig (a);
243
244 if (! error_state)
245 {
246 ComplexColumnVector lambda (a_eig.eigenvalues ());
247 ComplexMatrix Q (a_eig.eigenvectors ());
248
249 for (octave_idx_type i = 0; i < nr; i++)
250 lambda(i) = std::pow (lambda(i), b);
251
252 ComplexDiagMatrix D (lambda);
253
254 retval = ComplexMatrix (Q * D * Q.inverse ());
255 }
256 else
257 error ("xpow: matrix diagonalization failed");
258 }
259 }
260
261 return retval;
262 }
263
264 // -*- 6 -*-
265 octave_value
266 xpow (const Matrix& a, const Complex& b)
267 {
268 octave_value retval;
269
270 octave_idx_type nr = a.rows ();
271 octave_idx_type nc = a.cols ();
272
273 if (nr == 0 || nc == 0 || nr != nc)
274 error ("for A^b, A must be square");
275 else
276 {
277 EIG a_eig (a);
278
279 if (! error_state)
280 {
231 ComplexColumnVector lambda (a_eig.eigenvalues ()); 281 ComplexColumnVector lambda (a_eig.eigenvalues ());
232 ComplexMatrix Q (a_eig.eigenvectors ()); 282 ComplexMatrix Q (a_eig.eigenvectors ());
233 283
234 for (octave_idx_type i = 0; i < nr; i++) 284 for (octave_idx_type i = 0; i < nr; i++)
235 lambda (i) = std::pow (lambda (i), b); 285 lambda(i) = std::pow (lambda(i), b);
236 286
237 ComplexDiagMatrix D (lambda); 287 ComplexDiagMatrix D (lambda);
238 288
239 retval = ComplexMatrix (Q * D * Q.inverse ()); 289 retval = ComplexMatrix (Q * D * Q.inverse ());
240 } 290 }
241 } 291 else
242 292 error ("xpow: matrix diagonalization failed");
243 return retval;
244 }
245
246 // -*- 6 -*-
247 octave_value
248 xpow (const Matrix& a, const Complex& b)
249 {
250 octave_value retval;
251
252 octave_idx_type nr = a.rows ();
253 octave_idx_type nc = a.cols ();
254
255 if (nr == 0 || nc == 0 || nr != nc)
256 error ("for A^b, A must be square");
257 else
258 {
259 EIG a_eig (a);
260 ComplexColumnVector lambda (a_eig.eigenvalues ());
261 ComplexMatrix Q (a_eig.eigenvectors ());
262
263 for (octave_idx_type i = 0; i < nr; i++)
264 lambda (i) = std::pow (lambda (i), b);
265
266 ComplexDiagMatrix D (lambda);
267
268 retval = ComplexMatrix (Q * D * Q.inverse ());
269 } 293 }
270 294
271 return retval; 295 return retval;
272 } 296 }
273 297
297 if (nr == 0 || nc == 0 || nr != nc) 321 if (nr == 0 || nc == 0 || nr != nc)
298 error ("for x^A, A must be square"); 322 error ("for x^A, A must be square");
299 else 323 else
300 { 324 {
301 EIG b_eig (b); 325 EIG b_eig (b);
302 ComplexColumnVector lambda (b_eig.eigenvalues ()); 326
303 ComplexMatrix Q (b_eig.eigenvectors ()); 327 if (! error_state)
304 328 {
305 for (octave_idx_type i = 0; i < nr; i++) 329 ComplexColumnVector lambda (b_eig.eigenvalues ());
306 { 330 ComplexMatrix Q (b_eig.eigenvectors ());
307 Complex elt = lambda (i); 331
308 if (std::imag (elt) == 0.0) 332 for (octave_idx_type i = 0; i < nr; i++)
309 lambda (i) = std::pow (a, std::real (elt)); 333 {
310 else 334 Complex elt = lambda(i);
311 lambda (i) = std::pow (a, elt); 335 if (std::imag (elt) == 0.0)
312 } 336 lambda(i) = std::pow (a, std::real (elt));
313 ComplexDiagMatrix D (lambda); 337 else
314 338 lambda(i) = std::pow (a, elt);
315 retval = ComplexMatrix (Q * D * Q.inverse ()); 339 }
340 ComplexDiagMatrix D (lambda);
341
342 retval = ComplexMatrix (Q * D * Q.inverse ());
343 }
344 else
345 error ("xpow: matrix diagonalization failed");
316 } 346 }
317 347
318 return retval; 348 return retval;
319 } 349 }
320 350
339 if (nr == 0 || nc == 0 || nr != nc) 369 if (nr == 0 || nc == 0 || nr != nc)
340 error ("for x^A, A must be square"); 370 error ("for x^A, A must be square");
341 else 371 else
342 { 372 {
343 EIG b_eig (b); 373 EIG b_eig (b);
344 ComplexColumnVector lambda (b_eig.eigenvalues ()); 374
345 ComplexMatrix Q (b_eig.eigenvectors ()); 375 if (! error_state)
346 376 {
347 for (octave_idx_type i = 0; i < nr; i++) 377 ComplexColumnVector lambda (b_eig.eigenvalues ());
348 { 378 ComplexMatrix Q (b_eig.eigenvectors ());
349 Complex elt = lambda (i); 379
350 if (std::imag (elt) == 0.0) 380 for (octave_idx_type i = 0; i < nr; i++)
351 lambda (i) = std::pow (a, std::real (elt)); 381 {
352 else 382 Complex elt = lambda(i);
353 lambda (i) = std::pow (a, elt); 383 if (std::imag (elt) == 0.0)
354 } 384 lambda(i) = std::pow (a, std::real (elt));
355 ComplexDiagMatrix D (lambda); 385 else
356 386 lambda(i) = std::pow (a, elt);
357 retval = ComplexMatrix (Q * D * Q.inverse ()); 387 }
388 ComplexDiagMatrix D (lambda);
389
390 retval = ComplexMatrix (Q * D * Q.inverse ());
391 }
392 else
393 error ("xpow: matrix diagonalization failed");
358 } 394 }
359 395
360 return retval; 396 return retval;
361 } 397 }
362 398
423 } 459 }
424 } 460 }
425 else 461 else
426 { 462 {
427 EIG a_eig (a); 463 EIG a_eig (a);
464
465 if (! error_state)
466 {
467 ComplexColumnVector lambda (a_eig.eigenvalues ());
468 ComplexMatrix Q (a_eig.eigenvectors ());
469
470 for (octave_idx_type i = 0; i < nr; i++)
471 lambda(i) = std::pow (lambda(i), b);
472
473 ComplexDiagMatrix D (lambda);
474
475 retval = ComplexMatrix (Q * D * Q.inverse ());
476 }
477 else
478 error ("xpow: matrix diagonalization failed");
479 }
480 }
481
482 return retval;
483 }
484
485 // -*- 12 -*-
486 octave_value
487 xpow (const ComplexMatrix& a, const Complex& b)
488 {
489 octave_value retval;
490
491 octave_idx_type nr = a.rows ();
492 octave_idx_type nc = a.cols ();
493
494 if (nr == 0 || nc == 0 || nr != nc)
495 error ("for A^b, A must be square");
496 else
497 {
498 EIG a_eig (a);
499
500 if (! error_state)
501 {
428 ComplexColumnVector lambda (a_eig.eigenvalues ()); 502 ComplexColumnVector lambda (a_eig.eigenvalues ());
429 ComplexMatrix Q (a_eig.eigenvectors ()); 503 ComplexMatrix Q (a_eig.eigenvectors ());
430 504
431 for (octave_idx_type i = 0; i < nr; i++) 505 for (octave_idx_type i = 0; i < nr; i++)
432 lambda (i) = std::pow (lambda (i), b); 506 lambda(i) = std::pow (lambda(i), b);
433 507
434 ComplexDiagMatrix D (lambda); 508 ComplexDiagMatrix D (lambda);
435 509
436 retval = ComplexMatrix (Q * D * Q.inverse ()); 510 retval = ComplexMatrix (Q * D * Q.inverse ());
437 } 511 }
438 } 512 else
439 513 error ("xpow: matrix diagonalization failed");
440 return retval;
441 }
442
443 // -*- 12 -*-
444 octave_value
445 xpow (const ComplexMatrix& a, const Complex& b)
446 {
447 octave_value retval;
448
449 octave_idx_type nr = a.rows ();
450 octave_idx_type nc = a.cols ();
451
452 if (nr == 0 || nc == 0 || nr != nc)
453 error ("for A^b, A must be square");
454 else
455 {
456 EIG a_eig (a);
457 ComplexColumnVector lambda (a_eig.eigenvalues ());
458 ComplexMatrix Q (a_eig.eigenvectors ());
459
460 for (octave_idx_type i = 0; i < nr; i++)
461 lambda (i) = std::pow (lambda (i), b);
462
463 ComplexDiagMatrix D (lambda);
464
465 retval = ComplexMatrix (Q * D * Q.inverse ());
466 } 514 }
467 515
468 return retval; 516 return retval;
469 } 517 }
470 518