Mercurial > hg > octave-lyh

comparison liboctave/CMatrix.cc @ 7400:f9df7f7520e7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2008-01-18 08:25:24 by jwe]
author jwe
date Fri, 18 Jan 2008 08:25:24 +0000
parents 7da4a5262e2e
children a7a987b229b7
comparison
equal deleted inserted replaced
7399:c1a3d6c7d2fb 7400:f9df7f7520e7
2748 4.8562548562548563e-6, 2748 4.8562548562548563e-6,
2749 1.3875013875013875e-7, 2749 1.3875013875013875e-7,
2750 1.9270852604185938e-9, 2750 1.9270852604185938e-9,
2751 }; 2751 };
2752 2752
2753 static void
2754 solve_singularity_warning (double rcond)
2755 {
2756 (*current_liboctave_warning_handler)
2757 ("singular matrix encountered in expm calculation, rcond = %g",
2758 rcond);
2759 }
2760
2753 ComplexMatrix 2761 ComplexMatrix
2754 ComplexMatrix::expm (void) const 2762 ComplexMatrix::expm (void) const
2755 { 2763 {
2756 ComplexMatrix retval; 2764 ComplexMatrix retval;
2757 2765
2841 if (sqpow < 0) 2849 if (sqpow < 0)
2842 sqpow = 0; 2850 sqpow = 0;
2843 2851
2844 if (sqpow > 0) 2852 if (sqpow > 0)
2845 { 2853 {
2854 if (sqpow > 1023)
2855 sqpow = 1023;
2856
2846 double scale_factor = 1.0; 2857 double scale_factor = 1.0;
2847 for (octave_idx_type i = 0; i < sqpow; i++) 2858 for (octave_idx_type i = 0; i < sqpow; i++)
2848 scale_factor *= 2.0; 2859 scale_factor *= 2.0;
2849 2860
2850 m = m / scale_factor; 2861 m = m / scale_factor;
2887 dpp.elem (j, j) += 1.0; 2898 dpp.elem (j, j) += 1.0;
2888 } 2899 }
2889 2900
2890 // Compute pade approximation = inverse (dpp) * npp. 2901 // Compute pade approximation = inverse (dpp) * npp.
2891 2902
2892 retval = dpp.solve (npp); 2903 double rcond;
2893 2904 retval = dpp.solve (npp, info, rcond, solve_singularity_warning);
2905
2906 if (info < 0)
2907 return retval;
2908
2894 // Reverse preconditioning step 3: repeated squaring. 2909 // Reverse preconditioning step 3: repeated squaring.
2895 2910
2896 while (sqpow) 2911 while (sqpow)
2897 { 2912 {
2898 retval = retval * retval; 2913 retval = retval * retval;