42 #ifndef ANASAZI_SOLVER_UTILS_HPP 43 #define ANASAZI_SOLVER_UTILS_HPP 64 #include "Teuchos_ScalarTraits.hpp" 67 #include "Teuchos_BLAS.hpp" 68 #include "Teuchos_LAPACK.hpp" 69 #include "Teuchos_SerialDenseMatrix.hpp" 73 template<
class ScalarType,
class MV,
class OP>
77 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
78 typedef typename Teuchos::ScalarTraits<ScalarType> SCT;
95 static void permuteVectors(
const int n,
const std::vector<int> &perm, MV &Q, std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids = 0);
98 static void permuteVectors(
const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q);
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);
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);
174 static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
errorEquality(
const MV &X,
const MV &MX, Teuchos::RCP<const OP> M = Teuchos::null);
195 template<
class ScalarType,
class MV,
class OP>
207 template<
class ScalarType,
class MV,
class OP>
210 const std::vector<int> &perm,
212 std::vector<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids)
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();
223 TEUCHOS_TEST_FOR_EXCEPTION(n > MVT::GetNumberVecs(Q), std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): argument n larger than width of input multivector.");
229 for (i=0; i<n-1; i++) {
232 for (j=i; j<n; j++) {
233 if (permcopy[j] == i) {
237 TEUCHOS_TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): permutation index invalid.");
241 std::swap( permcopy[j], permcopy[i] );
247 for (i=n-2; i>=0; i--) {
254 std::swap( (*resids)[i], (*resids)[j] );
259 Teuchos::RCP<MV> tmpQ = MVT::CloneCopy( Q, index );
260 Teuchos::RCP<MV> tmpQj = MVT::CloneViewNonConst( Q, index );
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 );
271 template<
class ScalarType,
class MV,
class OP>
273 const std::vector<int> &perm,
274 Teuchos::SerialDenseMatrix<int,ScalarType> &Q)
278 Teuchos::BLAS<int,ScalarType> blas;
279 const int n = perm.size();
280 const int m = Q.numRows();
282 TEUCHOS_TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
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);
299 template<
class ScalarType,
class MV,
class OP>
302 const int n = MVT::GetNumberVecs(V);
303 const ScalarType ONE = SCT::one();
304 const ScalarType ZERO = SCT::zero();
307 if (MVT::GetNumberVecs(V) == 0 || MVT::GetGlobalLength(V) == 0 || k == 0) {
311 if (workMV == Teuchos::null) {
313 workMV = MVT::Clone(V,1);
315 else if (MVT::GetNumberVecs(*workMV) > 1) {
316 std::vector<int> first(1);
318 workMV = MVT::CloneViewNonConst(*workMV,first);
321 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,
"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
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.");
331 for (
int i=0; i<k; i++) {
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);
342 Teuchos::SerialDenseMatrix<int,ScalarType> v(Teuchos::Copy,H,n-i,1,i,i);
350 MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
354 Teuchos::SerialDenseMatrix<int,ScalarType> vT(v,Teuchos::CONJ_TRANS);
355 MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
357 actV = Teuchos::null;
368 template<
class ScalarType,
class MV,
class OP>
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)
419 Teuchos::LAPACK<int,ScalarType> lapack;
420 Teuchos::BLAS<int,ScalarType> blas;
425 if (size < nev || size < 0) {
428 if (KK.numCols() < size || KK.numRows() < size) {
431 if ((esType == 0 || esType == 1)) {
432 if (MM == Teuchos::null) {
435 else if (MM->numCols() < size || MM->numRows() < size) {
439 if (EV.numCols() < size || EV.numRows() < size) {
442 if (theta.size() < (
unsigned int) size) {
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);
454 std::vector<ScalarType> work(lwork);
455 std::vector<MagnitudeType> rwork(3*size-2);
458 std::vector<MagnitudeType> tt( size );
461 MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
463 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
464 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
466 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
467 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
475 for (rank = size; rank > 0; --rank) {
477 U = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );
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 ) );
487 lapack.HEGV(1,
'V',
'U', rank, KKcopy->values(), KKcopy->stride(),
488 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
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;
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));
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.");
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) {
526 maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
527 ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;
529 maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
530 ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
541 if ((maxNorm <= tol) && (maxOrth <= tol)) {
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 );
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 ) );
570 lapack.HEGV(1,
'V',
'U', size, KKcopy->values(), KKcopy->stride(),
571 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
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;
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;
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 );
607 KKcopy = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
611 lapack.HEEV(
'V',
'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
615 std::cerr << std::endl;
617 std::cerr <<
"Anasazi::SolverUtils::directSolver(): In DHEEV, argument " << -info <<
" has an illegal value\n";
620 std::cerr <<
"Anasazi::SolverUtils::directSolver(): In DHEEV, the algorithm failed to converge (" << info <<
").\n";
622 std::cerr << std::endl;
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 );
646 template<
class ScalarType,
class MV,
class OP>
647 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
654 MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
656 int xc = MVT::GetNumberVecs(X);
657 int mxc = MVT::GetNumberVecs(MX);
659 TEUCHOS_TEST_FOR_EXCEPTION(xc != mxc,std::invalid_argument,
"Anasazi::SolverUtils::errorEquality(): input multivecs have different number of columns.");
664 MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
665 std::vector<MagnitudeType> tmp( xc );
666 MVT::MvNorm(MX, tmp);
668 for (
int i = 0; i < xc; ++i) {
669 maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
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 );
679 MtimesX = MVT::CloneCopy(X);
681 MVT::MvAddMv( -1.0, MX, 1.0, *MtimesX, *MtimesX );
682 MVT::MvNorm( *MtimesX, tmp );
684 for (
int i = 0; i < xc; ++i) {
685 maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
688 return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
694 #endif // ANASAZI_SOLVER_UTILS_HPP 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'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...
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...