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