29 #ifndef ANASAZI_SOLVER_UTILS_HPP 30 #define ANASAZI_SOLVER_UTILS_HPP 51 #include "Teuchos_ScalarTraits.hpp" 54 #include "Teuchos_BLAS.hpp" 55 #include "Teuchos_LAPACK.hpp" 56 #include "Teuchos_SerialDenseMatrix.hpp" 60 template<
class ScalarType,
class MV,
class OP>
64 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
65 typedef typename Teuchos::ScalarTraits<ScalarType> SCT;
82 static void permuteVectors(
const int n,
const std::vector<int> &perm, MV &Q, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids = 0);
85 static void permuteVectors(
const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q);
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);
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);
161 static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
errorEquality(
const MV &X,
const MV &MX, Teuchos::RCP<const OP> M = Teuchos::null);
182 template<
class ScalarType,
class MV,
class OP>
194 template<
class ScalarType,
class MV,
class OP>
197 const std::vector<int> &perm,
199 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids)
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();
210 TEUCHOS_TEST_FOR_EXCEPTION(n >
MVT::GetNumberVecs(Q), std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): argument n larger than width of input multivector.");
216 for (i=0; i<n-1; i++) {
219 for (j=i; j<n; j++) {
220 if (permcopy[j] == i) {
224 TEUCHOS_TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): permutation index invalid.");
228 std::swap( permcopy[j], permcopy[i] );
234 for (i=n-2; i>=0; i--) {
241 std::swap( (*resids)[i], (*resids)[j] );
258 template<
class ScalarType,
class MV,
class OP>
260 const std::vector<int> &perm,
261 Teuchos::SerialDenseMatrix<int,ScalarType> &Q)
265 Teuchos::BLAS<int,ScalarType> blas;
266 const int n = perm.size();
267 const int m = Q.numRows();
269 TEUCHOS_TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
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);
286 template<
class ScalarType,
class MV,
class OP>
290 const ScalarType ONE = SCT::one();
291 const ScalarType ZERO = SCT::zero();
298 if (workMV == Teuchos::null) {
303 std::vector<int> first(1);
308 TEUCHOS_TEST_FOR_EXCEPTION(
MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,
"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
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.");
318 for (
int i=0; i<k; i++) {
321 std::vector<int> activeind(n-i);
322 for (
int j=0; j<n-i; j++) activeind[j] = j+i;
329 Teuchos::SerialDenseMatrix<int,ScalarType> v(Teuchos::Copy,H,n-i,1,i,i);
341 Teuchos::SerialDenseMatrix<int,ScalarType> vT(v,Teuchos::CONJ_TRANS);
344 actV = Teuchos::null;
355 template<
class ScalarType,
class MV,
class OP>
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)
406 Teuchos::LAPACK<int,ScalarType> lapack;
407 Teuchos::BLAS<int,ScalarType> blas;
412 if (size < nev || size < 0) {
415 if (KK.numCols() < size || KK.numRows() < size) {
418 if ((esType == 0 || esType == 1)) {
419 if (MM == Teuchos::null) {
422 else if (MM->numCols() < size || MM->numRows() < size) {
426 if (EV.numCols() < size || EV.numRows() < size) {
429 if (theta.size() < (
unsigned int) size) {
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);
441 std::vector<ScalarType> work(lwork);
442 std::vector<MagnitudeType> rwork(3*size-2);
445 std::vector<MagnitudeType> tt( size );
448 MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
450 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
451 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
453 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
454 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
462 for (rank = size; rank > 0; --rank) {
464 U = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );
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 ) );
474 lapack.HEGV(1,
'V',
'U', rank, KKcopy->values(), KKcopy->stride(),
475 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
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;
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));
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.");
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) {
513 maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
514 ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;
516 maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
517 ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
528 if ((maxNorm <= tol) && (maxOrth <= tol)) {
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 );
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 ) );
557 lapack.HEGV(1,
'V',
'U', size, KKcopy->values(), KKcopy->stride(),
558 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
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;
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;
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 );
594 KKcopy = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
598 lapack.HEEV(
'V',
'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
602 std::cerr << std::endl;
604 std::cerr <<
"Anasazi::SolverUtils::directSolver(): In DHEEV, argument " << -info <<
" has an illegal value\n";
607 std::cerr <<
"Anasazi::SolverUtils::directSolver(): In DHEEV, the algorithm failed to converge (" << info <<
").\n";
609 std::cerr << std::endl;
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 );
633 template<
class ScalarType,
class MV,
class OP>
634 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
641 MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
646 TEUCHOS_TEST_FOR_EXCEPTION(xc != mxc,std::invalid_argument,
"Anasazi::SolverUtils::errorEquality(): input multivecs have different number of columns.");
651 MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
652 std::vector<MagnitudeType> tmp( xc );
655 for (
int i = 0; i < xc; ++i) {
656 maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
659 std::vector<int> index( 1 );
660 Teuchos::RCP<MV> MtimesX;
661 if (M != Teuchos::null) {
671 for (
int i = 0; i < xc; ++i) {
672 maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
675 return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
681 #endif // ANASAZI_SOLVER_UTILS_HPP 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'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 .
SolverUtils()
Constructor.
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.