Mercurial > hg > octave-lyh

comparison liboctave/dMatrix.cc @ 8211:851803f7bb4d

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 08:44:55 +0200
parents 935be827eaf8
children 64cf956a109c
comparison
equal deleted inserted replaced
8210:a10397d26114 8211:851803f7bb4d
2657 2657
2658 // Reverse preconditioning step 2: inverse balancing. 2658 // Reverse preconditioning step 2: inverse balancing.
2659 // apply inverse scaling to computed exponential 2659 // apply inverse scaling to computed exponential
2660 for (octave_idx_type i = 0; i < nc; i++) 2660 for (octave_idx_type i = 0; i < nc; i++)
2661 for (octave_idx_type j = 0; j < nc; j++) 2661 for (octave_idx_type j = 0; j < nc; j++)
2662 retval(i,j) *= dscale(i) / dscale(j); 2662 retval(i,j) *= dscale(i) / dscale(j);
2663 2663
2664 OCTAVE_QUIT; 2664 OCTAVE_QUIT;
2665 2665
2666 // construct balancing permutation vector 2666 // construct balancing permutation vector
2667 Array<octave_idx_type> iperm (nc); 2667 Array<octave_idx_type> iperm (nc);
2668 for (octave_idx_type i = 0; i < nc; i++) 2668 for (octave_idx_type i = 0; i < nc; i++)
2669 iperm(i) = i; // identity permutation 2669 iperm(i) = i; // identity permutation
2670 2670
2671 // trailing permutations must be done in reverse order
2672 for (octave_idx_type i = nc - 1; i >= ihi; i--)
2673 {
2674 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
2675 octave_idx_type tmp = iperm(i);
2676 iperm(i) = iperm(swapidx);
2677 iperm(swapidx) = tmp;
2678 }
2679
2671 // leading permutations in forward order 2680 // leading permutations in forward order
2672 for (octave_idx_type i = 0; i < (ilo-1); i++) 2681 for (octave_idx_type i = 0; i < (ilo-1); i++)
2673 { 2682 {
2674 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1; 2683 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
2675 octave_idx_type tmp = iperm(i); 2684 octave_idx_type tmp = iperm(i);
2681 Array<octave_idx_type> invpvec (nc); 2690 Array<octave_idx_type> invpvec (nc);
2682 for (octave_idx_type i = 0; i < nc; i++) 2691 for (octave_idx_type i = 0; i < nc; i++)
2683 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method 2692 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method
2684 2693
2685 OCTAVE_QUIT; 2694 OCTAVE_QUIT;
2686 2695
2687 Matrix tmpMat = retval; 2696 Matrix tmpMat = retval;
2688 for (octave_idx_type i = 0; i < nc; i++) 2697 for (octave_idx_type i = 0; i < nc; i++)
2689 for (octave_idx_type j = 0; j < nc; j++) 2698 for (octave_idx_type j = 0; j < nc; j++)
2690 retval(i,j) = tmpMat(invpvec(i),invpvec(j)); 2699 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
2691 2700
2692 OCTAVE_QUIT; 2701 // Reverse preconditioning step 1: fix trace normalization.
2693
2694 for (octave_idx_type i = 0; i < nc; i++)
2695 iperm(i) = i; // identity permutation
2696
2697 // trailing permutations must be done in reverse order
2698 for (octave_idx_type i = nc - 1; i >= ihi; i--)
2699 {
2700 octave_idx_type swapidx = static_cast<octave_idx_type> (dpermute(i)) - 1;
2701 octave_idx_type tmp = iperm(i);
2702 iperm(i) = iperm(swapidx);
2703 iperm(swapidx) = tmp;
2704 }
2705
2706 // construct inverse balancing permutation vector
2707 for (octave_idx_type i = 0; i < nc; i++)
2708 invpvec(iperm(i)) = i; // Thanks to R. A. Lippert for this method
2709
2710 OCTAVE_QUIT;
2711
2712 tmpMat = retval;
2713 for (octave_idx_type i = 0; i < nc; i++)
2714 for (octave_idx_type j = 0; j < nc; j++)
2715 retval(i,j) = tmpMat(invpvec(i),invpvec(j));
2716
2717 // Reverse preconditioning step 1: fix trace normalization.
2718
2719 if (trshift > 0.0) 2702 if (trshift > 0.0)
2720 retval = exp (trshift) * retval; 2703 retval = exp (trshift) * retval;
2721 2704
2722 return retval; 2705 return retval;
2723 } 2706 }