AnasaziModalSolverUtils.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef ANASAZI_MODAL_SOLVER_UTILS_HPP
00030 #define ANASAZI_MODAL_SOLVER_UTILS_HPP
00031 
00048 #include "AnasaziConfigDefs.hpp"
00049 #include "AnasaziMultiVecTraits.hpp"
00050 #include "AnasaziOperatorTraits.hpp"
00051 #include "Teuchos_ScalarTraits.hpp"
00052 
00053 #include "AnasaziOutputManager.hpp"
00054 #include "Teuchos_BLAS.hpp"
00055 #include "Teuchos_LAPACK.hpp"
00056 #include "Teuchos_SerialDenseMatrix.hpp"
00057 
00058 namespace Anasazi {
00059 
00060   template<class ScalarType, class MV, class OP>
00061   class ModalSolverUtils 
00062   {  
00063   public:
00064     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00065     typedef typename Teuchos::ScalarTraits<ScalarType>  SCT;
00066     
00068 
00069 
00071 
00073     ModalSolverUtils( const Teuchos::RefCountPtr<OutputManager<ScalarType> > &om ); 
00074     
00076     virtual ~ModalSolverUtils() {};
00077     
00079     
00081 
00082     
00084     int sortScalars(int n, ScalarType *y, int *perm = 0) const;
00085     
00087     int sortScalars_Vectors(int n, ScalarType* lambda, MV* Q, std::vector<MagnitudeType>* resids = 0) const;
00088 
00090     void permuteVectors(const int n, const std::vector<int> &perm, MV &Q, std::vector<MagnitudeType>* resids = 0) const;
00091 
00093     void permuteVectors(const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q) const;
00094 
00096 
00098 
00099 
00101 
00122     int massOrthonormalize(MV &X, MV &MX, const OP *M, const MV &Q, int howMany,
00123                            int orthoType = 0, ScalarType kappa = 1.5625) const;
00124     
00125 
00127 
00153     int directSolver(int size, const Teuchos::SerialDenseMatrix<int,ScalarType> &KK, 
00154                      const Teuchos::SerialDenseMatrix<int,ScalarType> *MM,
00155                      Teuchos::SerialDenseMatrix<int,ScalarType> *EV,
00156                      std::vector<MagnitudeType>* theta,
00157                      int* nev, int esType = 0) const;
00159 
00161 
00162 
00164 
00166     MagnitudeType errorOrthogonality(const MV *X, const MV *R, const OP *M = 0) const;
00167     
00169 
00171     MagnitudeType errorOrthonormality(const MV *X, const OP *M = 0) const;
00172     
00174 
00176     MagnitudeType errorEquality(const MV *X, const MV *MX, const OP *M = 0) const;
00177     
00179     
00180   private:
00181 
00182     // Reference counted pointer to output manager used by eigensolver.
00183     Teuchos::RefCountPtr<OutputManager<ScalarType> > om_;
00184 
00186 
00187 
00188     typedef MultiVecTraits<ScalarType,MV> MVT;
00189     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00190 
00192   };
00193 
00194   //-----------------------------------------------------------------------------
00195   // 
00196   //  CONSTRUCTOR
00197   //
00198   //-----------------------------------------------------------------------------  
00199 
00200   template<class ScalarType, class MV, class OP>
00201   ModalSolverUtils<ScalarType, MV, OP>::ModalSolverUtils( const Teuchos::RefCountPtr<OutputManager<ScalarType> > &om ) 
00202     : om_(om)
00203   {}
00204 
00205   //-----------------------------------------------------------------------------
00206   // 
00207   //  SORTING METHODS
00208   //
00209   //-----------------------------------------------------------------------------
00210   
00211   template<class ScalarType, class MV, class OP>
00212   int ModalSolverUtils<ScalarType, MV, OP>::sortScalars( int n, ScalarType *y, int *perm) const 
00213   {
00214     // Sort a vector into increasing order of algebraic values
00215     //
00216     // Input:
00217     //
00218     // n    (integer ) = Size of the array (input)
00219     // y    (ScalarType* ) = Array of length n to be sorted (input/output)
00220     // perm (integer*) = Array of length n with the permutation (input/output)
00221     //                   Optional argument
00222     
00223     int i, j;
00224     int igap = n / 2;
00225     
00226     if (igap == 0) {
00227       if ((n > 0) && (perm != 0)) {
00228         perm[0] = 0;
00229       }
00230       return 0;
00231     }
00232     
00233     if (perm) {
00234       for (i = 0; i < n; ++i)
00235         perm[i] = i;
00236     }
00237     
00238     while (igap > 0) {
00239       for (i=igap; i<n; ++i) {
00240         for (j=i-igap; j>=0; j-=igap) {
00241           if (y[j] > y[j+igap]) {
00242             ScalarType tmpD = y[j];
00243             y[j] = y[j+igap];
00244             y[j+igap] = tmpD;
00245             if (perm) {
00246               int tmpI = perm[j];
00247               perm[j] = perm[j+igap];
00248               perm[j+igap] = tmpI;
00249             }
00250           }
00251           else {
00252             break;
00253           }
00254         }
00255       }
00256       igap = igap / 2;
00257     }
00258     
00259     return 0;
00260   }
00261   
00262   template<class ScalarType, class MV, class OP>
00263   int ModalSolverUtils<ScalarType, MV, OP>::sortScalars_Vectors( int n, ScalarType *lambda, MV *Q, 
00264                                                                  std::vector<MagnitudeType> *resids ) const
00265   {
00266     // This routines sorts the scalars (stored in lambda) in ascending order.
00267     // The associated vectors (stored in Q) are accordingly ordered.
00268     
00269     int info = 0;
00270     int i, j;
00271     std::vector<int> index(1);
00272     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00273     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00274     
00275     if ( n > MVT::GetNumberVecs( *Q ) ) { return -1; }
00276 
00277     int igap = n / 2;
00278     
00279     if ( Q ) {
00280       while (igap > 0) {
00281         for (i=igap; i < n; ++i) {
00282           for (j=i-igap; j>=0; j-=igap) {
00283             if (lambda[j] > lambda[j+igap]) {
00284 
00285               // Swap two scalars
00286               std::swap<ScalarType>(lambda[j],lambda[j+igap]);
00287 
00288               // Swap residuals (if they exist)
00289               if (resids) {
00290                 std::swap<MagnitudeType>((*resids)[j],(*resids)[j+igap]);
00291               }
00292 
00293               // Swap corresponding vectors
00294               index[0] = j;
00295               Teuchos::RefCountPtr<MV> tmpQ = MVT::CloneCopy( *Q, index );
00296               Teuchos::RefCountPtr<MV> tmpQj = MVT::CloneView( *Q, index );
00297               index[0] = j + igap;
00298               Teuchos::RefCountPtr<MV> tmpQgap = MVT::CloneView( *Q, index );
00299               MVT::MvAddMv( one, *tmpQgap, zero, *tmpQgap, *tmpQj );
00300               MVT::MvAddMv( one, *tmpQ, zero, *tmpQ, *tmpQgap );
00301             } 
00302             else {
00303               break;
00304             }
00305           }
00306         }
00307         igap = igap / 2;
00308       } // while (igap > 0)
00309     } // if ( Q )
00310     else {
00311       while (igap > 0) {
00312         for (i=igap; i < n; ++i) {
00313           for (j=i-igap; j>=0; j-=igap) {
00314             if (lambda[j] > lambda[j+igap]) {
00315               // Swap two scalars
00316               std::swap<ScalarType>(lambda[j],lambda[j+igap]);
00317               
00318               // Swap residuals (if they exist)
00319               if (resids) {
00320                 std::swap<MagnitudeType>((*resids)[j],(*resids)[j+igap]);
00321               }              
00322             } 
00323             else {
00324               break;
00325             }
00326           }
00327         }
00328         igap = igap / 2;
00329       } // while (igap > 0)
00330     } // if ( Q )
00331     
00332     return info;  
00333   }
00334 
00336   // permuteVectors for MV
00337   template<class ScalarType, class MV, class OP>
00338   void ModalSolverUtils<ScalarType, MV, OP>::permuteVectors(
00339               const int n,
00340               const std::vector<int> &perm, 
00341               MV &Q, 
00342               std::vector<MagnitudeType>* resids) const
00343   {
00344     // Permute the vectors according to the permutation vector \c perm, and
00345     // optionally the residual vector \c resids
00346     
00347     int i, j;
00348     std::vector<int> permcopy(perm), swapvec(n-1);
00349     std::vector<int> index(1);
00350     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00351     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00352 
00353     TEST_FOR_EXCEPTION(n > MVT::GetNumberVecs(Q), std::invalid_argument, "Anasazi::ModalSolverUtils::permuteVectors(): argument n larger than width of input multivector.");
00354 
00355     // We want to recover the elementary permutations (individual swaps) 
00356     // from the permutation vector. Do this by constructing the inverse
00357     // of the permutation, by sorting them to {1,2,...,n}, and recording
00358     // the elementary permutations of the inverse.
00359     for (i=0; i<n-1; i++) {
00360       //
00361       // find i in the permcopy vector
00362       for (j=i; j<n; j++) {
00363         if (permcopy[j] == i) {
00364           // found it at index j
00365           break;
00366         }
00367         TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument, "Anasazi::ModalSolverUtils::permuteVectors(): permutation index invalid.");
00368       }
00369       //
00370       // Swap two scalars
00371       std::swap<int>( permcopy[j], permcopy[i] );
00372 
00373       swapvec[i] = j;
00374     }
00375       
00376     // now apply the elementary permutations of the inverse in reverse order
00377     for (i=n-2; i>=0; i--) {
00378       j = swapvec[i];
00379       //
00380       // Swap (i,j)
00381       //
00382       // Swap residuals (if they exist)
00383       if (resids) {
00384         std::swap<MagnitudeType>(  (*resids)[i], (*resids)[j] );
00385       }
00386       //
00387       // Swap corresponding vectors
00388       index[0] = j;
00389       Teuchos::RefCountPtr<MV> tmpQ = MVT::CloneCopy( Q, index );
00390       Teuchos::RefCountPtr<MV> tmpQj = MVT::CloneView( Q, index );
00391       index[0] = i;
00392       Teuchos::RefCountPtr<MV> tmpQi = MVT::CloneView( Q, index );
00393       MVT::MvAddMv( one, *tmpQi, zero, *tmpQi, *tmpQj );
00394       MVT::MvAddMv( one, *tmpQ, zero, *tmpQ, *tmpQi );
00395     }
00396   }
00397 
00398 
00400   // permuteVectors for MV
00401   template<class ScalarType, class MV, class OP>
00402   void ModalSolverUtils<ScalarType, MV, OP>::permuteVectors(
00403               const std::vector<int> &perm, 
00404               Teuchos::SerialDenseMatrix<int,ScalarType> &Q) const 
00405   {
00406     // Permute the vectors in Q according to the permutation vector \c perm, and
00407     // optionally the residual vector \c resids
00408     Teuchos::BLAS<int,ScalarType> blas;
00409     const int n = perm.size();
00410     const int m = Q.numRows();
00411     
00412     TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument, "Anasazi::ModalSolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
00413 
00414     // Sort the primitive ritz vectors
00415     Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Q );
00416     for (int i=0; i<n; i++) {
00417       blas.COPY(m, copyQ[perm[i]], 1, Q[i], 1);
00418     }
00419   }
00420   
00421   //-----------------------------------------------------------------------------
00422   // 
00423   //  EIGENSOLVER PROJECTION METHODS
00424   //
00425   //-----------------------------------------------------------------------------
00426   
00427   template<class ScalarType, class MV, class OP>
00428   int ModalSolverUtils<ScalarType, MV, OP>::massOrthonormalize(MV &X, MV &MX, const OP *M, const MV &Q, 
00429                                                                int howMany, int orthoType, ScalarType kappa) const
00430   {
00431     // For the inner product defined by the operator M or the identity (M = 0)
00432     //   -> Orthogonalize X against Q
00433     //   -> Orthonormalize X 
00434     // Modify MX accordingly
00435     //
00436     // Note that when M is 0, MX is not referenced
00437     //
00438     // Parameter variables
00439     //
00440     // X  : Vectors to be transformed
00441     //
00442     // MX : Image of the block vector X by the mass matrix
00443     //
00444     // Q  : Vectors to orthogonalize against
00445     //
00446     // howMany : Number of vectors of X to be orthogonalized
00447     //           If this number is smaller than the total number of vectors in X,
00448     //           then it is assumed that the last "howMany" vectors are not orthogonal
00449     //           while the other vectors in X are orthogonal to Q and orthonormal.
00450     //
00451     // orthoType = 0 (default) > Performs both operations
00452     // orthoType = 1           > Performs Q^T M X = 0
00453     // orthoType = 2           > Performs X^T M X = I
00454     //
00455     // kappa= Coefficient determining when to perform a second Gram-Schmidt step
00456     //        Default value = 1.5625 = (1.25)^2 (as suggested in Parlett's book)
00457     //
00458     // Output the status of the computation
00459     //
00460     // info =   0 >> Success.
00461     //
00462     // info >   0 >> Indicate how many vectors have been tried to avoid rank deficiency for X 
00463     //
00464     // info =  -1 >> Failure >> X has zero columns
00465     //                       >> It happens when # col of X     > # rows of X 
00466     //                       >> It happens when # col of [Q X] > # rows of X 
00467     //                       >> It happens when no good random vectors could be found
00468     
00469     int i;
00470     int info = 0;
00471     ScalarType one     = SCT::one();
00472     MagnitudeType zero = SCT::magnitude(SCT::zero()); 
00473     ScalarType eps     = SCT::eps();
00474     
00475     // Orthogonalize X against Q
00476     //timeProj -= MyWatch.WallTime();
00477     if (orthoType != 2) {
00478       
00479       int xc = howMany;
00480       int xr = MVT::GetNumberVecs( X );
00481       int qc = MVT::GetNumberVecs( Q );
00482       
00483       std::vector<int> index(howMany);
00484       for (i=0; i<howMany; i++) {
00485         index[i] = xr - howMany + i;
00486       }
00487 
00488       Teuchos::RefCountPtr<MV> XX = MVT::CloneView( X, index );
00489       
00490       Teuchos::RefCountPtr<MV> MXX;
00491 
00492       if (M) {
00493         int mxr = MVT::GetNumberVecs( MX );
00494         for (i=0; i<howMany; i++) {
00495           index[i] = mxr - howMany + i;
00496         }
00497         MXX = MVT::CloneView( MX, index );
00498       } 
00499       else {
00500         MXX = MVT::CloneView( X, index );
00501       }
00502 
00503       // Perform the Gram-Schmidt transformation for a block of vectors
00504       
00505       // Compute the initial M-norms
00506       std::vector<ScalarType> oldDot( xc );
00507       MVT::MvDot( *XX, *MXX, &oldDot );
00508       
00509       // Define the product Q^T * (M*X)
00510       // Multiply Q' with MX
00511       Teuchos::SerialDenseMatrix<int,ScalarType> qTmx( qc, xc );
00512       MVT::MvTransMv( one, Q, *MXX, qTmx );
00513       
00514       // Multiply by Q and substract the result in X
00515       MVT::MvTimesMatAddMv( -one, Q, qTmx, one, *XX );
00516       
00517       // Update MX
00518       if (M) {
00519         if (qc >= xc) {
00520           OPT::Apply( *M, *XX, *MXX);
00521         }
00522         else {
00523           Teuchos::RefCountPtr<MV> MQ = MVT::Clone( Q, qc );
00524           OPT::Apply( *M, Q, *MQ );
00525           MVT::MvTimesMatAddMv( -one, *MQ, qTmx, one, *MXX );
00526         }  // if (qc >= xc)
00527       } // if (M)
00528       
00529       // Compute new M-norms
00530       std::vector<ScalarType> newDot(xc);
00531       MVT::MvDot( *XX, *MXX, &newDot );
00532       
00533       int j;
00534       for (j = 0; j < xc; ++j) {
00535         
00536         if ( SCT::magnitude(kappa*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
00537           
00538           // Apply another step of classical Gram-Schmidt
00539           MVT::MvTransMv( one, Q, *MXX, qTmx );
00540           
00541           MVT::MvTimesMatAddMv( -one, Q, qTmx, one, *XX );
00542           
00543           // Update MX
00544           if (M) {
00545             if (qc >= xc) {
00546               OPT::Apply( *M, *XX, *MXX);
00547             }
00548             else {
00549               Teuchos::RefCountPtr<MV> MQ = MVT::Clone( Q, qc );
00550               OPT::Apply( *M, Q, *MQ);
00551               MVT::MvTimesMatAddMv( -one, *MQ, qTmx, one, *MXX );
00552             } // if (qc >= xc)
00553           } // if (M)
00554           
00555           break;
00556         } // if (kappa*newDot[j] < oldDot[j])
00557       } // for (j = 0; j < xc; ++j)
00558       
00559     } // if (orthoType != 2)
00560 
00561     //timeProj += MyWatch.WallTime();
00562     
00563     // Orthonormalize X 
00564     //  timeNorm -= MyWatch.WallTime();
00565     if (orthoType != 1) {
00566       
00567       int j;
00568       int xc = MVT::GetNumberVecs( X );
00569       int xr = MVT::GetVecLength( X );
00570       int shift = (orthoType == 2) ? 0 : MVT::GetNumberVecs( Q );
00571       int mxc = (M) ? MVT::GetNumberVecs( MX ) : xc;
00572 
00573       std::vector<int> index( 1 );      
00574       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00575 
00576       for (j = 0; j < howMany; ++j) {
00577         
00578         int numX = xc - howMany + j;
00579         int numMX = mxc - howMany + j;
00580         
00581         // Put zero vectors in X when we are exceeding the space dimension
00582         if (numX + shift >= xr) {
00583           index[0] = numX;
00584           Teuchos::RefCountPtr<MV> XXj = MVT::CloneView( X, index );
00585           MVT::MvInit( *XXj, zero );
00586           if (M) {
00587             index[0] = numMX;
00588             Teuchos::RefCountPtr<MV> MXXj = MVT::CloneView( MX, index );
00589             MVT::MvInit( *MXXj, zero );
00590           }
00591           info = -1;
00592         }
00593         
00594         // Get a view of the vectors currently being worked on.
00595         index[0] = numX;
00596         Teuchos::RefCountPtr<MV> Xj = MVT::CloneView( X, index );
00597         index[0] = numMX;
00598         Teuchos::RefCountPtr<MV> MXj;
00599         if (M)
00600           MXj = MVT::CloneView( MX, index );
00601         else
00602           MXj = MVT::CloneView( X, index );
00603 
00604         // Get a view of the previous vectors.
00605         std::vector<int> prev_idx( numX );
00606         Teuchos::RefCountPtr<MV> prevXj;
00607 
00608         if (numX > 0) {
00609           for (i=0; i<numX; i++)
00610             prev_idx[i] = i;
00611           prevXj = MVT::CloneView( X, prev_idx );
00612         } 
00613 
00614         // Make storage for these Gram-Schmidt iterations.
00615         Teuchos::SerialDenseMatrix<int,ScalarType> product( numX, 1 );
00616 
00617         int numTrials;
00618         bool rankDef = true;
00619         for (numTrials = 0; numTrials < 10; ++numTrials) {
00620 
00621           //
00622           // Compute M-norm
00623           //
00624           MVT::MvDot( *Xj, *MXj, &oldDot );
00625           //
00626           // Save old MXj vector.
00627           //
00628           Teuchos::RefCountPtr<MV> oldMXj = MVT::CloneCopy( *MXj );
00629 
00630           if (numX > 0) {
00631             //
00632             // Apply the first step of Gram-Schmidt
00633             //
00634             MVT::MvTransMv( one, *prevXj, *MXj, product );
00635             //
00636             MVT::MvTimesMatAddMv( -one, *prevXj, product, one, *Xj );
00637             //
00638             if (M) {
00639               if (xc == mxc) {
00640                 Teuchos::RefCountPtr<MV> prevMXj = MVT::CloneView( MX, prev_idx );
00641                 MVT::MvTimesMatAddMv( -one, *prevMXj, product, one, *MXj );
00642               }
00643               else {
00644                 OPT::Apply( *M, *Xj, *MXj );
00645               }
00646             }
00647             //
00648             // Compute new M-norm
00649             //
00650             MVT::MvDot( *Xj, *MXj, &newDot );
00651             //
00652             // Check if a correction is needed.
00653             //
00654             if ( SCT::magnitude(kappa*newDot[0]) < SCT::magnitude(oldDot[0]) ) {
00655               //
00656               // Apply the second step of Gram-Schmidt
00657               //
00658               MVT::MvTransMv( one, *prevXj, *MXj, product );
00659               //
00660               MVT::MvTimesMatAddMv( -one, *prevXj, product, one, *Xj );
00661               //
00662               if (M) {
00663                 if (xc == mxc) {
00664                   Teuchos::RefCountPtr<MV> prevMXj = MVT::CloneView( MX, prev_idx );
00665                   MVT::MvTimesMatAddMv( -one, *prevMXj, product, one, *MXj );
00666                 }
00667                 else {
00668                   OPT::Apply( *M, *Xj, *MXj );
00669                 }
00670               }
00671             } // if (kappa*newDot[0] < oldDot[0])
00672             
00673           } // if (numX > 0)
00674             
00675           //
00676           // Compute M-norm with old MXj
00677           // NOTE:  Re-using newDot vector
00678           //
00679           MVT::MvDot( *Xj, *oldMXj, &newDot );
00680           
00681           // even if the vector Xj is complex, the inner product
00682           // Xj^H oldMXj should be, not only real, but positive.
00683           // If the real part isn't positive, then it suggests that 
00684           // Xj had very little energy outside of the previous vectors. Check this.
00685 
00686           if ( SCT::magnitude(newDot[0]) > SCT::magnitude(oldDot[0]*eps*eps) && SCT::real(newDot[0]) > zero ) {
00687             MVT::MvAddMv( one/Teuchos::ScalarTraits<ScalarType>::squareroot(newDot[0]), *Xj, zero, *Xj, *Xj );
00688             if (M) {
00689               MVT::MvAddMv( one/Teuchos::ScalarTraits<ScalarType>::squareroot(newDot[0]), *MXj, zero, *MXj, *MXj );
00690             }
00691             rankDef = false;
00692             break;
00693           }
00694           else {
00695             info += 1;
00696             MVT::MvRandom( *Xj );
00697             if (M) {
00698               OPT::Apply( *M, *Xj, *MXj );
00699             }
00700             if (orthoType == 0)
00701               massOrthonormalize(*Xj, *MXj, M, Q, 1, 1, kappa);
00702           } // if (norm > oldDot*eps*eps)
00703           
00704         }  // for (numTrials = 0; numTrials < 10; ++numTrials)
00705         
00706         if (rankDef == true) {
00707           MVT::MvInit( *Xj, zero );
00708           if (M)
00709             MVT::MvInit( *MXj, zero );
00710           info = -1;
00711           break;
00712         }
00713         
00714       } // for (j = 0; j < howMany; ++j)
00715         
00716     } // if (orthoType != 1)
00717     //timeNorm += MyWatch.WallTime();
00718 
00719     return info;
00720     
00721   }
00722   
00723   template<class ScalarType, class MV, class OP>
00724   int ModalSolverUtils<ScalarType, MV, OP>::directSolver(int size, const Teuchos::SerialDenseMatrix<int,ScalarType> &KK, 
00725                                                          const Teuchos::SerialDenseMatrix<int,ScalarType> *MM,
00726                                                          Teuchos::SerialDenseMatrix<int,ScalarType>* EV,
00727                                                          std::vector<MagnitudeType>* theta,
00728                                                          int* nev, int esType) const
00729   {
00730     // Routine for computing the first NEV generalized eigenpairs of the symmetric pencil (KK, MM)
00731     //
00732     // Parameter variables:
00733     //
00734     // size : Dimension of the eigenproblem (KK, MM)
00735     //
00736     // KK : Hermitian "stiffness" matrix 
00737     //
00738     // MM : Hermitian positive-definite "mass" matrix
00739     //
00740     // EV : Array to store the nev eigenvectors 
00741     //
00742     // theta : Array to store the eigenvalues (Size = nev )
00743     //
00744     // nev : Number of the smallest eigenvalues requested (input)
00745     //       Number of the smallest computed eigenvalues (output)
00746     //
00747     // esType : Flag to select the algorithm
00748     //
00749     // esType =  0 (default) Uses LAPACK routine (Cholesky factorization of MM)
00750     //                       with deflation of MM to get orthonormality of 
00751     //                       eigenvectors (S^T MM S = I)
00752     //
00753     // esType =  1           Uses LAPACK routine (Cholesky factorization of MM)
00754     //                       (no check of orthonormality)
00755     //
00756     // esType = 10           Uses LAPACK routine for simple eigenproblem on KK
00757     //                       (MM is not referenced in this case)
00758     //
00759     // Note: The code accesses only the upper triangular part of KK and MM.
00760     //
00761     // Return the integer info on the status of the computation
00762     //
00763     // info = 0 >> Success
00764     //
00765     // info = - 20 >> Failure in LAPACK routine
00766 
00767     // Define local arrays
00768 
00769     // Create blas/lapack objects.
00770     Teuchos::LAPACK<int,ScalarType> lapack;
00771     Teuchos::BLAS<int,ScalarType> blas;
00772 
00773     int i, j;
00774     int rank = 0;
00775     int info = 0;
00776 
00777     // Query LAPACK for the "optimal" block size for HEGV
00778     std::string lapack_name = "hetrd";
00779     std::string lapack_opts = "u";
00780     int NB = lapack.ILAENV(1, lapack_name, lapack_opts, size, -1, -1, -1);
00781     int lwork = size*(NB+1);
00782     std::vector<ScalarType> work(lwork);
00783     std::vector<MagnitudeType> rwork(3*size-2);
00784     // tt contains the eigenvalues from HEGV, which are necessarily real, and
00785     // HEGV expects this vector to be real as well
00786     std::vector<MagnitudeType> tt( size );
00787     typedef typename std::vector<MagnitudeType>::iterator MTIter;
00788 
00789     MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
00790     // MagnitudeType tol = 1e-12;
00791     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00792     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00793 
00794     Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
00795     Teuchos::RefCountPtr<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
00796 
00797     switch (esType) {
00798 
00799     default:
00800     case 0:
00801       //
00802       // Use LAPACK to compute the generalized eigenvectors
00803       //
00804       for (rank = size; rank > 0; --rank) {
00805         
00806         U = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );              
00807         //
00808         // Copy KK & MM
00809         //
00810         KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, rank, rank ) );
00811         MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
00812         //
00813         // Solve the generalized eigenproblem with LAPACK
00814         //
00815         info = 0;
00816         lapack.HEGV(1, 'V', 'U', rank, KKcopy->values(), KKcopy->stride(), 
00817                     MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
00818                     &rwork[0], &info);
00819         //
00820         // Treat error messages
00821         //
00822         if (info < 0) {
00823           cerr << endl;
00824           cerr << " In HEGV, argument " << -info << "has an illegal value.\n";
00825           cerr << endl;
00826           return -20;
00827         }
00828         if (info > 0) {
00829           if (info > rank)
00830             rank = info - rank;
00831           continue;
00832         }
00833         //
00834         // Check the quality of eigenvectors ( using mass-orthonormality )
00835         //
00836         MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
00837         for (i = 0; i < rank; ++i) {
00838           for (j = 0; j < i; ++j) {
00839             (*MMcopy)(i,j) = SCT::conjugate((*MM)(j,i));
00840           }
00841         }
00842         // U = 0*U + 1*MMcopy*KKcopy = MMcopy * KKcopy
00843         int ret = U->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*MMcopy,*KKcopy,zero);
00844         assert( ret == 0 );
00845         // MMcopy = 0*MMcopy + 1*KKcopy^H*U = KKcopy^H * MMcopy * KKcopy
00846         ret = MMcopy->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,one,*KKcopy,*U,zero);
00847         assert( ret == 0 );
00848         MagnitudeType maxNorm = SCT::magnitude(zero);
00849         MagnitudeType maxOrth = SCT::magnitude(zero);
00850         for (i = 0; i < rank; ++i) {
00851           for (j = i; j < rank; ++j) {
00852             if (j == i)
00853               maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
00854                       ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;            
00855             else 
00856               maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
00857                       ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
00858           }
00859         }
00860         /*        if (verbose > 4) {
00861         cout << " >> Local eigensolve >> Size: " << rank;
00862         cout.precision(2);
00863         cout.setf(ios::scientific, ios::floatfield);
00864         cout << " Normalization error: " << maxNorm;
00865         cout << " Orthogonality error: " << maxOrth;
00866         cout << endl;
00867         }*/
00868         if ((maxNorm <= tol) && (maxOrth <= tol)) {
00869           break;
00870         }
00871       } // for (rank = size; rank > 0; --rank)
00872       //
00873       // Copy the computed eigenvectors and eigenvalues
00874       // ( they may be less than the number requested because of deflation )
00875       //
00876       
00877       // cout << "directSolve    rank: " << rank << "\tsize: " << size << endl;
00878       
00879       *nev = (rank < *nev) ? rank : *nev;
00880       EV->putScalar( zero );
00881       std::copy<MTIter,MTIter>(tt.begin(),tt.begin()+*nev,theta->begin());
00882       for (i = 0; i < *nev; ++i) {
00883         blas.COPY( rank, (*KKcopy)[i], 1, (*EV)[i], 1 );
00884       }
00885       
00886       break;
00887 
00888     case 1:
00889       //
00890       // Use the Cholesky factorization of MM to compute the generalized eigenvectors
00891       //
00892       // Copy KK & MM
00893       //
00894       KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
00895       MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, size, size ) );
00896       //
00897       // Solve the generalized eigenproblem with LAPACK
00898       //
00899       info = 0;
00900       lapack.HEGV(1, 'V', 'U', size, KKcopy->values(), KKcopy->stride(), 
00901                   MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
00902                   &rwork[0], &info);
00903       //
00904       // Treat error messages
00905       //
00906       if (info < 0) {
00907         cerr << endl;
00908         cerr << " In HEGV, argument " << -info << "has an illegal value.\n";
00909         cerr << endl;
00910         return -20;
00911       }
00912       if (info > 0) {
00913         if (info > size)
00914           *nev = 0;
00915         else {
00916           cerr << endl;
00917           cerr << " In HEGV, DPOTRF or DHEEV returned an error code (" << info << ").\n";
00918           cerr << endl;
00919           return -20; 
00920         }
00921       }
00922       //
00923       // Copy the eigenvectors and eigenvalues
00924       //
00925       std::copy<MTIter,MTIter>(tt.begin(),tt.begin()+*nev,theta->begin());
00926       for (i = 0; i < *nev; ++i) {
00927         blas.COPY( size, (*KKcopy)[i], 1, (*EV)[i], 1 );
00928       }
00929       
00930       break;
00931 
00932     case 10:
00933       //
00934       // Simple eigenproblem
00935       //
00936       // Copy KK
00937       //
00938       KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
00939       //
00940       // Solve the generalized eigenproblem with LAPACK
00941       //
00942       lapack.HEEV('V', 'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
00943       //
00944       // Treat error messages
00945       if (info != 0) {
00946         cerr << endl;
00947         if (info < 0) 
00948           cerr << " In DHEEV, argument " << -info << " has an illegal value\n";
00949         else
00950           cerr << " In DHEEV, the algorithm failed to converge (" << info << ").\n";
00951         cerr << endl;
00952         info = -20;
00953         break;
00954       }
00955       //
00956       // Copy the eigenvectors
00957       //
00958       std::copy<MTIter,MTIter>(tt.begin(),tt.begin()+*nev,theta->begin());
00959       for (i = 0; i < *nev; ++i) {
00960         blas.COPY( size, (*KKcopy)[i], 1, (*EV)[i], 1 );
00961       }
00962       
00963       break;
00964       
00965     }
00966     
00967     // Clear the memory
00968     return info;
00969   }
00970 
00971 
00972   //-----------------------------------------------------------------------------
00973   // 
00974   //  SANITY CHECKING METHODS
00975   //
00976   //-----------------------------------------------------------------------------
00977 
00978   template<class ScalarType, class MV, class OP>
00979   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ModalSolverUtils<ScalarType, MV, OP>::errorOrthogonality(const MV *X, const MV *R, 
00980                                                                       const OP *M) const
00981   {
00982     // Return the maximum value of R_i^T * M * X_j / || MR_i || || X_j ||
00983     // When M is not specified, the identity is used.
00984     MagnitudeType maxDot = SCT::magnitude(SCT::zero());
00985     
00986     int xc = (X) ? MVT::GetNumberVecs( *X ) : 0;
00987     int rc = (R) ? MVT::GetNumberVecs( *R ) : 0;
00988     
00989     if (xc*rc == 0) {
00990       return maxDot;
00991     }
00992     
00993     int i, j;
00994     Teuchos::RefCountPtr<MV> MR;
00995     std::vector<MagnitudeType> normMR( rc );
00996     std::vector<MagnitudeType> normX( xc );
00997     if (M) {
00998       MR = MVT::Clone( *R, rc );
00999       OPT::Apply( *M, *R, *MR );
01000     }
01001     else {
01002       MR = MVT::CloneCopy( *R );
01003     }
01004     MVT::MvNorm( *MR, &normMR );
01005     MVT::MvNorm( *X, &normX );
01006 
01007     MagnitudeType dot = SCT::magnitude(SCT::zero());
01008     Teuchos::SerialDenseMatrix<int, ScalarType> xTMr( xc, rc );
01009     MVT::MvTransMv( 1.0, *X, *MR, xTMr );    
01010     for (i = 0; i < xc; ++i) {
01011       for (j = 0; j < rc; ++j) {
01012         dot = SCT::magnitude(xTMr(i,j)) / (normMR[j]*normX[i]);
01013         maxDot = (dot > maxDot) ? dot : maxDot;
01014       }
01015     }
01016     
01017     return maxDot;
01018   }
01019 
01020   template<class ScalarType, class MV, class OP>
01021   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ModalSolverUtils<ScalarType, MV, OP>::errorOrthonormality(const MV *X, const OP *M) const
01022   {
01023     // Return the maximum coefficient of the matrix X^T * M * X - I
01024     // When M is not specified, the identity is used.
01025     MagnitudeType maxDot = SCT::magnitude(SCT::zero());
01026     MagnitudeType one = SCT::magnitude(SCT::one());
01027     
01028     int xc = (X) ? MVT::GetNumberVecs( *X ) : 0;
01029     if (xc == 0) {
01030       return maxDot;
01031     }
01032     
01033     int i, j;
01034     std::vector<int> index( 1 );
01035     std::vector<ScalarType> dot( 1 );
01036     MagnitudeType tmpdot;
01037     Teuchos::RefCountPtr<MV> MXi;
01038     Teuchos::RefCountPtr<const MV> Xi;
01039 
01040     // Create space if there is a M matrix specified.
01041     if (M) {
01042       MXi = MVT::Clone( *X, 1 );
01043     }
01044 
01045     for (i = 0; i < xc; ++i) {
01046       index[0] = i;
01047       if (M) {
01048         Xi = MVT::CloneView( *X, index );
01049         OPT::Apply( *M, *Xi, *MXi );
01050       }
01051       else {
01052         MXi = MVT::CloneView( *(const_cast<MV *>(X)), index );
01053       }
01054       for (j = 0; j < xc; ++j) {
01055         index[0] = j;
01056         Xi = MVT::CloneView( *X, index );
01057         MVT::MvDot( *Xi, *MXi, &dot );
01058         tmpdot = (i == j) ? SCT::magnitude(dot[0] - one) : SCT::magnitude(dot[0]);
01059         maxDot = (tmpdot > maxDot) ? tmpdot : maxDot;
01060       }
01061     }
01062     
01063     return maxDot;    
01064   }
01065 
01066   template<class ScalarType, class MV, class OP>
01067   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType ModalSolverUtils<ScalarType, MV, OP>::errorEquality(const MV *X, const MV *MX, 
01068                                                                  const OP *M) const
01069   {
01070     // Return the maximum coefficient of the matrix M * X - MX
01071     // scaled by the maximum coefficient of MX.
01072     // When M is not specified, the identity is used.
01073     
01074     MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
01075     
01076     int xc = (X) ? MVT::GetNumberVecs( *X ) : 0;
01077     int mxc = (MX) ? MVT::GetNumberVecs( *MX ) : 0;
01078     
01079     if ((xc != mxc) || (xc*mxc == 0)) {
01080       return maxDiff;
01081     }
01082     
01083     int i;
01084     MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
01085     std::vector<MagnitudeType> tmp( xc );
01086     MVT::MvNorm( *MX, &tmp );
01087 
01088     for (i = 0; i < xc; ++i) {
01089       maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
01090     }
01091 
01092     std::vector<int> index( 1 );
01093     Teuchos::RefCountPtr<MV> MtimesX; 
01094     if (M) {
01095       MtimesX = MVT::Clone( *X, xc );
01096       OPT::Apply( *M, *X, *MtimesX );
01097     }
01098     else {
01099       MtimesX = MVT::CloneCopy( *(const_cast<MV *>(X)) );
01100     }
01101     MVT::MvAddMv( -1.0, *MX, 1.0, *MtimesX, *MtimesX );
01102     MVT::MvNorm( *MtimesX, &tmp );
01103    
01104     for (i = 0; i < xc; ++i) {
01105       maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
01106     }
01107     
01108     return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
01109     
01110   }    
01111   
01112 } // end namespace Anasazi
01113 
01114 #endif // ANASAZI_MODAL_SOLVER_UTILS_HPP
01115 

Generated on Thu Sep 18 12:31:37 2008 for Anasazi by doxygen 1.3.9.1