Mercurial > hg > octave-nkf

comparison liboctave/CMatrix.cc @ 11875:f53f57d2ee51 release-3-0-x

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
11874:757ca614e918 11875:f53f57d2ee51
2962 // Done in two steps: inverse scaling, then inverse permutation 2962 // Done in two steps: inverse scaling, then inverse permutation
2963 2963
2964 // inverse scaling (diagonal transformation) 2964 // inverse scaling (diagonal transformation)
2965 for (octave_idx_type i = 0; i < nc; i++) 2965 for (octave_idx_type i = 0; i < nc; i++)
2966 for (octave_idx_type j = 0; j < nc; j++) 2966 for (octave_idx_type j = 0; j < nc; j++)
2967 retval(i,j) *= dscale(i) / dscale(j); 2967 retval(i,j) *= dscale(i) / dscale(j);
2968 2968
2969 OCTAVE_QUIT; 2969 OCTAVE_QUIT;
2970 2970
2971 // construct balancing permutation vector 2971 // construct balancing permutation vector
2972 Array<octave_idx_type> iperm (nc); 2972 Array<octave_idx_type> iperm (nc);
2973 for (octave_idx_type i = 0; i < nc; i++) 2973 for (octave_idx_type i = 0; i < nc; i++)
2974 iperm(i) = i; // initialize to identity permutation 2974 iperm(i) = i; // identity permutation
2975 2975
2976 // leading permutations in forward order 2976 // trailing permutations must be done in reverse order
2977 for (octave_idx_type i = 0; i < (ilo-1); i++) 2977 for (octave_idx_type i = nc - 1; i >= ihi; i--)
2978 { 2978 {
2979 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1; 2979 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
2980 octave_idx_type tmp = iperm(i); 2980 octave_idx_type tmp = iperm(i);
2981 iperm(i) = iperm(swapidx); 2981 iperm(i) = iperm(swapidx);
2982 iperm(swapidx) = tmp; 2982 iperm(swapidx) = tmp;
2983 } 2983 }
2984 2984
2985 // leading permutations in forward order
2986 for (octave_idx_type i = 0; i < (ilo-1); i++)
2987 {
2988 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
2989 octave_idx_type tmp = iperm(i);
2990 iperm(i) = iperm (swapidx);
2991 iperm(swapidx) = tmp;
2992 }
2993
2985 // construct inverse balancing permutation vector 2994 // construct inverse balancing permutation vector
2986 Array<octave_idx_type> invpvec (nc); 2995 Array<octave_idx_type> invpvec (nc);
2987 for (octave_idx_type i = 0; i < nc; i++) 2996 for (octave_idx_type i = 0; i < nc; i++)
2988 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method 2997 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method
2989 2998
2992 ComplexMatrix tmpMat = retval; 3001 ComplexMatrix tmpMat = retval;
2993 for (octave_idx_type i = 0; i < nc; i++) 3002 for (octave_idx_type i = 0; i < nc; i++)
2994 for (octave_idx_type j = 0; j < nc; j++) 3003 for (octave_idx_type j = 0; j < nc; j++)
2995 retval(i,j) = tmpMat(invpvec(i),invpvec(j)); 3004 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
2996 3005
2997 OCTAVE_QUIT; 3006 // Reverse preconditioning step 1: fix trace normalization.
2998
2999 for (octave_idx_type i = 0; i < nc; i++)
3000 iperm(i) = i; // initialize to identity permutation
3001
3002 // trailing permutations must be done in reverse order
3003 for (octave_idx_type i = nc - 1; i >= ihi; i--)
3004 {
3005 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
3006 octave_idx_type tmp = iperm(i);
3007 iperm(i) = iperm(swapidx);
3008 iperm(swapidx) = tmp;
3009 }
3010
3011 // construct inverse balancing permutation vector
3012 for (octave_idx_type i = 0; i < nc; i++)
3013 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method
3014
3015 OCTAVE_QUIT;
3016
3017 tmpMat = retval;
3018 for (octave_idx_type i = 0; i < nc; i++)
3019 for (octave_idx_type j = 0; j < nc; j++)
3020 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
3021
3022 // Reverse preconditioning step 1: fix trace normalization.
3023 3007
3024 return exp (trshift) * retval; 3008 return exp (trshift) * retval;
3025 } 3009 }
3026 3010
3027 // column vector by row vector -> matrix operations 3011 // column vector by row vector -> matrix operations