42 #ifndef ANASAZI_SOLVER_UTILS_HPP 43 #define ANASAZI_SOLVER_UTILS_HPP 73 template<
class ScalarType,
class MV,
class OP>
165 int &nev,
int esType = 0);
195 template<
class ScalarType,
class MV,
class OP>
207 template<
class ScalarType,
class MV,
class OP>
210 const std::vector<int> &perm,
218 std::vector<int> permcopy(perm), swapvec(n-1);
219 std::vector<int> index(1);
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] );
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,
279 const int n = perm.size();
282 TEUCHOS_TEST_FOR_EXCEPTION(n != Q.
numCols(), std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
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.");
326 TEUCHOS_TEST_FOR_EXCEPTION( (
int)tau.size() != k, std::invalid_argument,
"Anasazi::SolverUtils::applyHouse(): tau must have at least k entries.");
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;
350 MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
355 MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
357 actV = Teuchos::null;
368 template<
class ScalarType,
class MV,
class OP>
375 int &nev,
int esType)
425 if (size < nev || size < 0) {
431 if ((esType == 0 || esType == 1)) {
432 if (MM == Teuchos::null) {
435 else if (MM->numCols() < size || MM->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()));
475 for (rank = size; rank > 0; --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;
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));
516 std::logic_error,
"Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
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;
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 );
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 );
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>
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 );
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
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
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...
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
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
OrdinalType numRows() const
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...
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...
void COPY(const OrdinalType &n, const ScalarType *x, const OrdinalType &incx, ScalarType *y, const OrdinalType &incy) const
OrdinalType numCols() const