Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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"
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,
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,
275  {
276  // Permute the vectors in Q according to the permutation vector \c perm, and
277  // optionally the residual vector \c resids
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
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
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
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,
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.
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 
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 
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
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
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>
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 
OrdinalType ILAENV(const OrdinalType &ispec, const std::string &NAME, const std::string &OPTS, const OrdinalType &N1=-1, const OrdinalType &N2=-1, const OrdinalType &N3=-1, const OrdinalType &N4=-1) const
Declaration of basic traits for the multivector type.
virtual ~SolverUtils()
Destructor.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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...
Anasazi&#39;s templated, static class providing utilities for the solvers.
void COPY(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
Abstract class definition for Anasazi Output Managers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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...
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
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...
OrdinalType numCols() const
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...
void HEEV(const char &JOBZ, const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, MagnitudeType *W, ScalarType *WORK, const OrdinalType &lwork, MagnitudeType *RWORK, OrdinalType *info) const
void HEGV(const OrdinalType &itype, const char &JOBZ, const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, ScalarType *B, const OrdinalType &ldb, MagnitudeType *W, ScalarType *WORK, const OrdinalType &lwork, MagnitudeType *RWORK, OrdinalType *info) const
OrdinalType numRows() const