Mercurial > octave

comparison src/tc-extras.cc @ 51:7806354a10d3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1993-08-11 20:48:00 by jwe]
author jwe
date Wed, 11 Aug 1993 20:50:40 +0000
parents 78fd87e624cb
children 7849db4b6dbc
comparison
equal deleted inserted replaced
50:6028dcac27ef 51:7806354a10d3
324 break; 324 break;
325 } 325 }
326 return retval; 326 return retval;
327 } 327 }
328 328
329 // XXX FIXME XXX -- the next three functions should really be just 329 // XXX FIXME XXX -- the next two functions (and expm) should really be just
330 // one... 330 // one...
331
332 tree_constant *
333 matrix_exp (tree_constant& a)
334 {
335 tree_constant *retval = new tree_constant [2];
336
337 tree_constant tmp = a.make_numeric ();;
338
339 if (tmp.rows () == 0 || tmp.columns () == 0)
340 {
341 int flag = user_pref.propagate_empty_matrices;
342 if (flag != 0)
343 {
344 if (flag < 0)
345 gripe_empty_arg ("expm", 0);
346 Matrix m;
347 retval = new tree_constant [2];
348 retval[0] = tree_constant (m);
349 return retval;
350 }
351 else
352 gripe_empty_arg ("expm", 1);
353 }
354
355 switch (tmp.const_type ())
356 {
357 case tree_constant_rep::matrix_constant:
358 {
359 Matrix m = tmp.matrix_value ();
360
361 int nr = m.rows ();
362 int nc = m.columns ();
363
364 if (nr == 0 || nc == 0 || nr != nc)
365 gripe_square_matrix_required ("expm");
366 else
367 {
368 EIG m_eig (m);
369 ComplexColumnVector lambda (m_eig.eigenvalues ());
370 ComplexMatrix Q (m_eig.eigenvectors ());
371
372 for (int i = 0; i < nr; i++)
373 {
374 Complex elt = lambda.elem (i);
375 if (imag (elt) == 0.0)
376 lambda.elem (i) = exp (real (elt));
377 else
378 lambda.elem (i) = exp (elt);
379 }
380
381 ComplexDiagMatrix D (lambda);
382 ComplexMatrix result = Q * D * Q.inverse ();
383
384 retval[0] = tree_constant (result);
385 }
386 }
387 break;
388 case tree_constant_rep::complex_matrix_constant:
389 {
390 ComplexMatrix m = tmp.complex_matrix_value ();
391
392 int nr = m.rows ();
393 int nc = m.columns ();
394
395 if (nr == 0 || nc == 0 || nr != nc)
396 gripe_square_matrix_required ("expm");
397 else
398 {
399 EIG m_eig (m);
400 ComplexColumnVector lambda (m_eig.eigenvalues ());
401 ComplexMatrix Q (m_eig.eigenvectors ());
402
403 for (int i = 0; i < nr; i++)
404 {
405 Complex elt = lambda.elem (i);
406 if (imag (elt) == 0.0)
407 lambda.elem (i) = exp (real (elt));
408 else
409 lambda.elem (i) = exp (elt);
410 }
411
412 ComplexDiagMatrix D (lambda);
413 ComplexMatrix result = Q * D * Q.inverse ();
414
415 retval[0] = tree_constant (result);
416 }
417 }
418 break;
419 case tree_constant_rep::scalar_constant:
420 {
421 double d = tmp.double_value ();
422 retval[0] = tree_constant (exp (d));
423 }
424 break;
425 case tree_constant_rep::complex_scalar_constant:
426 {
427 Complex c = tmp.complex_value ();
428 retval[0] = tree_constant (exp (c));
429 }
430 break;
431 default:
432 break;
433 }
434 return retval;
435 }
436 331
437 tree_constant * 332 tree_constant *
438 matrix_log (tree_constant& a) 333 matrix_log (tree_constant& a)
439 { 334 {
440 tree_constant *retval = new tree_constant [2]; 335 tree_constant *retval = new tree_constant [2];