Mercurial > hg > octave-nkf

comparison liboctave/fMatrix.cc @ 8392:c187f0e3a7ee

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use m-file implementation for expm
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 10 Dec 2008 11:55:23 +0100
parents 25bc2d31e1bf
children 5114ea5a41b5
comparison
equal deleted inserted replaced
8391:343f0fbca6eb 8392:c187f0e3a7ee
2570 solve_singularity_warning (float rcon) 2570 solve_singularity_warning (float rcon)
2571 { 2571 {
2572 (*current_liboctave_warning_handler) 2572 (*current_liboctave_warning_handler)
2573 ("singular matrix encountered in expm calculation, rcond = %g", 2573 ("singular matrix encountered in expm calculation, rcond = %g",
2574 rcon); 2574 rcon);
2575 }
2576
2577 FloatMatrix
2578 FloatMatrix::expm (void) const
2579 {
2580 FloatMatrix retval;
2581
2582 FloatMatrix m = *this;
2583
2584 if (numel () == 1)
2585 return FloatMatrix (1, 1, exp (m(0)));
2586
2587 octave_idx_type nc = columns ();
2588
2589 // Preconditioning step 1: trace normalization to reduce dynamic
2590 // range of poles, but avoid making stable eigenvalues unstable.
2591
2592 // trace shift value
2593 volatile float trshift = 0.0;
2594
2595 for (octave_idx_type i = 0; i < nc; i++)
2596 trshift += m.elem (i, i);
2597
2598 trshift /= nc;
2599
2600 if (trshift > 0.0)
2601 {
2602 for (octave_idx_type i = 0; i < nc; i++)
2603 m.elem (i, i) -= trshift;
2604 }
2605
2606 // Preconditioning step 2: balancing; code follows development
2607 // in AEPBAL
2608
2609 float *p_m = m.fortran_vec ();
2610
2611 octave_idx_type info, ilo, ihi, ilos, ihis;
2612 Array<float> dpermute (nc);
2613 Array<float> dscale (nc);
2614
2615 // permutation first
2616 char job = 'P';
2617 F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
2618 nc, p_m, nc, ilo, ihi,
2619 dpermute.fortran_vec (), info
2620 F77_CHAR_ARG_LEN (1)));
2621
2622 // then scaling
2623 job = 'S';
2624 F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
2625 nc, p_m, nc, ilos, ihis,
2626 dscale.fortran_vec (), info
2627 F77_CHAR_ARG_LEN (1)));
2628
2629 // Preconditioning step 3: scaling.
2630
2631 FloatColumnVector work(nc);
2632 float inf_norm;
2633
2634 F77_XFCN (xslange, XSLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
2635 nc, nc, m.fortran_vec (), nc,
2636 work.fortran_vec (), inf_norm
2637 F77_CHAR_ARG_LEN (1)));
2638
2639 octave_idx_type sqpow = static_cast<octave_idx_type> (inf_norm > 0.0
2640 ? (1.0 + log (inf_norm) / log (2.0))
2641 : 0.0);
2642
2643 // Check whether we need to square at all.
2644
2645 if (sqpow < 0)
2646 sqpow = 0;
2647
2648 if (sqpow > 0)
2649 {
2650 if (sqpow > 1023)
2651 sqpow = 1023;
2652
2653 float scale_factor = 1.0;
2654 for (octave_idx_type i = 0; i < sqpow; i++)
2655 scale_factor *= 2.0;
2656
2657 m = m / scale_factor;
2658 }
2659
2660 // npp, dpp: pade' approx polynomial matrices.
2661
2662 FloatMatrix npp (nc, nc, 0.0);
2663 float *pnpp = npp.fortran_vec ();
2664 FloatMatrix dpp = npp;
2665 float *pdpp = dpp.fortran_vec ();
2666
2667 // Now powers a^8 ... a^1.
2668
2669 octave_idx_type minus_one_j = -1;
2670 for (octave_idx_type j = 7; j >= 0; j--)
2671 {
2672 for (octave_idx_type i = 0; i < nc; i++)
2673 {
2674 octave_idx_type k = i * nc + i;
2675 pnpp[k] += padec[j];
2676 pdpp[k] += minus_one_j * padec[j];
2677 }
2678
2679 npp = m * npp;
2680 pnpp = npp.fortran_vec ();
2681
2682 dpp = m * dpp;
2683 pdpp = dpp.fortran_vec ();
2684
2685 minus_one_j *= -1;
2686 }
2687
2688 // Zero power.
2689
2690 dpp = -dpp;
2691 for (octave_idx_type j = 0; j < nc; j++)
2692 {
2693 npp.elem (j, j) += 1.0;
2694 dpp.elem (j, j) += 1.0;
2695 }
2696
2697 // Compute pade approximation = inverse (dpp) * npp.
2698
2699 float rcon;
2700 retval = dpp.solve (npp, info, rcon, solve_singularity_warning);
2701
2702 if (info < 0)
2703 return retval;
2704
2705 // Reverse preconditioning step 3: repeated squaring.
2706
2707 while (sqpow)
2708 {
2709 retval = retval * retval;
2710 sqpow--;
2711 }
2712
2713 // Reverse preconditioning step 2: inverse balancing.
2714 // apply inverse scaling to computed exponential
2715 for (octave_idx_type i = 0; i < nc; i++)
2716 for (octave_idx_type j = 0; j < nc; j++)
2717 retval(i,j) *= dscale(i) / dscale(j);
2718
2719 OCTAVE_QUIT;
2720
2721 // construct balancing permutation vector
2722 Array<octave_idx_type> iperm (nc);
2723 for (octave_idx_type i = 0; i < nc; i++)
2724 iperm(i) = i; // identity permutation
2725
2726 // trailing permutations must be done in reverse order
2727 for (octave_idx_type i = nc - 1; i >= ihi; i--)
2728 {
2729 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
2730 octave_idx_type tmp = iperm(i);
2731 iperm(i) = iperm(swapidx);
2732 iperm(swapidx) = tmp;
2733 }
2734
2735 // leading permutations in forward order
2736 for (octave_idx_type i = 0; i < (ilo-1); i++)
2737 {
2738 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
2739 octave_idx_type tmp = iperm(i);
2740 iperm(i) = iperm (swapidx);
2741 iperm(swapidx) = tmp;
2742 }
2743
2744 // construct inverse balancing permutation vector
2745 Array<octave_idx_type> invpvec (nc);
2746 for (octave_idx_type i = 0; i < nc; i++)
2747 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method
2748
2749 OCTAVE_QUIT;
2750
2751 FloatMatrix tmpMat = retval;
2752 for (octave_idx_type i = 0; i < nc; i++)
2753 for (octave_idx_type j = 0; j < nc; j++)
2754 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
2755
2756 // Reverse preconditioning step 1: fix trace normalization.
2757
2758 if (trshift > 0.0)
2759 retval = expf (trshift) * retval;
2760
2761 return retval;
2762 } 2575 }
2763 2576
2764 FloatMatrix& 2577 FloatMatrix&
2765 FloatMatrix::operator += (const FloatDiagMatrix& a) 2578 FloatMatrix::operator += (const FloatDiagMatrix& a)
2766 { 2579 {