AnasaziSolverUtils.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_SOLVER_UTILS_HPP
00030 #define ANASAZI_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 SolverUtils 
00062   {  
00063   public:
00064     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00065     typedef typename Teuchos::ScalarTraits<ScalarType>  SCT;
00066     
00068 
00069 
00071     SolverUtils();
00072     
00074     virtual ~SolverUtils() {};
00075     
00077     
00079 
00080     
00082     static void permuteVectors(const int n, const std::vector<int> &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids = 0);
00083 
00085     static void permuteVectors(const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q);
00086 
00088 
00090 
00091 
00093 
00114     static void applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix<int,ScalarType> &H, const std::vector<ScalarType> &tau, Teuchos::RCP<MV> workMV = Teuchos::null);
00115 
00117 
00119 
00120 
00122 
00148     static int directSolver(int size, const Teuchos::SerialDenseMatrix<int,ScalarType> &KK, 
00149                      Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
00150                      Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
00151                      std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
00152                      int &nev, int esType = 0);
00154 
00156 
00157 
00159 
00161     static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M = Teuchos::null);
00162     
00164     
00165   private:
00166 
00168 
00169 
00170     typedef MultiVecTraits<ScalarType,MV> MVT;
00171     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00172 
00174   };
00175 
00176   //-----------------------------------------------------------------------------
00177   // 
00178   //  CONSTRUCTOR
00179   //
00180   //-----------------------------------------------------------------------------  
00181 
00182   template<class ScalarType, class MV, class OP>
00183   SolverUtils<ScalarType, MV, OP>::SolverUtils() {}
00184 
00185 
00186   //-----------------------------------------------------------------------------
00187   // 
00188   //  SORTING METHODS
00189   //
00190   //-----------------------------------------------------------------------------
00191   
00193   // permuteVectors for MV
00194   template<class ScalarType, class MV, class OP>
00195   void SolverUtils<ScalarType, MV, OP>::permuteVectors(
00196               const int n,
00197               const std::vector<int> &perm, 
00198               MV &Q, 
00199               std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids)
00200   {
00201     // Permute the vectors according to the permutation vector \c perm, and
00202     // optionally the residual vector \c resids
00203     
00204     int i, j;
00205     std::vector<int> permcopy(perm), swapvec(n-1);
00206     std::vector<int> index(1);
00207     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00208     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00209 
00210     TEST_FOR_EXCEPTION(n > MVT::GetNumberVecs(Q), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): argument n larger than width of input multivector.");
00211 
00212     // We want to recover the elementary permutations (individual swaps) 
00213     // from the permutation vector. Do this by constructing the inverse
00214     // of the permutation, by sorting them to {1,2,...,n}, and recording
00215     // the elementary permutations of the inverse.
00216     for (i=0; i<n-1; i++) {
00217       //
00218       // find i in the permcopy vector
00219       for (j=i; j<n; j++) {
00220         if (permcopy[j] == i) {
00221           // found it at index j
00222           break;
00223         }
00224         TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): permutation index invalid.");
00225       }
00226       //
00227       // Swap two scalars
00228       std::swap( permcopy[j], permcopy[i] );
00229 
00230       swapvec[i] = j;
00231     }
00232       
00233     // now apply the elementary permutations of the inverse in reverse order
00234     for (i=n-2; i>=0; i--) {
00235       j = swapvec[i];
00236       //
00237       // Swap (i,j)
00238       //
00239       // Swap residuals (if they exist)
00240       if (resids) {
00241         std::swap(  (*resids)[i], (*resids)[j] );
00242       }
00243       //
00244       // Swap corresponding vectors
00245       index[0] = j;
00246       Teuchos::RCP<MV> tmpQ = MVT::CloneCopy( Q, index );
00247       Teuchos::RCP<MV> tmpQj = MVT::CloneViewNonConst( Q, index );
00248       index[0] = i;
00249       Teuchos::RCP<MV> tmpQi = MVT::CloneViewNonConst( Q, index );
00250       MVT::MvAddMv( one, *tmpQi, zero, *tmpQi, *tmpQj );
00251       MVT::MvAddMv( one, *tmpQ, zero, *tmpQ, *tmpQi );
00252     }
00253   }
00254 
00255 
00257   // permuteVectors for MV
00258   template<class ScalarType, class MV, class OP>
00259   void SolverUtils<ScalarType, MV, OP>::permuteVectors(
00260               const std::vector<int> &perm, 
00261               Teuchos::SerialDenseMatrix<int,ScalarType> &Q)
00262   {
00263     // Permute the vectors in Q according to the permutation vector \c perm, and
00264     // optionally the residual vector \c resids
00265     Teuchos::BLAS<int,ScalarType> blas;
00266     const int n = perm.size();
00267     const int m = Q.numRows();
00268     
00269     TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
00270 
00271     // Sort the primitive ritz vectors
00272     Teuchos::SerialDenseMatrix<int,ScalarType> copyQ( Q );
00273     for (int i=0; i<n; i++) {
00274       blas.COPY(m, copyQ[perm[i]], 1, Q[i], 1);
00275     }
00276   }
00277 
00278 
00279   //-----------------------------------------------------------------------------
00280   // 
00281   //  BASIS UPDATE METHODS
00282   //
00283   //-----------------------------------------------------------------------------
00284 
00285   // apply householder reflectors to multivector
00286   template<class ScalarType, class MV, class OP>
00287   void SolverUtils<ScalarType, MV, OP>::applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix<int,ScalarType> &H, const std::vector<ScalarType> &tau, Teuchos::RCP<MV> workMV) {
00288 
00289     const int n = MVT::GetNumberVecs(V);
00290     const ScalarType ONE = SCT::one();
00291     const ScalarType ZERO = SCT::zero();
00292 
00293     // early exit if V has zero-size or if k==0
00294     if (MVT::GetNumberVecs(V) == 0 || MVT::GetVecLength(V) == 0 || k == 0) {
00295       return;
00296     }
00297 
00298     if (workMV == Teuchos::null) {
00299       // user did not give us any workspace; allocate some
00300       workMV = MVT::Clone(V,1);
00301     }
00302     else if (MVT::GetNumberVecs(*workMV) > 1) {
00303       std::vector<int> first(1);
00304       first[0] = 0;
00305       workMV = MVT::CloneViewNonConst(*workMV,first);
00306     }
00307     else {
00308       TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
00309     }
00310     // Q = H_1 ... H_k is square, with as many rows as V has vectors
00311     // however, H need only have k columns, one each for the k reflectors.
00312     TEST_FOR_EXCEPTION( H.numCols() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): H must have at least k columns.");
00313     TEST_FOR_EXCEPTION( (int)tau.size() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): tau must have at least k entries.");
00314     TEST_FOR_EXCEPTION( H.numRows() != MVT::GetNumberVecs(V), std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): Size of H,V are inconsistent.");
00315 
00316     // perform the loop
00317     // flops: Sum_{i=0:k-1} 4 m (n-i) == 4mnk - 2m(k^2- k)
00318     for (int i=0; i<k; i++) {
00319       // apply V H_i+1 = V - tau_i+1 (V v_i+1) v_i+1^T
00320       // because of the structure of v_i+1, this transform does not affect the first i columns of V
00321       std::vector<int> activeind(n-i);
00322       for (int j=0; j<n-i; j++) activeind[j] = j+i;
00323       Teuchos::RCP<MV> actV = MVT::CloneViewNonConst(V,activeind);
00324 
00325       // note, below H_i, v_i and tau_i are mathematical objects which use 1-based indexing
00326       // while H, v and tau are data structures using 0-based indexing
00327 
00328       // get v_i+1: i-th column of H
00329       Teuchos::SerialDenseMatrix<int,ScalarType> v(Teuchos::Copy,H,n-i,1,i,i);
00330       // v_i+1(1:i) = 0: this isn't part of v
00331       // e_i+1^T v_i+1 = 1 = v(0)
00332       v(0,0) = ONE;
00333 
00334       // compute -tau_i V v_i
00335       // tau_i+1 is tau[i]
00336       // flops: 2 m n-i
00337       MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
00338 
00339       // perform V = V + workMV v_i^T
00340       // flops: 2 m n-i
00341       Teuchos::SerialDenseMatrix<int,ScalarType> vT(v,Teuchos::CONJ_TRANS);
00342       MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
00343 
00344       actV = Teuchos::null;
00345     }
00346   }
00347 
00348   
00349   //-----------------------------------------------------------------------------
00350   // 
00351   //  EIGENSOLVER PROJECTION METHODS
00352   //
00353   //-----------------------------------------------------------------------------
00354   
00355   template<class ScalarType, class MV, class OP>
00356   int SolverUtils<ScalarType, MV, OP>::directSolver(
00357       int size, 
00358       const Teuchos::SerialDenseMatrix<int,ScalarType> &KK, 
00359       Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
00360       Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
00361       std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
00362       int &nev, int esType)
00363   {
00364     // Routine for computing the first NEV generalized eigenpairs of the symmetric pencil (KK, MM)
00365     //
00366     // Parameter variables:
00367     //
00368     // size : Dimension of the eigenproblem (KK, MM)
00369     //
00370     // KK : Hermitian "stiffness" matrix 
00371     //
00372     // MM : Hermitian positive-definite "mass" matrix
00373     //
00374     // EV : Matrix to store the nev eigenvectors 
00375     //
00376     // theta : Array to store the eigenvalues (Size = nev )
00377     //
00378     // nev : Number of the smallest eigenvalues requested (input)
00379     //       Number of the smallest computed eigenvalues (output)
00380     //       Routine may compute and return more or less eigenvalues than requested.
00381     //
00382     // esType : Flag to select the algorithm
00383     //
00384     // esType =  0 (default) Uses LAPACK routine (Cholesky factorization of MM)
00385     //                       with deflation of MM to get orthonormality of 
00386     //                       eigenvectors (S^T MM S = I)
00387     //
00388     // esType =  1           Uses LAPACK routine (Cholesky factorization of MM)
00389     //                       (no check of orthonormality)
00390     //
00391     // esType = 10           Uses LAPACK routine for simple eigenproblem on KK
00392     //                       (MM is not referenced in this case)
00393     //
00394     // Note: The code accesses only the upper triangular part of KK and MM.
00395     //
00396     // Return the integer info on the status of the computation
00397     //
00398     // info = 0 >> Success
00399     //
00400     // info < 0 >> error in the info-th argument
00401     // info = - 20 >> Failure in LAPACK routine
00402 
00403     // Define local arrays
00404 
00405     // Create blas/lapack objects.
00406     Teuchos::LAPACK<int,ScalarType> lapack;
00407     Teuchos::BLAS<int,ScalarType> blas;
00408 
00409     int rank = 0;
00410     int info = 0;
00411     
00412     if (size < nev || size < 0) {
00413       return -1;
00414     }
00415     if (KK.numCols() < size || KK.numRows() < size) {
00416       return -2;
00417     }
00418     if ((esType == 0 || esType == 1)) {
00419       if (MM == Teuchos::null) {
00420         return -3;
00421       }
00422       else if (MM->numCols() < size || MM->numRows() < size) {
00423         return -3;
00424       }
00425     }
00426     if (EV.numCols() < size || EV.numRows() < size) {
00427       return -4;
00428     }
00429     if (theta.size() < (unsigned int) size) {
00430       return -5;
00431     }
00432     if (nev <= 0) {
00433       return -6;
00434     }
00435 
00436     // Query LAPACK for the "optimal" block size for HEGV
00437     std::string lapack_name = "hetrd";
00438     std::string lapack_opts = "u";
00439     int NB = lapack.ILAENV(1, lapack_name, lapack_opts, size, -1, -1, -1);
00440     int lwork = size*(NB+1);
00441     std::vector<ScalarType> work(lwork);
00442     std::vector<MagnitudeType> rwork(3*size-2);
00443     // tt contains the eigenvalues from HEGV, which are necessarily real, and
00444     // HEGV expects this vector to be real as well
00445     std::vector<MagnitudeType> tt( size );
00446     typedef typename std::vector<MagnitudeType>::iterator MTIter;
00447 
00448     MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
00449     // MagnitudeType tol = 1e-12;
00450     ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00451     ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00452 
00453     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
00454     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
00455 
00456     switch (esType) {
00457       default:
00458       case 0:
00459         //
00460         // Use LAPACK to compute the generalized eigenvectors
00461         //
00462         for (rank = size; rank > 0; --rank) {
00463 
00464           U = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );              
00465           //
00466           // Copy KK & MM
00467           //
00468           KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, rank, rank ) );
00469           MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
00470           //
00471           // Solve the generalized eigenproblem with LAPACK
00472           //
00473           info = 0;
00474           lapack.HEGV(1, 'V', 'U', rank, KKcopy->values(), KKcopy->stride(), 
00475               MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
00476               &rwork[0], &info);
00477           //
00478           // Treat error messages
00479           //
00480           if (info < 0) {
00481             std::cerr << std::endl;
00482             std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
00483             std::cerr << std::endl;
00484             return -20;
00485           }
00486           if (info > 0) {
00487             if (info > rank)
00488               rank = info - rank;
00489             continue;
00490           }
00491           //
00492           // Check the quality of eigenvectors ( using mass-orthonormality )
00493           //
00494           MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
00495           for (int i = 0; i < rank; ++i) {
00496             for (int j = 0; j < i; ++j) {
00497               (*MMcopy)(i,j) = SCT::conjugate((*MM)(j,i));
00498             }
00499           }
00500           // U = 0*U + 1*MMcopy*KKcopy = MMcopy * KKcopy
00501           TEST_FOR_EXCEPTION( 
00502               U->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*MMcopy,*KKcopy,zero) != 0,
00503               std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
00504           // MMcopy = 0*MMcopy + 1*KKcopy^H*U = KKcopy^H * MMcopy * KKcopy
00505           TEST_FOR_EXCEPTION( 
00506               MMcopy->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,one,*KKcopy,*U,zero) != 0,
00507               std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
00508           MagnitudeType maxNorm = SCT::magnitude(zero);
00509           MagnitudeType maxOrth = SCT::magnitude(zero);
00510           for (int i = 0; i < rank; ++i) {
00511             for (int j = i; j < rank; ++j) {
00512               if (j == i)
00513                 maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
00514                   ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;            
00515               else 
00516                 maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
00517                   ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
00518             }
00519           }
00520           /*        if (verbose > 4) {
00521                     std::cout << " >> Local eigensolve >> Size: " << rank;
00522                     std::cout.precision(2);
00523                     std::cout.setf(std::ios::scientific, std::ios::floatfield);
00524                     std::cout << " Normalization error: " << maxNorm;
00525                     std::cout << " Orthogonality error: " << maxOrth;
00526                     std::cout << endl;
00527                     }*/
00528           if ((maxNorm <= tol) && (maxOrth <= tol)) {
00529             break;
00530           }
00531         } // for (rank = size; rank > 0; --rank)
00532         //
00533         // Copy the computed eigenvectors and eigenvalues
00534         // ( they may be less than the number requested because of deflation )
00535         //
00536         // std::cout << "directSolve    rank: " << rank << "\tsize: " << size << endl;
00537         nev = (rank < nev) ? rank : nev;
00538         EV.putScalar( zero );
00539         std::copy(tt.begin(),tt.begin()+nev,theta.begin());
00540         for (int i = 0; i < nev; ++i) {
00541           blas.COPY( rank, (*KKcopy)[i], 1, EV[i], 1 );
00542         }
00543         break;
00544 
00545       case 1:
00546         //
00547         // Use the Cholesky factorization of MM to compute the generalized eigenvectors
00548         //
00549         // Copy KK & MM
00550         //
00551         KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
00552         MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, size, size ) );
00553         //
00554         // Solve the generalized eigenproblem with LAPACK
00555         //
00556         info = 0;
00557         lapack.HEGV(1, 'V', 'U', size, KKcopy->values(), KKcopy->stride(), 
00558             MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
00559             &rwork[0], &info);
00560         //
00561         // Treat error messages
00562         //
00563         if (info < 0) {
00564           std::cerr << std::endl;
00565           std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
00566           std::cerr << std::endl;
00567           return -20;
00568         }
00569         if (info > 0) {
00570           if (info > size)
00571             nev = 0;
00572           else {
00573             std::cerr << std::endl;
00574             std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, DPOTRF or DHEEV returned an error code (" << info << ").\n";
00575             std::cerr << std::endl;
00576             return -20; 
00577           }
00578         }
00579         //
00580         // Copy the eigenvectors and eigenvalues
00581         //
00582         std::copy(tt.begin(),tt.begin()+nev,theta.begin());
00583         for (int i = 0; i < nev; ++i) {
00584           blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
00585         }
00586         break;
00587 
00588       case 10:
00589         //
00590         // Simple eigenproblem
00591         //
00592         // Copy KK
00593         //
00594         KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
00595         //
00596         // Solve the generalized eigenproblem with LAPACK
00597         //
00598         lapack.HEEV('V', 'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
00599         //
00600         // Treat error messages
00601         if (info != 0) {
00602           std::cerr << std::endl;
00603           if (info < 0) {
00604             std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, argument " << -info << " has an illegal value\n";
00605           }
00606           else {
00607             std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, the algorithm failed to converge (" << info << ").\n";
00608           }
00609           std::cerr << std::endl;
00610           info = -20;
00611           break;
00612         }
00613         //
00614         // Copy the eigenvectors
00615         //
00616         std::copy(tt.begin(),tt.begin()+nev,theta.begin());
00617         for (int i = 0; i < nev; ++i) {
00618           blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
00619         }
00620         break;
00621     }
00622 
00623     return info;
00624   }
00625 
00626 
00627   //-----------------------------------------------------------------------------
00628   // 
00629   //  SANITY CHECKING METHODS
00630   //
00631   //-----------------------------------------------------------------------------
00632 
00633   template<class ScalarType, class MV, class OP>
00634   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00635   SolverUtils<ScalarType, MV, OP>::errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M)
00636   {
00637     // Return the maximum coefficient of the matrix M * X - MX
00638     // scaled by the maximum coefficient of MX.
00639     // When M is not specified, the identity is used.
00640     
00641     MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
00642     
00643     int xc = MVT::GetNumberVecs(X);
00644     int mxc = MVT::GetNumberVecs(MX);
00645     
00646     TEST_FOR_EXCEPTION(xc != mxc,std::invalid_argument,"Anasazi::SolverUtils::errorEquality(): input multivecs have different number of columns.");
00647     if (xc == 0) {
00648       return maxDiff;
00649     }
00650     
00651     MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
00652     std::vector<MagnitudeType> tmp( xc );
00653     MVT::MvNorm(MX, tmp);
00654 
00655     for (int i = 0; i < xc; ++i) {
00656       maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
00657     }
00658 
00659     std::vector<int> index( 1 );
00660     Teuchos::RCP<MV> MtimesX; 
00661     if (M != Teuchos::null) {
00662       MtimesX = MVT::Clone( X, xc );
00663       OPT::Apply( *M, X, *MtimesX );
00664     }
00665     else {
00666       MtimesX = MVT::CloneCopy(X);
00667     }
00668     MVT::MvAddMv( -1.0, MX, 1.0, *MtimesX, *MtimesX );
00669     MVT::MvNorm( *MtimesX, tmp );
00670    
00671     for (int i = 0; i < xc; ++i) {
00672       maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
00673     }
00674     
00675     return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
00676     
00677   }    
00678   
00679 } // end namespace Anasazi
00680 
00681 #endif // ANASAZI_SOLVER_UTILS_HPP
00682 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 09:56:59 2011 for Anasazi by  doxygen 1.6.3