Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziSolverUtils.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef ANASAZI_SOLVER_UTILS_HPP
43#define ANASAZI_SOLVER_UTILS_HPP
44
61#include "AnasaziConfigDefs.hpp"
64#include "Teuchos_ScalarTraits.hpp"
65
67#include "Teuchos_BLAS.hpp"
68#include "Teuchos_LAPACK.hpp"
69#include "Teuchos_SerialDenseMatrix.hpp"
70
71namespace Anasazi {
72
73 template<class ScalarType, class MV, class OP>
75 {
76 public:
77 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
78 typedef typename Teuchos::ScalarTraits<ScalarType> SCT;
79
81
82
85
87 virtual ~SolverUtils() {};
88
90
92
93
95 static void permuteVectors(const int n, const std::vector<int> &perm, MV &Q, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids = 0);
96
98 static void permuteVectors(const std::vector<int> &perm, Teuchos::SerialDenseMatrix<int,ScalarType> &Q);
99
101
103
104
106
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);
128
130
132
133
135
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);
167
169
170
172
174 static typename Teuchos::ScalarTraits<ScalarType>::magnitudeType errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M = Teuchos::null);
175
177
178 private:
179
181
182
183 typedef MultiVecTraits<ScalarType,MV> MVT;
184 typedef OperatorTraits<ScalarType,MV,OP> OPT;
185
187 };
188
189 //-----------------------------------------------------------------------------
190 //
191 // CONSTRUCTOR
192 //
193 //-----------------------------------------------------------------------------
194
195 template<class ScalarType, class MV, class OP>
197
198
199 //-----------------------------------------------------------------------------
200 //
201 // SORTING METHODS
202 //
203 //-----------------------------------------------------------------------------
204
206 // permuteVectors for MV
207 template<class ScalarType, class MV, class OP>
209 const int n,
210 const std::vector<int> &perm,
211 MV &Q,
212 std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >* resids)
213 {
214 // Permute the vectors according to the permutation vector \c perm, and
215 // optionally the residual vector \c resids
216
217 int i, j;
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();
222
223 TEUCHOS_TEST_FOR_EXCEPTION(n > MVT::GetNumberVecs(Q), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): argument n larger than width of input multivector.");
224
225 // We want to recover the elementary permutations (individual swaps)
226 // from the permutation vector. Do this by constructing the inverse
227 // of the permutation, by sorting them to {1,2,...,n}, and recording
228 // the elementary permutations of the inverse.
229 for (i=0; i<n-1; i++) {
230 //
231 // find i in the permcopy vector
232 for (j=i; j<n; j++) {
233 if (permcopy[j] == i) {
234 // found it at index j
235 break;
236 }
237 TEUCHOS_TEST_FOR_EXCEPTION(j == n-1, std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): permutation index invalid.");
238 }
239 //
240 // Swap two scalars
241 std::swap( permcopy[j], permcopy[i] );
242
243 swapvec[i] = j;
244 }
245
246 // now apply the elementary permutations of the inverse in reverse order
247 for (i=n-2; i>=0; i--) {
248 j = swapvec[i];
249 //
250 // Swap (i,j)
251 //
252 // Swap residuals (if they exist)
253 if (resids) {
254 std::swap( (*resids)[i], (*resids)[j] );
255 }
256 //
257 // Swap corresponding vectors
258 index[0] = j;
259 Teuchos::RCP<MV> tmpQ = MVT::CloneCopy( Q, index );
260 Teuchos::RCP<MV> tmpQj = MVT::CloneViewNonConst( Q, index );
261 index[0] = i;
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 );
265 }
266 }
267
268
270 // permuteVectors for MV
271 template<class ScalarType, class MV, class OP>
273 const std::vector<int> &perm,
274 Teuchos::SerialDenseMatrix<int,ScalarType> &Q)
275 {
276 // Permute the vectors in Q according to the permutation vector \c perm, and
277 // optionally the residual vector \c resids
278 Teuchos::BLAS<int,ScalarType> blas;
279 const int n = perm.size();
280 const int m = Q.numRows();
281
282 TEUCHOS_TEST_FOR_EXCEPTION(n != Q.numCols(), std::invalid_argument, "Anasazi::SolverUtils::permuteVectors(): size of permutation vector not equal to number of columns.");
283
284 // Sort the primitive ritz vectors
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);
288 }
289 }
290
291
292 //-----------------------------------------------------------------------------
293 //
294 // BASIS UPDATE METHODS
295 //
296 //-----------------------------------------------------------------------------
297
298 // apply householder reflectors to multivector
299 template<class ScalarType, class MV, class OP>
300 void SolverUtils<ScalarType, MV, OP>::applyHouse(int k, MV &V, const Teuchos::SerialDenseMatrix<int,ScalarType> &H, const std::vector<ScalarType> &tau, Teuchos::RCP<MV> workMV) {
301
302 const int n = MVT::GetNumberVecs(V);
303 const ScalarType ONE = SCT::one();
304 const ScalarType ZERO = SCT::zero();
305
306 // early exit if V has zero-size or if k==0
307 if (MVT::GetNumberVecs(V) == 0 || MVT::GetGlobalLength(V) == 0 || k == 0) {
308 return;
309 }
310
311 if (workMV == Teuchos::null) {
312 // user did not give us any workspace; allocate some
313 workMV = MVT::Clone(V,1);
314 }
315 else if (MVT::GetNumberVecs(*workMV) > 1) {
316 std::vector<int> first(1);
317 first[0] = 0;
318 workMV = MVT::CloneViewNonConst(*workMV,first);
319 }
320 else {
321 TEUCHOS_TEST_FOR_EXCEPTION(MVT::GetNumberVecs(*workMV) < 1,std::invalid_argument,"Anasazi::SolverUtils::applyHouse(): work multivector was empty.");
322 }
323 // Q = H_1 ... H_k is square, with as many rows as V has vectors
324 // however, H need only have k columns, one each for the k reflectors.
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.");
328
329 // perform the loop
330 // flops: Sum_{i=0:k-1} 4 m (n-i) == 4mnk - 2m(k^2- k)
331 for (int i=0; i<k; i++) {
332 // apply V H_i+1 = V - tau_i+1 (V v_i+1) v_i+1^T
333 // because of the structure of v_i+1, this transform does not affect the first i columns of V
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);
337
338 // note, below H_i, v_i and tau_i are mathematical objects which use 1-based indexing
339 // while H, v and tau are data structures using 0-based indexing
340
341 // get v_i+1: i-th column of H
342 Teuchos::SerialDenseMatrix<int,ScalarType> v(Teuchos::Copy,H,n-i,1,i,i);
343 // v_i+1(1:i) = 0: this isn't part of v
344 // e_i+1^T v_i+1 = 1 = v(0)
345 v(0,0) = ONE;
346
347 // compute -tau_i V v_i
348 // tau_i+1 is tau[i]
349 // flops: 2 m n-i
350 MVT::MvTimesMatAddMv(-tau[i],*actV,v,ZERO,*workMV);
351
352 // perform V = V + workMV v_i^T
353 // flops: 2 m n-i
354 Teuchos::SerialDenseMatrix<int,ScalarType> vT(v,Teuchos::CONJ_TRANS);
355 MVT::MvTimesMatAddMv(ONE,*workMV,vT,ONE,*actV);
356
357 actV = Teuchos::null;
358 }
359 }
360
361
362 //-----------------------------------------------------------------------------
363 //
364 // EIGENSOLVER PROJECTION METHODS
365 //
366 //-----------------------------------------------------------------------------
367
368 template<class ScalarType, class MV, class OP>
370 int size,
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)
376 {
377 // Routine for computing the first NEV generalized eigenpairs of the symmetric pencil (KK, MM)
378 //
379 // Parameter variables:
380 //
381 // size : Dimension of the eigenproblem (KK, MM)
382 //
383 // KK : Hermitian "stiffness" matrix
384 //
385 // MM : Hermitian positive-definite "mass" matrix
386 //
387 // EV : Matrix to store the nev eigenvectors
388 //
389 // theta : Array to store the eigenvalues (Size = nev )
390 //
391 // nev : Number of the smallest eigenvalues requested (input)
392 // Number of the smallest computed eigenvalues (output)
393 // Routine may compute and return more or less eigenvalues than requested.
394 //
395 // esType : Flag to select the algorithm
396 //
397 // esType = 0 (default) Uses LAPACK routine (Cholesky factorization of MM)
398 // with deflation of MM to get orthonormality of
399 // eigenvectors (S^T MM S = I)
400 //
401 // esType = 1 Uses LAPACK routine (Cholesky factorization of MM)
402 // (no check of orthonormality)
403 //
404 // esType = 10 Uses LAPACK routine for simple eigenproblem on KK
405 // (MM is not referenced in this case)
406 //
407 // Note: The code accesses only the upper triangular part of KK and MM.
408 //
409 // Return the integer info on the status of the computation
410 //
411 // info = 0 >> Success
412 //
413 // info < 0 >> error in the info-th argument
414 // info = - 20 >> Failure in LAPACK routine
415
416 // Define local arrays
417
418 // Create blas/lapack objects.
419 Teuchos::LAPACK<int,ScalarType> lapack;
420 Teuchos::BLAS<int,ScalarType> blas;
421
422 int rank = 0;
423 int info = 0;
424
425 if (size < nev || size < 0) {
426 return -1;
427 }
428 if (KK.numCols() < size || KK.numRows() < size) {
429 return -2;
430 }
431 if ((esType == 0 || esType == 1)) {
432 if (MM == Teuchos::null) {
433 return -3;
434 }
435 else if (MM->numCols() < size || MM->numRows() < size) {
436 return -3;
437 }
438 }
439 if (EV.numCols() < size || EV.numRows() < size) {
440 return -4;
441 }
442 if (theta.size() < (unsigned int) size) {
443 return -5;
444 }
445 if (nev <= 0) {
446 return -6;
447 }
448
449 // Query LAPACK for the "optimal" block size for HEGV
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); // For HEEV, lwork should be NB+2, instead of NB+1
454 std::vector<ScalarType> work(lwork);
455 std::vector<MagnitudeType> rwork(3*size-2);
456 // tt contains the eigenvalues from HEGV, which are necessarily real, and
457 // HEGV expects this vector to be real as well
458 std::vector<MagnitudeType> tt( size );
459 //typedef typename std::vector<MagnitudeType>::iterator MTIter; // unused
460
461 MagnitudeType tol = SCT::magnitude(SCT::squareroot(SCT::eps()));
462 // MagnitudeType tol = 1e-12;
463 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
464 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
465
466 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > KKcopy, MMcopy;
467 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > U;
468
469 switch (esType) {
470 default:
471 case 0:
472 //
473 // Use LAPACK to compute the generalized eigenvectors
474 //
475 for (rank = size; rank > 0; --rank) {
476
477 U = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(rank,rank) );
478 //
479 // Copy KK & MM
480 //
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 ) );
483 //
484 // Solve the generalized eigenproblem with LAPACK
485 //
486 info = 0;
487 lapack.HEGV(1, 'V', 'U', rank, KKcopy->values(), KKcopy->stride(),
488 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
489 &rwork[0], &info);
490 //
491 // Treat error messages
492 //
493 if (info < 0) {
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;
497 return -20;
498 }
499 if (info > 0) {
500 if (info > rank)
501 rank = info - rank;
502 continue;
503 }
504 //
505 // Check the quality of eigenvectors ( using mass-orthonormality )
506 //
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));
511 }
512 }
513 // U = 0*U + 1*MMcopy*KKcopy = MMcopy * KKcopy
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.");
517 // MMcopy = 0*MMcopy + 1*KKcopy^H*U = KKcopy^H * MMcopy * KKcopy
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) {
525 if (j == i)
526 maxNorm = SCT::magnitude((*MMcopy)(i,j) - one) > maxNorm
527 ? SCT::magnitude((*MMcopy)(i,j) - one) : maxNorm;
528 else
529 maxOrth = SCT::magnitude((*MMcopy)(i,j)) > maxOrth
530 ? SCT::magnitude((*MMcopy)(i,j)) : maxOrth;
531 }
532 }
533 /* if (verbose > 4) {
534 std::cout << " >> Local eigensolve >> Size: " << rank;
535 std::cout.precision(2);
536 std::cout.setf(std::ios::scientific, std::ios::floatfield);
537 std::cout << " Normalization error: " << maxNorm;
538 std::cout << " Orthogonality error: " << maxOrth;
539 std::cout << endl;
540 }*/
541 if ((maxNorm <= tol) && (maxOrth <= tol)) {
542 break;
543 }
544 } // for (rank = size; rank > 0; --rank)
545 //
546 // Copy the computed eigenvectors and eigenvalues
547 // ( they may be less than the number requested because of deflation )
548 //
549 // std::cout << "directSolve rank: " << rank << "\tsize: " << size << endl;
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 );
555 }
556 break;
557
558 case 1:
559 //
560 // Use the Cholesky factorization of MM to compute the generalized eigenvectors
561 //
562 // Copy KK & MM
563 //
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 ) );
566 //
567 // Solve the generalized eigenproblem with LAPACK
568 //
569 info = 0;
570 lapack.HEGV(1, 'V', 'U', size, KKcopy->values(), KKcopy->stride(),
571 MMcopy->values(), MMcopy->stride(), &tt[0], &work[0], lwork,
572 &rwork[0], &info);
573 //
574 // Treat error messages
575 //
576 if (info < 0) {
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;
580 return -20;
581 }
582 if (info > 0) {
583 if (info > size)
584 nev = 0;
585 else {
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;
589 return -20;
590 }
591 }
592 //
593 // Copy the eigenvectors and eigenvalues
594 //
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 );
598 }
599 break;
600
601 case 10:
602 //
603 // Simple eigenproblem
604 //
605 // Copy KK
606 //
607 KKcopy = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, KK, size, size ) );
608 //
609 // Solve the generalized eigenproblem with LAPACK
610 //
611 lapack.HEEV('V', 'U', size, KKcopy->values(), KKcopy->stride(), &tt[0], &work[0], lwork, &rwork[0], &info);
612 //
613 // Treat error messages
614 if (info != 0) {
615 std::cerr << std::endl;
616 if (info < 0) {
617 std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, argument " << -info << " has an illegal value\n";
618 }
619 else {
620 std::cerr << "Anasazi::SolverUtils::directSolver(): In DHEEV, the algorithm failed to converge (" << info << ").\n";
621 }
622 std::cerr << std::endl;
623 info = -20;
624 break;
625 }
626 //
627 // Copy the eigenvectors
628 //
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 );
632 }
633 break;
634 }
635
636 return info;
637 }
638
639
640 //-----------------------------------------------------------------------------
641 //
642 // SANITY CHECKING METHODS
643 //
644 //-----------------------------------------------------------------------------
645
646 template<class ScalarType, class MV, class OP>
647 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
648 SolverUtils<ScalarType, MV, OP>::errorEquality(const MV &X, const MV &MX, Teuchos::RCP<const OP> M)
649 {
650 // Return the maximum coefficient of the matrix M * X - MX
651 // scaled by the maximum coefficient of MX.
652 // When M is not specified, the identity is used.
653
654 MagnitudeType maxDiff = SCT::magnitude(SCT::zero());
655
656 int xc = MVT::GetNumberVecs(X);
657 int mxc = MVT::GetNumberVecs(MX);
658
659 TEUCHOS_TEST_FOR_EXCEPTION(xc != mxc,std::invalid_argument,"Anasazi::SolverUtils::errorEquality(): input multivecs have different number of columns.");
660 if (xc == 0) {
661 return maxDiff;
662 }
663
664 MagnitudeType maxCoeffX = SCT::magnitude(SCT::zero());
665 std::vector<MagnitudeType> tmp( xc );
666 MVT::MvNorm(MX, tmp);
667
668 for (int i = 0; i < xc; ++i) {
669 maxCoeffX = (tmp[i] > maxCoeffX) ? tmp[i] : maxCoeffX;
670 }
671
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 );
677 }
678 else {
679 MtimesX = MVT::CloneCopy(X);
680 }
681 MVT::MvAddMv( -1.0, MX, 1.0, *MtimesX, *MtimesX );
682 MVT::MvNorm( *MtimesX, tmp );
683
684 for (int i = 0; i < xc; ++i) {
685 maxDiff = (tmp[i] > maxDiff) ? tmp[i] : maxDiff;
686 }
687
688 return (maxCoeffX == 0.0) ? maxDiff : maxDiff/maxCoeffX;
689
690 }
691
692} // end namespace Anasazi
693
694#endif // ANASAZI_SOLVER_UTILS_HPP
695
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Abstract class definition for Anasazi Output Managers.
Anasazi's templated, static class providing utilities for the solvers.
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::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.
virtual ~SolverUtils()
Destructor.
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,...
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.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.