Mercurial > hg > octave-jordi

comparison liboctave/CMatrix.cc @ 6958:a18c784ae599

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-04 19:21:23 by dbateman]
author dbateman
date Thu, 04 Oct 2007 19:21:23 +0000
parents c05fbb1b7e1f
children 93c65f2a5668
comparison
equal deleted inserted replaced
6957:768a19157591 6958:a18c784ae599
2625 { 2625 {
2626 ComplexMatrix retval; 2626 ComplexMatrix retval;
2627 2627
2628 ComplexMatrix m = *this; 2628 ComplexMatrix m = *this;
2629 2629
2630 if (numel () == 1)
2631 return ComplexMatrix (1, 1, exp (m(0)));
2632
2633 octave_idx_type nc = columns (); 2630 octave_idx_type nc = columns ();
2634 2631
2635 // Preconditioning step 1: trace normalization to reduce dynamic 2632 // Preconditioning step 1: trace normalization to reduce dynamic
2636 // range of poles, but avoid making stable eigenvalues unstable. 2633 // range of poles, but avoid making stable eigenvalues unstable.
2637 2634
2642 trshift += m.elem (i, i); 2639 trshift += m.elem (i, i);
2643 2640
2644 trshift /= nc; 2641 trshift /= nc;
2645 2642
2646 if (trshift.real () < 0.0) 2643 if (trshift.real () < 0.0)
2647 trshift = trshift.imag (); 2644 {
2645 trshift = trshift.imag ();
2646 if (trshift.real () > 709.0)
2647 trshift = 709.0;
2648 }
2648 2649
2649 for (octave_idx_type i = 0; i < nc; i++) 2650 for (octave_idx_type i = 0; i < nc; i++)
2650 m.elem (i, i) -= trshift; 2651 m.elem (i, i) -= trshift;
2651 2652
2652 // Preconditioning step 2: eigenvalue balancing. 2653 // Preconditioning step 2: eigenvalue balancing.
2720 } 2721 }
2721 2722
2722 // npp, dpp: pade' approx polynomial matrices. 2723 // npp, dpp: pade' approx polynomial matrices.
2723 2724
2724 ComplexMatrix npp (nc, nc, 0.0); 2725 ComplexMatrix npp (nc, nc, 0.0);
2726 Complex *pnpp = npp.fortran_vec ();
2725 ComplexMatrix dpp = npp; 2727 ComplexMatrix dpp = npp;
2728 Complex *pdpp = dpp.fortran_vec ();
2726 2729
2727 // Now powers a^8 ... a^1. 2730 // Now powers a^8 ... a^1.
2728 2731
2729 int minus_one_j = -1; 2732 int minus_one_j = -1;
2730 for (octave_idx_type j = 7; j >= 0; j--) 2733 for (octave_idx_type j = 7; j >= 0; j--)
2731 { 2734 {
2732 npp = m * npp + m * padec[j]; 2735 for (octave_idx_type i = 0; i < nc; i++)
2733 dpp = m * dpp + m * (minus_one_j * padec[j]); 2736 {
2737 octave_idx_type k = i * nc + i;
2738 pnpp [k] = pnpp [k] + padec [j];
2739 pdpp [k] = pdpp [k] + minus_one_j * padec [j];
2740 }
2741 npp = m * npp;
2742 dpp = m * dpp;
2734 minus_one_j *= -1; 2743 minus_one_j *= -1;
2735 } 2744 }
2736 2745
2737 // Zero power. 2746 // Zero power.
2738 2747