29 #ifndef ANASAZI_SOLVER_UTILS_HPP 30 #define ANASAZI_SOLVER_UTILS_HPP 60 template<
class ScalarType,
class MV,
class OP>
152 int &nev,
int esType = 0);
182 template<
class ScalarType,
class MV,
class OP>
194 template<
class ScalarType,
class MV,
class OP>
197 const std::vector<int> &perm,
205 std::vector<int> permcopy(perm), swapvec(n-1);
206 std::vector<int> index(1);
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] );
250 MVT::MvAddMv( one, *tmpQi, zero, *tmpQi, *tmpQj );
251 MVT::MvAddMv( one, *tmpQ, zero, *tmpQ, *tmpQi );
258 template<
class ScalarType,
class MV,
class OP>
260 const std::vector<int> &perm,
266 const int n = perm.size();
269 TEUCHOS_TEST_FOR_EXCEPTION(n != Q.
numCols(), std::invalid_argument,
"Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
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>
289 const int n = MVT::GetNumberVecs(V);
290 const ScalarType ONE = SCT::one();
291 const ScalarType ZERO = SCT::zero();
294 if (MVT::GetNumberVecs(V) == 0 || MVT::GetGlobalLength(V) == 0 || k == 0) {
298 if (workMV == Teuchos::null) {
300 workMV = MVT::Clone(V,1);
302 else if (MVT::GetNumberVecs(*workMV) > 1) {
303 std::vector<int> first(1);
305 workMV = MVT::CloneViewNonConst(*workMV,first);
308 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,
"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
313 TEUCHOS_TEST_FOR_EXCEPTION( (
int)tau.size() != k, std::invalid_argument,
"Anasazi::SolverUtils::applyHouse(): tau must have at least k entries.");
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;
337 MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
342 MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
344 actV = Teuchos::null;
355 template<
class ScalarType,
class MV,
class OP>
362 int &nev,
int esType)
412 if (size < nev || size < 0) {
418 if ((esType == 0 || esType == 1)) {
419 if (MM == Teuchos::null) {
422 else if (MM->numCols() < size || MM->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()));
462 for (rank = size; rank > 0; --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;
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));
503 std::logic_error,
"Anasazi::SolverUtils::directSolver() call to Teuchos::SerialDenseMatrix::multiply() returned an error.");
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;
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 );
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 );
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>
641 MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
643 int xc = MVT::GetNumberVecs(X);
644 int mxc = MVT::GetNumberVecs(MX);
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 );
653 MVT::MvNorm(MX, tmp);
655 for (
int i = 0; i < xc; ++i) {
656 maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
659 std::vector<int> index( 1 );
661 if (M != Teuchos::null) {
662 MtimesX = MVT::Clone( X, xc );
663 OPT::Apply( *M, X, *MtimesX );
666 MtimesX = MVT::CloneCopy(X);
668 MVT::MvAddMv( -1.0, MX, 1.0, *MtimesX, *MtimesX );
669 MVT::MvNorm( *MtimesX, tmp );
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
Declaration of basic traits for the multivector type.
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
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...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi's templated, static class providing utilities for the solvers.
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
Abstract class definition for Anasazi Output Managers.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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...
void COPY(const OrdinalType n, const ScalarType *x, const OrdinalType incx, ScalarType *y, const OrdinalType incy) const
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...
OrdinalType numCols() 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