Mercurial > octave-nkf
changeset 8211:851803f7bb4d
improve inverse preconditioning according to Marco Caliari
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Sun, 12 Oct 2008 08:44:55 +0200 |
| parents | a10397d26114 |
| children | ebf6f6a0f9a7 |
| files | liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc liboctave/fCMatrix.cc liboctave/fMatrix.cc |
| diffstat | 5 files changed, 58 insertions(+), 115 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc Fri Oct 10 15:10:58 2008 -0400 +++ b/liboctave/CMatrix.cc Sun Oct 12 08:44:55 2008 +0200 @@ -3027,21 +3027,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; } @@ -3057,32 +3066,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 Fri Oct 10 15:10:58 2008 -0400 +++ b/liboctave/ChangeLog Sun Oct 12 08:44:55 2008 +0200 @@ -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 Fri Oct 10 15:10:58 2008 -0400 +++ b/liboctave/dMatrix.cc Sun Oct 12 08:44:55 2008 +0200 @@ -2659,7 +2659,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; @@ -2668,6 +2668,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++) { @@ -2683,39 +2692,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;
--- a/liboctave/fCMatrix.cc Fri Oct 10 15:10:58 2008 -0400 +++ b/liboctave/fCMatrix.cc Sun Oct 12 08:44:55 2008 +0200 @@ -3020,21 +3020,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; } @@ -3050,32 +3059,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/fMatrix.cc Fri Oct 10 15:10:58 2008 -0400 +++ b/liboctave/fMatrix.cc Sun Oct 12 08:44:55 2008 +0200 @@ -2658,7 +2658,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; @@ -2667,6 +2667,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++) { @@ -2682,38 +2691,13 @@ invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method OCTAVE_QUIT; - + FloatMatrix 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 = expf (trshift) * retval;