Mercurial > hg > octave-lyh

comparison liboctave/fCMatrix.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
2916 solve_singularity_warning (float rcon) 2916 solve_singularity_warning (float rcon)
2917 { 2917 {
2918 (*current_liboctave_warning_handler) 2918 (*current_liboctave_warning_handler)
2919 ("singular matrix encountered in expm calculation, rcond = %g", 2919 ("singular matrix encountered in expm calculation, rcond = %g",
2920 rcon); 2920 rcon);
2921 }
2922
2923 FloatComplexMatrix
2924 FloatComplexMatrix::expm (void) const
2925 {
2926 FloatComplexMatrix retval;
2927
2928 FloatComplexMatrix m = *this;
2929
2930 octave_idx_type nc = columns ();
2931
2932 // Preconditioning step 1: trace normalization to reduce dynamic
2933 // range of poles, but avoid making stable eigenvalues unstable.
2934
2935 // trace shift value
2936 FloatComplex trshift = 0.0;
2937
2938 for (octave_idx_type i = 0; i < nc; i++)
2939 trshift += m.elem (i, i);
2940
2941 trshift /= nc;
2942
2943 if (trshift.real () < 0.0)
2944 {
2945 trshift = trshift.imag ();
2946 if (trshift.real () > 709.0)
2947 trshift = 709.0;
2948 }
2949
2950 for (octave_idx_type i = 0; i < nc; i++)
2951 m.elem (i, i) -= trshift;
2952
2953 // Preconditioning step 2: eigenvalue balancing.
2954 // code follows development in AEPBAL
2955
2956 FloatComplex *mp = m.fortran_vec ();
2957
2958 octave_idx_type info, ilo, ihi,ilos,ihis;
2959 Array<float> dpermute (nc);
2960 Array<float> dscale (nc);
2961
2962 // FIXME -- should pass job as a parameter in expm
2963
2964 // Permute first
2965 char job = 'P';
2966 F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
2967 nc, mp, nc, ilo, ihi,
2968 dpermute.fortran_vec (), info
2969 F77_CHAR_ARG_LEN (1)));
2970
2971 // then scale
2972 job = 'S';
2973 F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
2974 nc, mp, nc, ilos, ihis,
2975 dscale.fortran_vec (), info
2976 F77_CHAR_ARG_LEN (1)));
2977
2978 // Preconditioning step 3: scaling.
2979
2980 FloatColumnVector work (nc);
2981 float inf_norm;
2982
2983 F77_XFCN (xclange, XCLANGE, (F77_CONST_CHAR_ARG2 ("I", 1),
2984 nc, nc, m.fortran_vec (), nc,
2985 work.fortran_vec (), inf_norm
2986 F77_CHAR_ARG_LEN (1)));
2987
2988 int sqpow = (inf_norm > 0.0
2989 ? static_cast<int> (1.0 + log (inf_norm) / log (2.0)) : 0);
2990
2991 // Check whether we need to square at all.
2992
2993 if (sqpow < 0)
2994 sqpow = 0;
2995
2996 if (sqpow > 0)
2997 {
2998 if (sqpow > 1023)
2999 sqpow = 1023;
3000
3001 float scale_factor = 1.0;
3002 for (octave_idx_type i = 0; i < sqpow; i++)
3003 scale_factor *= 2.0;
3004
3005 m = m / scale_factor;
3006 }
3007
3008 // npp, dpp: pade' approx polynomial matrices.
3009
3010 FloatComplexMatrix npp (nc, nc, 0.0);
3011 FloatComplex *pnpp = npp.fortran_vec ();
3012 FloatComplexMatrix dpp = npp;
3013 FloatComplex *pdpp = dpp.fortran_vec ();
3014
3015 // Now powers a^8 ... a^1.
3016
3017 int minus_one_j = -1;
3018 for (octave_idx_type j = 7; j >= 0; j--)
3019 {
3020 for (octave_idx_type i = 0; i < nc; i++)
3021 {
3022 octave_idx_type k = i * nc + i;
3023 pnpp[k] += padec[j];
3024 pdpp[k] += minus_one_j * padec[j];
3025 }
3026
3027 npp = m * npp;
3028 pnpp = npp.fortran_vec ();
3029
3030 dpp = m * dpp;
3031 pdpp = dpp.fortran_vec ();
3032
3033 minus_one_j *= -1;
3034 }
3035
3036 // Zero power.
3037
3038 dpp = -dpp;
3039 for (octave_idx_type j = 0; j < nc; j++)
3040 {
3041 npp.elem (j, j) += 1.0;
3042 dpp.elem (j, j) += 1.0;
3043 }
3044
3045 // Compute pade approximation = inverse (dpp) * npp.
3046
3047 float rcon;
3048 retval = dpp.solve (npp, info, rcon, solve_singularity_warning);
3049
3050 if (info < 0)
3051 return retval;
3052
3053 // Reverse preconditioning step 3: repeated squaring.
3054
3055 while (sqpow)
3056 {
3057 retval = retval * retval;
3058 sqpow--;
3059 }
3060
3061 // Reverse preconditioning step 2: inverse balancing.
3062 // Done in two steps: inverse scaling, then inverse permutation
3063
3064 // inverse scaling (diagonal transformation)
3065 for (octave_idx_type i = 0; i < nc; i++)
3066 for (octave_idx_type j = 0; j < nc; j++)
3067 retval(i,j) *= dscale(i) / dscale(j);
3068
3069 OCTAVE_QUIT;
3070
3071 // construct balancing permutation vector
3072 Array<octave_idx_type> iperm (nc);
3073 for (octave_idx_type i = 0; i < nc; i++)
3074 iperm(i) = i; // identity permutation
3075
3076 // trailing permutations must be done in reverse order
3077 for (octave_idx_type i = nc - 1; i >= ihi; i--)
3078 {
3079 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
3080 octave_idx_type tmp = iperm(i);
3081 iperm(i) = iperm(swapidx);
3082 iperm(swapidx) = tmp;
3083 }
3084
3085 // leading permutations in forward order
3086 for (octave_idx_type i = 0; i < (ilo-1); i++)
3087 {
3088 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
3089 octave_idx_type tmp = iperm(i);
3090 iperm(i) = iperm (swapidx);
3091 iperm(swapidx) = tmp;
3092 }
3093
3094 // construct inverse balancing permutation vector
3095 Array<octave_idx_type> invpvec (nc);
3096 for (octave_idx_type i = 0; i < nc; i++)
3097 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method
3098
3099 OCTAVE_QUIT;
3100
3101 FloatComplexMatrix tmpMat = retval;
3102 for (octave_idx_type i = 0; i < nc; i++)
3103 for (octave_idx_type j = 0; j < nc; j++)
3104 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
3105
3106 // Reverse preconditioning step 1: fix trace normalization.
3107
3108 return exp (trshift) * retval;
3109 } 2921 }
3110 2922
3111 // column vector by row vector -> matrix operations 2923 // column vector by row vector -> matrix operations
3112 2924
3113 FloatComplexMatrix 2925 FloatComplexMatrix