Mercurial > hg > octave-nkf
diff liboctave/CMatrix.cc @ 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 | 0b9c56b6bf0e |
children |
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; }