Mercurial > hg > octave-nkf
changeset 3468:a2dc6de198f9
[project @ 2000-01-21 22:13:13 by jwe]
| author | jwe |
|---|---|
| date | Fri, 21 Jan 2000 22:13:13 +0000 |
| parents | 63378e4a34f2 |
| children | fe0c38ca9d82 |
| files | liboctave/CMatrix.cc liboctave/ChangeLog liboctave/dMatrix.cc |
| diffstat | 3 files changed, 75 insertions(+), 63 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CMatrix.cc +++ b/liboctave/CMatrix.cc @@ -1606,11 +1606,13 @@ Complex *mp = m.fortran_vec (); int info, ilo, ihi,ilos,ihis; - Array<double> dpermute(nc); - Array<double> dscale(nc); - - // FIXME: should pass job as a parameter in expm - char job = 'P'; // Permute first + Array<double> dpermute (nc); + Array<double> dscale (nc); + + // XXX FIXME XXX -- should pass job as a parameter in expm + + // Permute first + char job = 'P'; F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi, dpermute.fortran_vec(), info, 1L, 1L)); @@ -1620,7 +1622,8 @@ return retval; } - job = 'S'; // then scale + // then scale + job = 'S'; F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilos, ihis, dscale.fortran_vec(), info, 1L, 1L)); @@ -1701,43 +1704,42 @@ // Done in two steps: inverse scaling, then inverse permutation // inverse scaling (diagonal transformation) - int ii, jj; - for(ii = 0; ii < nc ; ii++) - for(jj=0; jj < nc ; jj++) - retval(ii,jj) *= dscale(ii)/dscale(jj); + for (int i = 0; i < nc; i++) + for (int j = 0; j < nc; j++) + retval(i,j) *= dscale(i) / dscale(j); // construct balancing permutation vector - Array<int> ipermute(nc); - for(ii=0 ; ii < nc ; ii++) - ipermute(ii) = ii; // initialize to identity permutation + Array<int> ipermute (nc); + for (int i = 0; i < nc; i++) + ipermute(i) = i; // initialize to identity permutation // leading permutations in forward order - for( ii = 0 ; ii < (ilo-1) ; ii++) - { - int swapidx = ( (int) dpermute(ii) ) -1; - int tmp = ipermute(ii); - ipermute(ii) = ipermute( swapidx ); - ipermute(swapidx) = tmp; - } + for (int i = 0; i < (ilo-1); i++) + { + int swapidx = static_cast<int> (dpermute(i)) - 1; + int tmp = ipermute(i); + ipermute(i) = ipermute(swapidx); + ipermute(swapidx) = tmp; + } // trailing permutations must be done in reverse order - for( ii = nc-1 ; ii >= ihi ; ii--) - { - int swapidx = ( (int) dpermute(ii) ) -1; - int tmp = ipermute(ii); - ipermute(ii) = ipermute( swapidx ); - ipermute(swapidx) = tmp; - } + for (int i = nc - 1; i >= ihi; i--) + { + int swapidx = static_cast<int> (dpermute(i)) - 1; + int tmp = ipermute(i); + ipermute(i) = ipermute(swapidx); + ipermute(swapidx) = tmp; + } // construct inverse balancing permutation vector - Array<int> invpvec(nc); - for( ii = 0 ; ii < nc ; ii++) - invpvec(ipermute(ii)) = ii ; // Thanks to R. A. Lippert for this method + Array<int> invpvec (nc); + for (int i = 0; i < nc; i++) + invpvec(ipermute(i)) = i; // Thanks to R. A. Lippert for this method ComplexMatrix tmpMat = retval; - for( ii = 0 ; ii < nc ; ii ++) - for( jj= 0 ; jj < nc ; jj++ ) - retval(ii,jj) = tmpMat(invpvec(ii),invpvec(jj)); + for (int i = 0; i < nc; i++) + for (int j = 0; j < nc; j++) + retval(i,j) = tmpMat(invpvec(i),invpvec(j)); // Reverse preconditioning step 1: fix trace normalization.
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,10 @@ +2000-01-21 John W. Eaton <jwe@bevo.che.wisc.edu> + + * CMatrix.cc (ComplexMatrix::expm): Apply permutation and scaling + operations directly in step 2 and reverse step 2. + * dMatrix.cc (Matrix::expm): Apply permutation and scaling + operations directly in step 2 and reverse step 2. + 2000-01-20 John W. Eaton <jwe@bevo.che.wisc.edu> * oct-time.h, oct-time.cc (octave_strptime): New class.
--- a/liboctave/dMatrix.cc +++ b/liboctave/dMatrix.cc @@ -1356,16 +1356,18 @@ double *p_m = m.fortran_vec (); - int info, ilo, ihi,ilos,ihis, ii, jj; - Array<double> dpermute(nc); - Array<double> dscale(nc); + int info, ilo, ihi, ilos, ihis; + Array<double> dpermute (nc); + Array<double> dscale (nc); double *dp; - char job = 'P'; // permutation first + // permutation first + char job = 'P'; dp = dpermute.fortran_vec(); F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilo, ihi, dp, info, 1L, 1L)); - job = 'S'; // then scaling + // then scaling + job = 'S'; dp = dscale.fortran_vec(); F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilos, ihis, dp, info, 1L, 1L)); @@ -1445,41 +1447,42 @@ // Reverse preconditioning step 2: inverse balancing. // apply inverse scaling to computed exponential - for(ii = 0; ii < nc ; ii++) - for(jj=0; jj < nc ; jj++) - retval(ii,jj) *= dscale(ii)/dscale(jj); + for (int i = 0; i < nc; i++) + for (int j = 0; j < nc; j++) + retval(i,j) *= dscale(i) / dscale(j); // construct balancing permutation vector - Array<int> ipermute(nc); - for(ii=0 ; ii < nc ; ii++) ipermute(ii) = ii; // identity permutation + Array<int> ipermute (nc); + for (int i = 0; i < nc; i++) + ipermute(i) = i; // identity permutation // leading permutations in forward order - for( ii = 0 ; ii < (ilo-1) ; ii++) - { - int swapidx = ( (int) dpermute(ii) ) -1; - int tmp = ipermute(ii); - ipermute(ii) = ipermute( swapidx ); - ipermute(swapidx) = tmp; - } + for (int i = 0; i < (ilo-1); i++) + { + int swapidx = static_cast<int> (dpermute(i)) - 1; + int tmp = ipermute(i); + ipermute(i) = ipermute (swapidx); + ipermute(swapidx) = tmp; + } // trailing permutations must be done in reverse order - for( ii = nc-1 ; ii >= ihi ; ii--) - { - int swapidx = ( (int) dpermute(ii) ) -1; - int tmp = ipermute(ii); - ipermute(ii) = ipermute( swapidx ); - ipermute(swapidx) = tmp; - } + for (int i = nc - 1; i >= ihi; i--) + { + int swapidx = static_cast<int> (dpermute(i)) - 1; + int tmp = ipermute(i); + ipermute(i) = ipermute(swapidx); + ipermute(swapidx) = tmp; + } // construct inverse balancing permutation vector - Array<int> invpvec(nc); - for( ii = 0 ; ii < nc ; ii++) - invpvec(ipermute(ii)) = ii ; // Thanks to R. A. Lippert for this method + Array<int> invpvec (nc); + for (int i = 0; i < nc; i++) + invpvec(ipermute(i)) = i; // Thanks to R. A. Lippert for this method Matrix tmpMat = retval; - for( ii = 0 ; ii < nc ; ii ++) - for( jj= 0 ; jj < nc ; jj++ ) - retval(ii,jj) = tmpMat(invpvec(ii),invpvec(jj)); + for (int i = 0; i < nc; i++) + for (int j = 0; j < nc; j++) + retval(i,j) = tmpMat(invpvec(i),invpvec(j)); // Reverse preconditioning step 1: fix trace normalization.