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 
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 
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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...
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
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.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
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 Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
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...
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.