Anasazi  Version of the Day
AnasaziSolverUtils.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef ANASAZI_SOLVER_UTILS_HPP
43 #define ANASAZI_SOLVER_UTILS_HPP
44 
61 #include "AnasaziConfigDefs.hpp"
64 #include "Teuchos_ScalarTraits.hpp"
65 
66 #include "AnasaziOutputManager.hpp"
67 #include "Teuchos_BLAS.hpp"
68 #include "Teuchos_LAPACK.hpp"
69 #include "Teuchos_SerialDenseMatrix.hpp"
70 
71 namespace Anasazi {
72 
73  template<class ScalarType, class MV, class OP>
75  {
76  public:
77  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
78  typedef typename Teuchos::ScalarTraits<ScalarType> SCT;
79 
81 
82 
84  SolverUtils();
85 
87  virtual ~SolverUtils() {};
88 
90 
92 
93 
95  static void permuteVectors(const int n, const std::vector<int> &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids = 0);
96 
98  static void permuteVectors(const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q);
99 
101 
103 
104 
106 
127  static void applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix<int,ScalarType> &H, const std::vector<ScalarType> &tau, Teuchos::RCP<MV> workMV = Teuchos::null);
128 
130 
132 
133 
135 
161  static int directSolver(int size, const Teuchos::SerialDenseMatrix<int,ScalarType> &KK,
162  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
163  Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
164  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
165  int &nev, int esType = 0);
167 
169 
170 
172 
174  static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M = Teuchos::null);
175 
177 
178  private:
179 
181 
182 
183  typedef MultiVecTraits<ScalarType,MV> MVT;
185 
187  };
188 
189  //-----------------------------------------------------------------------------
190  //
191  // CONSTRUCTOR
192  //
193  //-----------------------------------------------------------------------------
194 
195  template<class ScalarType, class MV, class OP>
197 
198 
199  //-----------------------------------------------------------------------------
200  //
201  // SORTING METHODS
202  //
203  //-----------------------------------------------------------------------------
204 
206  // permuteVectors for MV
207  template<class ScalarType, class MV, class OP>
209  const int n,
210  const std::vector<int> &perm,
211  MV &Q,
212  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids)
213  {
214  // Permute the vectors according to the permutation vector \c perm, and
215  // optionally the residual vector \c resids
216 
217  int i, j;
218  std::vector<int> permcopy(perm), swapvec(n-1);
219  std::vector<int> index(1);
220  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
221  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
222 
223  TEUCHOS_TEST_FOR_EXCEPTION(n > MVT::GetNumberVecs(Q), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): argument n larger than width of input multivector.");
224 
225  // We want to recover the elementary permutations (individual swaps)
226  // from the permutation vector. Do this by constructing the inverse
227  // of the permutation, by sorting them to {1,2,...,n}, and recording
228  // the elementary permutations of the inverse.
229  for (i=0; i<n-1; i++) {
230  //
231  // find i in the permcopy vector
232  for (j=i; j<n; j++) {
233  if (permcopy[j] == i) {
234  // found it at index j
235  break;
236  }
237  TEUCHOS_TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): permutation index invalid.");
238  }
239  //
240  // Swap two scalars
241  std::swap( permcopy[j], permcopy[i] );
242 
243  swapvec[i] = j;
244  }
245 
246  // now apply the elementary permutations of the inverse in reverse order
247  for (i=n-2; i>=0; i--) {
248  j = swapvec[i];
249  //
250  // Swap (i,j)
251  //
252  // Swap residuals (if they exist)
253  if (resids) {
254  std::swap( (*resids)[i], (*resids)[j] );
255  }
256  //
257  // Swap corresponding vectors
258  index[0] = j;
259  Teuchos::RCP<MV> tmpQ = MVT::CloneCopy( Q, index );
260  Teuchos::RCP<MV> tmpQj = MVT::CloneViewNonConst( Q, index );
261  index[0] = i;
262  Teuchos::RCP<MV> tmpQi = MVT::CloneViewNonConst( Q, index );
263  MVT::MvAddMv( one, *tmpQi, zero, *tmpQi, *tmpQj );
264  MVT::MvAddMv( one, *tmpQ, zero, *tmpQ, *tmpQi );
265  }
266  }
267 
268 
270  // permuteVectors for MV
271  template<class ScalarType, class MV, class OP>
273  const std::vector<int> &perm,
274  Teuchos::SerialDenseMatrix<int,ScalarType> &Q)
275  {
276  // Permute the vectors in Q according to the permutation vector \c perm, and
277  // optionally the residual vector \c resids
278  Teuchos::BLAS<int,ScalarType> blas;
279  const int n = perm.size();
280  const int m = Q.numRows();
281 
282  TEUCHOS_TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
283 
284  // Sort the primitive ritz vectors
285  Teuchos::SerialDenseMatrix<int,ScalarType> copyQ(Teuchos::Copy, Q);
286  for (int i=0; i<n; i++) {
287  blas.COPY(m, copyQ[perm[i]], 1, Q[i], 1);
288  }
289  }
290 
291 
292  //-----------------------------------------------------------------------------
293  //
294  // BASIS UPDATE METHODS
295  //
296  //-----------------------------------------------------------------------------
297 
298  // apply householder reflectors to multivector
299  template<class ScalarType, class MV, class OP>
300  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) {
301 
302  const int n = MVT::GetNumberVecs(V);
303  const ScalarType ONE = SCT::one();
304  const ScalarType ZERO = SCT::zero();
305 
306  // early exit if V has zero-size or if k==0
307  if (MVT::GetNumberVecs(V) == 0 || MVT::GetGlobalLength(V) == 0 || k == 0) {
308  return;
309  }
310 
311  if (workMV == Teuchos::null) {
312  // user did not give us any workspace; allocate some
313  workMV = MVT::Clone(V,1);
314  }
315  else if (MVT::GetNumberVecs(*workMV) > 1) {
316  std::vector<int> first(1);
317  first[0] = 0;
318  workMV = MVT::CloneViewNonConst(*workMV,first);
319  }
320  else {
321  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
322  }
323  // Q = H_1 ... H_k is square, with as many rows as V has vectors
324  // however, H need only have k columns, one each for the k reflectors.
325  TEUCHOS_TEST_FOR_EXCEPTION( H.numCols() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): H must have at least k columns.");
326  TEUCHOS_TEST_FOR_EXCEPTION( (int)tau.size() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): tau must have at least k entries.");
327  TEUCHOS_TEST_FOR_EXCEPTION( H.numRows() != MVT::GetNumberVecs(V), std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): Size of H,V are inconsistent.");
328 
329  // perform the loop
330  // flops: Sum_{i=0:k-1} 4 m (n-i) == 4mnk - 2m(k^2- k)
331  for (int i=0; i<k; i++) {
332  // apply V H_i+1 = V - tau_i+1 (V v_i+1) v_i+1^T
333  // because of the structure of v_i+1, this transform does not affect the first i columns of V
334  std::vector<int> activeind(n-i);
335  for (int j=0; j<n-i; j++) activeind[j] = j+i;
336  Teuchos::RCP<MV> actV = MVT::CloneViewNonConst(V,activeind);
337 
338  // note, below H_i, v_i and tau_i are mathematical objects which use 1-based indexing
339  // while H, v and tau are data structures using 0-based indexing
340 
341  // get v_i+1: i-th column of H
342  Teuchos::SerialDenseMatrix<int,ScalarType> v(Teuchos::Copy,H,n-i,1,i,i);
343  // v_i+1(1:i) = 0: this isn't part of v
344  // e_i+1^T v_i+1 = 1 = v(0)
345  v(0,0) = ONE;
346 
347  // compute -tau_i V v_i
348  // tau_i+1 is tau[i]
349  // flops: 2 m n-i
350  MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
351 
352  // perform V = V + workMV v_i^T
353  // flops: 2 m n-i
354  Teuchos::SerialDenseMatrix<int,ScalarType> vT(v,Teuchos::CONJ_TRANS);
355  MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
356 
357  actV = Teuchos::null;
358  }
359  }
360 
361 
362  //-----------------------------------------------------------------------------
363  //
364  // EIGENSOLVER PROJECTION METHODS
365  //
366  //-----------------------------------------------------------------------------
367 
368  template<class ScalarType, class MV, class OP>
370  int size,
371  const Teuchos::SerialDenseMatrix<int,ScalarType> &KK,
372  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
373  Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
374  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
375  int &nev, int esType)
376  {
377  // Routine for computing the first NEV generalized eigenpairs of the symmetric pencil (KK, MM)
378  //
379  // Parameter variables:
380  //
381  // size : Dimension of the eigenproblem (KK, MM)
382  //
383  // KK : Hermitian "stiffness" matrix
384  //
385  // MM : Hermitian positive-definite "mass" matrix
386  //
387  // EV : Matrix to store the nev eigenvectors
388  //
389  // theta : Array to store the eigenvalues (Size = nev )
390  //
391  // nev : Number of the smallest eigenvalues requested (input)
392  // Number of the smallest computed eigenvalues (output)
393  // Routine may compute and return more or less eigenvalues than requested.
394  //
395  // esType : Flag to select the algorithm
396  //
397  // esType = 0 (default) Uses LAPACK routine (Cholesky factorization of MM)
398  // with deflation of MM to get orthonormality of
399  // eigenvectors (S^T MM S = I)
400  //
401  // esType = 1 Uses LAPACK routine (Cholesky factorization of MM)
402  // (no check of orthonormality)
403  //
404  // esType = 10 Uses LAPACK routine for simple eigenproblem on KK
405  // (MM is not referenced in this case)
406  //
407  // Note: The code accesses only the upper triangular part of KK and MM.
408  //
409  // Return the integer info on the status of the computation
410  //
411  // info = 0 >> Success
412  //
413  // info < 0 >> error in the info-th argument
414  // info = - 20 >> Failure in LAPACK routine
415 
416  // Define local arrays
417 
418  // Create blas/lapack objects.
419  Teuchos::LAPACK<int,ScalarType> lapack;
420  Teuchos::BLAS<int,ScalarType> blas;
421 
422  int rank = 0;
423  int info = 0;
424 
425  if (size < nev || size < 0) {
426  return -1;
427  }
428  if (KK.numCols() < size || KK.numRows() < size) {
429  return -2;
430  }
431  if ((esType == 0 || esType == 1)) {
432  if (MM == Teuchos::null) {
433  return -3;
434  }
435  else if (MM->numCols() < size || MM->numRows() < size) {
436  return -3;
437  }
438  }
439  if (EV.numCols() < size || EV.numRows() < size) {
440  return -4;
441  }
442  if (theta.size() < (unsigned int) size) {
443  return -5;
444  }
445  if (nev <= 0) {
446  return -6;
447  }
448 
449  // Query LAPACK for the "optimal" block size for HEGV
450  std::string lapack_name = "hetrd";
451  std::string lapack_opts = "u";
452  int NB = lapack.ILAENV(1, lapack_name, lapack_opts, size, -1, -1, -1);
453  int lwork = size*(NB+2); // For HEEV, lwork should be NB+2, instead of NB+1
454  std::vector<ScalarType> work(lwork);
455  std::vector<MagnitudeType> rwork(3*size-2);
456  // tt contains the eigenvalues from HEGV, which are necessarily real, and
457  // HEGV expects this vector to be real as well
458  std::vector<MagnitudeType> tt( size );
459  //typedef typename std::vector<MagnitudeType>::iterator MTIter; // unused
460 
461  MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
462  // MagnitudeType tol = 1e-12;
463  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
464  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
465 
466  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
467  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
468 
469  switch (esType) {
470  default:
471  case 0:
472  //
473  // Use LAPACK to compute the generalized eigenvectors
474  //
475  for (rank = size; rank > 0; --rank) {
476 
477  U = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );
478  //
479  // Copy KK & MM
480  //
481  KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, rank, rank ) );
482  MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
483  //
484  // Solve the generalized eigenproblem with LAPACK
485  //
486  info = 0;
487  lapack.HEGV(1, 'V', 'U', rank, KKcopy->values(), KKcopy->stride(),
488  MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
489  &rwork[0], &info);
490  //
491  // Treat error messages
492  //
493  if (info < 0) {
494  std::cerr << std::endl;
495  std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
496  std::cerr << std::endl;
497  return -20;
498  }
499  if (info > 0) {
500  if (info > rank)
501  rank = info - rank;
502  continue;
503  }
504  //
505  // Check the quality of eigenvectors ( using mass-orthonormality )
506  //
507  MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
508  for (int i = 0; i < rank; ++i) {
509  for (int j = 0; j < i; ++j) {
510  (*MMcopy)(i,j) = SCT::conjugate((*MM)(j,i));
511  }
512  }
513  // U = 0*U + 1*MMcopy*KKcopy = MMcopy * KKcopy
514  TEUCHOS_TEST_FOR_EXCEPTION(
515  U->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*MMcopy,*KKcopy,zero) != 0,
516  std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
517  // MMcopy = 0*MMcopy + 1*KKcopy^H*U = KKcopy^H * MMcopy * KKcopy
518  TEUCHOS_TEST_FOR_EXCEPTION(
519  MMcopy->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,one,*KKcopy,*U,zero) != 0,
520  std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
521  MagnitudeType maxNorm = SCT::magnitude(zero);
522  MagnitudeType maxOrth = SCT::magnitude(zero);
523  for (int i = 0; i < rank; ++i) {
524  for (int j = i; j < rank; ++j) {
525  if (j == i)
526  maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
527  ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;
528  else
529  maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
530  ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
531  }
532  }
533  /* if (verbose > 4) {
534  std::cout << " >> Local eigensolve >> Size: " << rank;
535  std::cout.precision(2);
536  std::cout.setf(std::ios::scientific, std::ios::floatfield);
537  std::cout << " Normalization error: " << maxNorm;
538  std::cout << " Orthogonality error: " << maxOrth;
539  std::cout << endl;
540  }*/
541  if ((maxNorm <= tol) && (maxOrth <= tol)) {
542  break;
543  }
544  } // for (rank = size; rank > 0; --rank)
545  //
546  // Copy the computed eigenvectors and eigenvalues
547  // ( they may be less than the number requested because of deflation )
548  //
549  // std::cout << "directSolve rank: " << rank << "\tsize: " << size << endl;
550  nev = (rank < nev) ? rank : nev;
551  EV.putScalar( zero );
552  std::copy(tt.begin(),tt.begin()+nev,theta.begin());
553  for (int i = 0; i < nev; ++i) {
554  blas.COPY( rank, (*KKcopy)[i], 1, EV[i], 1 );
555  }
556  break;
557 
558  case 1:
559  //
560  // Use the Cholesky factorization of MM to compute the generalized eigenvectors
561  //
562  // Copy KK & MM
563  //
564  KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
565  MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, size, size ) );
566  //
567  // Solve the generalized eigenproblem with LAPACK
568  //
569  info = 0;
570  lapack.HEGV(1, 'V', 'U', size, KKcopy->values(), KKcopy->stride(),
571  MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
572  &rwork[0], &info);
573  //
574  // Treat error messages
575  //
576  if (info < 0) {
577  std::cerr << std::endl;
578  std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
579  std::cerr << std::endl;
580  return -20;
581  }
582  if (info > 0) {
583  if (info > size)
584  nev = 0;
585  else {
586  std::cerr << std::endl;
587  std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, DPOTRF or DHEEV returned an error code (" << info << ").\n";
588  std::cerr << std::endl;
589  return -20;
590  }
591  }
592  //
593  // Copy the eigenvectors and eigenvalues
594  //
595  std::copy(tt.begin(),tt.begin()+nev,theta.begin());
596  for (int i = 0; i < nev; ++i) {
597  blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
598  }
599  break;
600 
601  case 10:
602  //
603  // Simple eigenproblem
604  //
605  // Copy KK
606  //
607  KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
608  //
609  // Solve the generalized eigenproblem with LAPACK
610  //
611  lapack.HEEV('V', 'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
612  //
613  // Treat error messages
614  if (info != 0) {
615  std::cerr << std::endl;
616  if (info < 0) {
617  std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, argument " << -info << " has an illegal value\n";
618  }
619  else {
620  std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, the algorithm failed to converge (" << info << ").\n";
621  }
622  std::cerr << std::endl;
623  info = -20;
624  break;
625  }
626  //
627  // Copy the eigenvectors
628  //
629  std::copy(tt.begin(),tt.begin()+nev,theta.begin());
630  for (int i = 0; i < nev; ++i) {
631  blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
632  }
633  break;
634  }
635 
636  return info;
637  }
638 
639 
640  //-----------------------------------------------------------------------------
641  //
642  // SANITY CHECKING METHODS
643  //
644  //-----------------------------------------------------------------------------
645 
646  template<class ScalarType, class MV, class OP>
647  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
648  SolverUtils<ScalarType, MV, OP>::errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M)
649  {
650  // Return the maximum coefficient of the matrix M * X - MX
651  // scaled by the maximum coefficient of MX.
652  // When M is not specified, the identity is used.
653 
654  MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
655 
656  int xc = MVT::GetNumberVecs(X);
657  int mxc = MVT::GetNumberVecs(MX);
658 
659  TEUCHOS_TEST_FOR_EXCEPTION(xc != mxc,std::invalid_argument,"Anasazi::SolverUtils::errorEquality(): input multivecs have different number of columns.");
660  if (xc == 0) {
661  return maxDiff;
662  }
663 
664  MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
665  std::vector<MagnitudeType> tmp( xc );
666  MVT::MvNorm(MX, tmp);
667 
668  for (int i = 0; i < xc; ++i) {
669  maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
670  }
671 
672  std::vector<int> index( 1 );
673  Teuchos::RCP<MV> MtimesX;
674  if (M != Teuchos::null) {
675  MtimesX = MVT::Clone( X, xc );
676  OPT::Apply( *M, X, *MtimesX );
677  }
678  else {
679  MtimesX = MVT::CloneCopy(X);
680  }
681  MVT::MvAddMv( -1.0, MX, 1.0, *MtimesX, *MtimesX );
682  MVT::MvNorm( *MtimesX, tmp );
683 
684  for (int i = 0; i < xc; ++i) {
685  maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
686  }
687 
688  return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
689 
690  }
691 
692 } // end namespace Anasazi
693 
694 #endif // ANASAZI_SOLVER_UTILS_HPP
695 
Declaration of basic traits for the multivector type.
virtual ~SolverUtils()
Destructor.
Virtual base class which defines basic traits for the operator type.
static Teuchos::ScalarTraits< ScalarType >::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP< const OP > M=Teuchos::null)
Return the maximum coefficient of the matrix scaled by the maximum coefficient of MX...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi&#39;s templated, static class providing utilities for the solvers.
Abstract class definition for Anasazi Output Managers.
static void applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix< int, ScalarType > &H, const std::vector< ScalarType > &tau, Teuchos::RCP< MV > workMV=Teuchos::null)
Apply a sequence of Householder reflectors (from GEQRF) to a multivector, using minimal workspace...
static int directSolver(int size, const Teuchos::SerialDenseMatrix< int, ScalarType > &KK, Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > MM, Teuchos::SerialDenseMatrix< int, ScalarType > &EV, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &theta, int &nev, int esType=0)
Routine for computing the first NEV generalized eigenpairs of the Hermitian pencil (KK...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
static void permuteVectors(const int n, const std::vector< int > &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > *resids=0)
Permute the vectors in a multivector according to the permutation vector perm, and optionally the res...