Mercurial > hg > octave-terminal
changeset 11875:f53f57d2ee51 release-3-0-x
improve inverse preconditioning according to Marco Caliari
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Sun, 12 Oct 2008 10:31:14 +0200 |
parents | 757ca614e918 |
children | b1122c225f6c |
files | liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc |
diffstat | 3 files changed, 33 insertions(+), 58 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -2964,21 +2964,30 @@ // inverse scaling (diagonal transformation) for (octave_idx_type i = 0; i < nc; i++) for (octave_idx_type j = 0; j < nc; j++) - retval(i,j) *= dscale(i) / dscale(j); + retval(i,j) *= dscale(i) / dscale(j); OCTAVE_QUIT; // construct balancing permutation vector Array<octave_idx_type> iperm (nc); for (octave_idx_type i = 0; i < nc; i++) - iperm(i) = i; // initialize to identity permutation + iperm(i) = i; // identity permutation + + // trailing permutations must be done in reverse order + for (octave_idx_type i = nc - 1; i >= ihi; i--) + { + octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1; + octave_idx_type tmp = iperm(i); + iperm(i) = iperm(swapidx); + iperm(swapidx) = tmp; + } // leading permutations in forward order for (octave_idx_type i = 0; i < (ilo-1); i++) { octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1; octave_idx_type tmp = iperm(i); - iperm(i) = iperm(swapidx); + iperm(i) = iperm (swapidx); iperm(swapidx) = tmp; } @@ -2994,32 +3003,7 @@ for (octave_idx_type j = 0; j < nc; j++) retval(i,j) = tmpMat(invpvec(i),invpvec(j)); - OCTAVE_QUIT; - - for (octave_idx_type i = 0; i < nc; i++) - iperm(i) = i; // initialize to identity permutation - - // trailing permutations must be done in reverse order - for (octave_idx_type i = nc - 1; i >= ihi; i--) - { - octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1; - octave_idx_type tmp = iperm(i); - iperm(i) = iperm(swapidx); - iperm(swapidx) = tmp; - } - - // construct inverse balancing permutation vector - for (octave_idx_type i = 0; i < nc; i++) - invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method - - OCTAVE_QUIT; - - tmpMat = retval; - for (octave_idx_type i = 0; i < nc; i++) - for (octave_idx_type j = 0; j < nc; j++) - retval(i,j) = tmpMat(invpvec(i),invpvec(j)); - - // Reverse preconditioning step 1: fix trace normalization. + // Reverse preconditioning step 1: fix trace normalization. return exp (trshift) * retval; }
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,11 @@ +2008-10-12 Jaroslav Hajek <highegg@gmail.com> + + * CSparse.cc (ComplexMatrix::expm): Improve inverse preconditioning + according to Marco Caliari. + * dSparse.cc (Matrix::expm): Likewise. + * fCSparse.cc (FloatComplexMatrix::expm): Likewise. + * fSparse.cc (FloatMatrix::expm): Likewise. + 2008-10-10 Jaroslav Hajek <highegg@gmail.com> * sparse-util.h (SparseCholPrint): Change char * argument to const
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -2589,7 +2589,7 @@ // apply inverse scaling to computed exponential for (octave_idx_type i = 0; i < nc; i++) for (octave_idx_type j = 0; j < nc; j++) - retval(i,j) *= dscale(i) / dscale(j); + retval(i,j) *= dscale(i) / dscale(j); OCTAVE_QUIT; @@ -2598,6 +2598,15 @@ for (octave_idx_type i = 0; i < nc; i++) iperm(i) = i; // identity permutation + // trailing permutations must be done in reverse order + for (octave_idx_type i = nc - 1; i >= ihi; i--) + { + octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1; + octave_idx_type tmp = iperm(i); + iperm(i) = iperm(swapidx); + iperm(swapidx) = tmp; + } + // leading permutations in forward order for (octave_idx_type i = 0; i < (ilo-1); i++) { @@ -2613,39 +2622,13 @@ invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method OCTAVE_QUIT; - + Matrix tmpMat = retval; for (octave_idx_type i = 0; i < nc; i++) for (octave_idx_type j = 0; j < nc; j++) retval(i,j) = tmpMat(invpvec(i),invpvec(j)); - OCTAVE_QUIT; - - for (octave_idx_type i = 0; i < nc; i++) - iperm(i) = i; // identity permutation - - // trailing permutations must be done in reverse order - for (octave_idx_type i = nc - 1; i >= ihi; i--) - { - octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1; - octave_idx_type tmp = iperm(i); - iperm(i) = iperm(swapidx); - iperm(swapidx) = tmp; - } - - // construct inverse balancing permutation vector - for (octave_idx_type i = 0; i < nc; i++) - invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method - - OCTAVE_QUIT; - - tmpMat = retval; - for (octave_idx_type i = 0; i < nc; i++) - for (octave_idx_type j = 0; j < nc; j++) - retval(i,j) = tmpMat(invpvec(i),invpvec(j)); - - // Reverse preconditioning step 1: fix trace normalization. - + // Reverse preconditioning step 1: fix trace normalization. if (trshift > 0.0) retval = exp (trshift) * retval;