# HG changeset patch # User Jaroslav Hajek # Date 1223793895 -7200 # Node ID 851803f7bb4d9521f80c57478d74d5eb4e4a3e17 # Parent a10397d26114998bca6c7c5570eb2feb40f77b91 improve inverse preconditioning according to Marco Caliari diff -r a10397d26114 -r 851803f7bb4d liboctave/CMatrix.cc --- 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 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 (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 (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 (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; } diff -r a10397d26114 -r 851803f7bb4d liboctave/ChangeLog --- 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 + + * 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 * sparse-util.h (SparseCholPrint): Change char * argument to const diff -r a10397d26114 -r 851803f7bb4d liboctave/dMatrix.cc --- 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 (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 (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; diff -r a10397d26114 -r 851803f7bb4d liboctave/fCMatrix.cc --- 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 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 (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 (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 (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; } diff -r a10397d26114 -r 851803f7bb4d liboctave/fMatrix.cc --- 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 (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 (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;