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);
183 typedef MultiVecTraits<ScalarType,MV> MVT;
184 typedef OperatorTraits<ScalarType,MV,OP> OPT;
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 );
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;
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 );
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;