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, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef ANASAZI_SOLVER_UTILS_HPP
30 #define ANASAZI_SOLVER_UTILS_HPP
31 
48 #include "AnasaziConfigDefs.hpp"
51 #include "Teuchos_ScalarTraits.hpp"
52 
53 #include "AnasaziOutputManager.hpp"
54 #include "Teuchos_BLAS.hpp"
55 #include "Teuchos_LAPACK.hpp"
56 #include "Teuchos_SerialDenseMatrix.hpp"
57 
58 namespace Anasazi {
59 
60  template<class ScalarType, class MV, class OP>
62  {
63  public:
64  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
65  typedef typename Teuchos::ScalarTraits<ScalarType> SCT;
66 
68 
69 
71  SolverUtils();
72 
74  virtual ~SolverUtils() {};
75 
77 
79 
80 
82  static void permuteVectors(const int n, const std::vector<int> &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids = 0);
83 
85  static void permuteVectors(const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q);
86 
88 
90 
91 
93 
114  static void applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix<int,ScalarType> &H, const std::vector<ScalarType> &tau, Teuchos::RCP<MV> workMV = Teuchos::null);
115 
117 
119 
120 
122 
148  static int directSolver(int size, const Teuchos::SerialDenseMatrix<int,ScalarType> &KK,
149  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
150  Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
151  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
152  int &nev, int esType = 0);
154 
156 
157 
159 
161  static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M = Teuchos::null);
162 
164 
165  private:
166 
168 
169 
170  typedef MultiVecTraits<ScalarType,MV> MVT;
172 
174  };
175 
176  //-----------------------------------------------------------------------------
177  //
178  // CONSTRUCTOR
179  //
180  //-----------------------------------------------------------------------------
181 
182  template<class ScalarType, class MV, class OP>
184 
185 
186  //-----------------------------------------------------------------------------
187  //
188  // SORTING METHODS
189  //
190  //-----------------------------------------------------------------------------
191 
193  // permuteVectors for MV
194  template<class ScalarType, class MV, class OP>
196  const int n,
197  const std::vector<int> &perm,
198  MV &Q,
199  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids)
200  {
201  // Permute the vectors according to the permutation vector \c perm, and
202  // optionally the residual vector \c resids
203 
204  int i, j;
205  std::vector<int> permcopy(perm), swapvec(n-1);
206  std::vector<int> index(1);
207  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
208  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
209 
210  TEUCHOS_TEST_FOR_EXCEPTION(n > MVT::GetNumberVecs(Q), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): argument n larger than width of input multivector.");
211 
212  // We want to recover the elementary permutations (individual swaps)
213  // from the permutation vector. Do this by constructing the inverse
214  // of the permutation, by sorting them to {1,2,...,n}, and recording
215  // the elementary permutations of the inverse.
216  for (i=0; i<n-1; i++) {
217  //
218  // find i in the permcopy vector
219  for (j=i; j<n; j++) {
220  if (permcopy[j] == i) {
221  // found it at index j
222  break;
223  }
224  TEUCHOS_TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): permutation index invalid.");
225  }
226  //
227  // Swap two scalars
228  std::swap( permcopy[j], permcopy[i] );
229 
230  swapvec[i] = j;
231  }
232 
233  // now apply the elementary permutations of the inverse in reverse order
234  for (i=n-2; i>=0; i--) {
235  j = swapvec[i];
236  //
237  // Swap (i,j)
238  //
239  // Swap residuals (if they exist)
240  if (resids) {
241  std::swap( (*resids)[i], (*resids)[j] );
242  }
243  //
244  // Swap corresponding vectors
245  index[0] = j;
246  Teuchos::RCP<MV> tmpQ = MVT::CloneCopy( Q, index );
247  Teuchos::RCP<MV> tmpQj = MVT::CloneViewNonConst( Q, index );
248  index[0] = i;
249  Teuchos::RCP<MV> tmpQi = MVT::CloneViewNonConst( Q, index );
250  MVT::MvAddMv( one, *tmpQi, zero, *tmpQi, *tmpQj );
251  MVT::MvAddMv( one, *tmpQ, zero, *tmpQ, *tmpQi );
252  }
253  }
254 
255 
257  // permuteVectors for MV
258  template<class ScalarType, class MV, class OP>
260  const std::vector<int> &perm,
261  Teuchos::SerialDenseMatrix<int,ScalarType> &Q)
262  {
263  // Permute the vectors in Q according to the permutation vector \c perm, and
264  // optionally the residual vector \c resids
265  Teuchos::BLAS<int,ScalarType> blas;
266  const int n = perm.size();
267  const int m = Q.numRows();
268 
269  TEUCHOS_TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
270 
271  // Sort the primitive ritz vectors
272  Teuchos::SerialDenseMatrix<int,ScalarType> copyQ(Teuchos::Copy, Q);
273  for (int i=0; i<n; i++) {
274  blas.COPY(m, copyQ[perm[i]], 1, Q[i], 1);
275  }
276  }
277 
278 
279  //-----------------------------------------------------------------------------
280  //
281  // BASIS UPDATE METHODS
282  //
283  //-----------------------------------------------------------------------------
284 
285  // apply householder reflectors to multivector
286  template<class ScalarType, class MV, class OP>
287  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) {
288 
289  const int n = MVT::GetNumberVecs(V);
290  const ScalarType ONE = SCT::one();
291  const ScalarType ZERO = SCT::zero();
292 
293  // early exit if V has zero-size or if k==0
294  if (MVT::GetNumberVecs(V) == 0 || MVT::GetGlobalLength(V) == 0 || k == 0) {
295  return;
296  }
297 
298  if (workMV == Teuchos::null) {
299  // user did not give us any workspace; allocate some
300  workMV = MVT::Clone(V,1);
301  }
302  else if (MVT::GetNumberVecs(*workMV) > 1) {
303  std::vector<int> first(1);
304  first[0] = 0;
305  workMV = MVT::CloneViewNonConst(*workMV,first);
306  }
307  else {
308  TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
309  }
310  // Q = H_1 ... H_k is square, with as many rows as V has vectors
311  // however, H need only have k columns, one each for the k reflectors.
312  TEUCHOS_TEST_FOR_EXCEPTION( H.numCols() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): H must have at least k columns.");
313  TEUCHOS_TEST_FOR_EXCEPTION( (int)tau.size() != k, std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): tau must have at least k entries.");
314  TEUCHOS_TEST_FOR_EXCEPTION( H.numRows() != MVT::GetNumberVecs(V), std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): Size of H,V are inconsistent.");
315 
316  // perform the loop
317  // flops: Sum_{i=0:k-1} 4 m (n-i) == 4mnk - 2m(k^2- k)
318  for (int i=0; i<k; i++) {
319  // apply V H_i+1 = V - tau_i+1 (V v_i+1) v_i+1^T
320  // because of the structure of v_i+1, this transform does not affect the first i columns of V
321  std::vector<int> activeind(n-i);
322  for (int j=0; j<n-i; j++) activeind[j] = j+i;
323  Teuchos::RCP<MV> actV = MVT::CloneViewNonConst(V,activeind);
324 
325  // note, below H_i, v_i and tau_i are mathematical objects which use 1-based indexing
326  // while H, v and tau are data structures using 0-based indexing
327 
328  // get v_i+1: i-th column of H
329  Teuchos::SerialDenseMatrix<int,ScalarType> v(Teuchos::Copy,H,n-i,1,i,i);
330  // v_i+1(1:i) = 0: this isn't part of v
331  // e_i+1^T v_i+1 = 1 = v(0)
332  v(0,0) = ONE;
333 
334  // compute -tau_i V v_i
335  // tau_i+1 is tau[i]
336  // flops: 2 m n-i
337  MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
338 
339  // perform V = V + workMV v_i^T
340  // flops: 2 m n-i
341  Teuchos::SerialDenseMatrix<int,ScalarType> vT(v,Teuchos::CONJ_TRANS);
342  MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
343 
344  actV = Teuchos::null;
345  }
346  }
347 
348 
349  //-----------------------------------------------------------------------------
350  //
351  // EIGENSOLVER PROJECTION METHODS
352  //
353  //-----------------------------------------------------------------------------
354 
355  template<class ScalarType, class MV, class OP>
357  int size,
358  const Teuchos::SerialDenseMatrix<int,ScalarType> &KK,
359  Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > MM,
360  Teuchos::SerialDenseMatrix<int,ScalarType> &EV,
361  std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType > &theta,
362  int &nev, int esType)
363  {
364  // Routine for computing the first NEV generalized eigenpairs of the symmetric pencil (KK, MM)
365  //
366  // Parameter variables:
367  //
368  // size : Dimension of the eigenproblem (KK, MM)
369  //
370  // KK : Hermitian "stiffness" matrix
371  //
372  // MM : Hermitian positive-definite "mass" matrix
373  //
374  // EV : Matrix to store the nev eigenvectors
375  //
376  // theta : Array to store the eigenvalues (Size = nev )
377  //
378  // nev : Number of the smallest eigenvalues requested (input)
379  // Number of the smallest computed eigenvalues (output)
380  // Routine may compute and return more or less eigenvalues than requested.
381  //
382  // esType : Flag to select the algorithm
383  //
384  // esType = 0 (default) Uses LAPACK routine (Cholesky factorization of MM)
385  // with deflation of MM to get orthonormality of
386  // eigenvectors (S^T MM S = I)
387  //
388  // esType = 1 Uses LAPACK routine (Cholesky factorization of MM)
389  // (no check of orthonormality)
390  //
391  // esType = 10 Uses LAPACK routine for simple eigenproblem on KK
392  // (MM is not referenced in this case)
393  //
394  // Note: The code accesses only the upper triangular part of KK and MM.
395  //
396  // Return the integer info on the status of the computation
397  //
398  // info = 0 >> Success
399  //
400  // info < 0 >> error in the info-th argument
401  // info = - 20 >> Failure in LAPACK routine
402 
403  // Define local arrays
404 
405  // Create blas/lapack objects.
406  Teuchos::LAPACK<int,ScalarType> lapack;
407  Teuchos::BLAS<int,ScalarType> blas;
408 
409  int rank = 0;
410  int info = 0;
411 
412  if (size < nev || size < 0) {
413  return -1;
414  }
415  if (KK.numCols() < size || KK.numRows() < size) {
416  return -2;
417  }
418  if ((esType == 0 || esType == 1)) {
419  if (MM == Teuchos::null) {
420  return -3;
421  }
422  else if (MM->numCols() < size || MM->numRows() < size) {
423  return -3;
424  }
425  }
426  if (EV.numCols() < size || EV.numRows() < size) {
427  return -4;
428  }
429  if (theta.size() < (unsigned int) size) {
430  return -5;
431  }
432  if (nev <= 0) {
433  return -6;
434  }
435 
436  // Query LAPACK for the "optimal" block size for HEGV
437  std::string lapack_name = "hetrd";
438  std::string lapack_opts = "u";
439  int NB = lapack.ILAENV(1, lapack_name, lapack_opts, size, -1, -1, -1);
440  int lwork = size*(NB+2); // For HEEV, lwork should be NB+2, instead of NB+1
441  std::vector<ScalarType> work(lwork);
442  std::vector<MagnitudeType> rwork(3*size-2);
443  // tt contains the eigenvalues from HEGV, which are necessarily real, and
444  // HEGV expects this vector to be real as well
445  std::vector<MagnitudeType> tt( size );
446  //typedef typename std::vector<MagnitudeType>::iterator MTIter; // unused
447 
448  MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
449  // MagnitudeType tol = 1e-12;
450  ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
451  ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
452 
453  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
454  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
455 
456  switch (esType) {
457  default:
458  case 0:
459  //
460  // Use LAPACK to compute the generalized eigenvectors
461  //
462  for (rank = size; rank > 0; --rank) {
463 
464  U = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );
465  //
466  // Copy KK & MM
467  //
468  KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, rank, rank ) );
469  MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
470  //
471  // Solve the generalized eigenproblem with LAPACK
472  //
473  info = 0;
474  lapack.HEGV(1, 'V', 'U', rank, KKcopy->values(), KKcopy->stride(),
475  MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
476  &rwork[0], &info);
477  //
478  // Treat error messages
479  //
480  if (info < 0) {
481  std::cerr << std::endl;
482  std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
483  std::cerr << std::endl;
484  return -20;
485  }
486  if (info > 0) {
487  if (info > rank)
488  rank = info - rank;
489  continue;
490  }
491  //
492  // Check the quality of eigenvectors ( using mass-orthonormality )
493  //
494  MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, rank, rank ) );
495  for (int i = 0; i < rank; ++i) {
496  for (int j = 0; j < i; ++j) {
497  (*MMcopy)(i,j) = SCT::conjugate((*MM)(j,i));
498  }
499  }
500  // U = 0*U + 1*MMcopy*KKcopy = MMcopy * KKcopy
501  TEUCHOS_TEST_FOR_EXCEPTION(
502  U->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*MMcopy,*KKcopy,zero) != 0,
503  std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
504  // MMcopy = 0*MMcopy + 1*KKcopy^H*U = KKcopy^H * MMcopy * KKcopy
505  TEUCHOS_TEST_FOR_EXCEPTION(
506  MMcopy->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,one,*KKcopy,*U,zero) != 0,
507  std::logic_error, "Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
508  MagnitudeType maxNorm = SCT::magnitude(zero);
509  MagnitudeType maxOrth = SCT::magnitude(zero);
510  for (int i = 0; i < rank; ++i) {
511  for (int j = i; j < rank; ++j) {
512  if (j == i)
513  maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
514  ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;
515  else
516  maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
517  ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
518  }
519  }
520  /* if (verbose > 4) {
521  std::cout << " >> Local eigensolve >> Size: " << rank;
522  std::cout.precision(2);
523  std::cout.setf(std::ios::scientific, std::ios::floatfield);
524  std::cout << " Normalization error: " << maxNorm;
525  std::cout << " Orthogonality error: " << maxOrth;
526  std::cout << endl;
527  }*/
528  if ((maxNorm <= tol) && (maxOrth <= tol)) {
529  break;
530  }
531  } // for (rank = size; rank > 0; --rank)
532  //
533  // Copy the computed eigenvectors and eigenvalues
534  // ( they may be less than the number requested because of deflation )
535  //
536  // std::cout << "directSolve rank: " << rank << "\tsize: " << size << endl;
537  nev = (rank < nev) ? rank : nev;
538  EV.putScalar( zero );
539  std::copy(tt.begin(),tt.begin()+nev,theta.begin());
540  for (int i = 0; i < nev; ++i) {
541  blas.COPY( rank, (*KKcopy)[i], 1, EV[i], 1 );
542  }
543  break;
544 
545  case 1:
546  //
547  // Use the Cholesky factorization of MM to compute the generalized eigenvectors
548  //
549  // Copy KK & MM
550  //
551  KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
552  MMcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *MM, size, size ) );
553  //
554  // Solve the generalized eigenproblem with LAPACK
555  //
556  info = 0;
557  lapack.HEGV(1, 'V', 'U', size, KKcopy->values(), KKcopy->stride(),
558  MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
559  &rwork[0], &info);
560  //
561  // Treat error messages
562  //
563  if (info < 0) {
564  std::cerr << std::endl;
565  std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, argument " << -info << "has an illegal value.\n";
566  std::cerr << std::endl;
567  return -20;
568  }
569  if (info > 0) {
570  if (info > size)
571  nev = 0;
572  else {
573  std::cerr << std::endl;
574  std::cerr << "Anasazi::SolverUtils::directSolver(): In HEGV, DPOTRF or DHEEV returned an error code (" << info << ").\n";
575  std::cerr << std::endl;
576  return -20;
577  }
578  }
579  //
580  // Copy the eigenvectors and eigenvalues
581  //
582  std::copy(tt.begin(),tt.begin()+nev,theta.begin());
583  for (int i = 0; i < nev; ++i) {
584  blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
585  }
586  break;
587 
588  case 10:
589  //
590  // Simple eigenproblem
591  //
592  // Copy KK
593  //
594  KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
595  //
596  // Solve the generalized eigenproblem with LAPACK
597  //
598  lapack.HEEV('V', 'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
599  //
600  // Treat error messages
601  if (info != 0) {
602  std::cerr << std::endl;
603  if (info < 0) {
604  std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, argument " << -info << " has an illegal value\n";
605  }
606  else {
607  std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, the algorithm failed to converge (" << info << ").\n";
608  }
609  std::cerr << std::endl;
610  info = -20;
611  break;
612  }
613  //
614  // Copy the eigenvectors
615  //
616  std::copy(tt.begin(),tt.begin()+nev,theta.begin());
617  for (int i = 0; i < nev; ++i) {
618  blas.COPY( size, (*KKcopy)[i], 1, EV[i], 1 );
619  }
620  break;
621  }
622 
623  return info;
624  }
625 
626 
627  //-----------------------------------------------------------------------------
628  //
629  // SANITY CHECKING METHODS
630  //
631  //-----------------------------------------------------------------------------
632 
633  template<class ScalarType, class MV, class OP>
634  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
635  SolverUtils<ScalarType, MV, OP>::errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M)
636  {
637  // Return the maximum coefficient of the matrix M * X - MX
638  // scaled by the maximum coefficient of MX.
639  // When M is not specified, the identity is used.
640 
641  MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
642 
643  int xc = MVT::GetNumberVecs(X);
644  int mxc = MVT::GetNumberVecs(MX);
645 
646  TEUCHOS_TEST_FOR_EXCEPTION(xc != mxc,std::invalid_argument,"Anasazi::SolverUtils::errorEquality(): input multivecs have different number of columns.");
647  if (xc == 0) {
648  return maxDiff;
649  }
650 
651  MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
652  std::vector<MagnitudeType> tmp( xc );
653  MVT::MvNorm(MX, tmp);
654 
655  for (int i = 0; i < xc; ++i) {
656  maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
657  }
658 
659  std::vector<int> index( 1 );
660  Teuchos::RCP<MV> MtimesX;
661  if (M != Teuchos::null) {
662  MtimesX = MVT::Clone( X, xc );
663  OPT::Apply( *M, X, *MtimesX );
664  }
665  else {
666  MtimesX = MVT::CloneCopy(X);
667  }
668  MVT::MvAddMv( -1.0, MX, 1.0, *MtimesX, *MtimesX );
669  MVT::MvNorm( *MtimesX, tmp );
670 
671  for (int i = 0; i < xc; ++i) {
672  maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
673  }
674 
675  return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
676 
677  }
678 
679 } // end namespace Anasazi
680 
681 #endif // ANASAZI_SOLVER_UTILS_HPP
682 
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...