Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziMVOPTester.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Anasazi: Block Eigensolvers Package
4//
5// Copyright 2004 NTESS and the Anasazi contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9//
10#ifndef ANASAZI_MVOPTESTER_HPP
11#define ANASAZI_MVOPTESTER_HPP
12
13// Assumptions that I have made:
14// * I assume/verify that a multivector must have at least one vector. This seems
15// to be consistent with Epetra_MultiVec.
16// * I do not assume that an operator is deterministic; I do assume that the
17// operator, applied to 0, will return 0.
18
27#include "AnasaziConfigDefs.hpp"
28#include "AnasaziTypes.hpp"
29
30#include "AnasaziMultiVecTraits.hpp"
31#include "AnasaziOperatorTraits.hpp"
32#include "AnasaziOutputManager.hpp"
33
34#include "Teuchos_SetScientific.hpp"
35#include "Teuchos_RCP.hpp"
36#include "Teuchos_as.hpp"
37
38namespace Anasazi {
39
45 template< class ScalarType, class MV >
46 bool TestMultiVecTraits(
47 const Teuchos::RCP<OutputManager<ScalarType> > &om,
48 const Teuchos::RCP<const MV> &A ) {
49
50 using std::endl;
51 using Teuchos::SetScientific;
52
53 /* MVT Contract:
54
55 Clone(MV,int)
56 CloneCopy(MV)
57 CloneCopy(MV,vector<int>)
58 USER: will request positive number of vectors
59 MV: will return a multivector with exactly the number of
60 requested vectors.
61 vectors are the same dimension as the cloned MV
62
63
64 CloneView(MV,vector<int>) [const and non-const]
65 USER: There is no assumed communication between creation and
66 destruction of a view. I.e., after a view is created, changes to the
67 source multivector are not reflected in the view. Likewise, until
68 destruction of the view, changes in the view are not reflected in the
69 source multivector.
70
71 GetGlobalLength
72 MV: will always be positive (MV cannot have zero vectors)
73
74 GetNumberVecs
75 MV: will always be positive (MV cannot have zero vectors)
76
77 MvAddMv
78 USER: multivecs will be of the same dimension and same number of vecs
79 MV: input vectors will not be modified
80 performing C=0*A+1*B will assign B to C exactly
81
82 MvTimesMatAddMv
83 USER: multivecs and serialdensematrix will be of the proper shape
84 MV: input arguments will not be modified
85 following arithmetic relations hold exactly:
86 A*I = A
87 0*B = B
88 1*B = B
89
90 MvTransMv
91 USER: SerialDenseMatrix will be large enough to hold results.
92 MV: SerialDenseMatrix will not be resized.
93 Inner products will satisfy |a'*b| <= |a|*|b|
94 alpha == 0 => SerialDenseMatrix == 0
95
96 MvDot
97 USER: Results vector will be large enough for results.
98 Both multivectors will have the same number of vectors.
99 (Epetra crashes, otherwise.)
100 MV: Inner products will satisfy |a'*b| <= |a|*|b|
101 Results vector will not be resized.
102
103 MvNorm
104 MV: vector norm is always non-negative, and zero
105 only for zero vectors.
106 results vector should not be resized
107
108 SetBlock
109 USER: indices will be distinct
110 MV: assigns copies of the vectors to the specified
111 locations, leaving the other vectors untouched.
112
113 MvRandom
114 MV: Generate zero vector with "zero" probability
115 Don't gen the same vectors twice.
116
117 MvInit
118 MV: Init(alpha) sets all elements to alpha
119
120 MvScale (two versions)
121 MV: scales multivector values
122
123 MvPrint
124 MV: routine does not modify vectors (not tested here)
125 *********************************************************************/
126
127 typedef MultiVecTraits<ScalarType, MV> MVT;
128 typedef Teuchos::ScalarTraits<ScalarType> SCT;
129 typedef typename SCT::magnitudeType MagType;
130
131 const ScalarType one = SCT::one();
132 const ScalarType zero = SCT::zero();
133 const MagType zero_mag = Teuchos::ScalarTraits<MagType>::zero();
134 const MagType tol = SCT::eps()*100;
135
136 // Don't change these two without checking the initialization of ind below
137 const int numvecs = 10;
138 const int numvecs_2 = 5;
139
140 std::vector<int> ind(numvecs_2);
141
142 /* Initialize indices for selected copies/views
143 The MVT specialization should not assume that
144 these are ordered or even distinct.
145 Also retrieve the edges.
146
147 However, to spice things up, grab the first vector,
148 last vector, and choose the others randomly.
149 */
150 TEUCHOS_TEST_FOR_EXCEPT(numvecs_2 != 5);
151 ind[0] = 0;
152 ind[1] = 5;
153 ind[2] = 2;
154 ind[3] = 2;
155 ind[4] = 9;
156
157 /*********** GetNumberVecs() *****************************************
158 Verify:
159 1) This number should be strictly positive
160 *********************************************************************/
161 if ( MVT::GetNumberVecs(*A) <= 0 ) {
162 om->stream(Warnings)
163 << "*** ERROR *** MultiVectorTraits::GetNumberVecs()." << endl
164 << "Returned <= 0." << endl;
165 return false;
166 }
167
168 /*********** GetGlobalLength() ***************************************
169 Verify:
170 1) This number should be strictly positive
171 *********************************************************************/
172 if ( MVT::GetGlobalLength(*A) <= 0 ) {
173 om->stream(Warnings)
174 << "*** ERROR *** MultiVectorTraitsExt::GetGlobalLength()" << endl
175 << "Returned <= 0." << endl;
176 return false;
177 }
178
179 /*********** Clone() and MvNorm() ************************************
180 Verify:
181 1) Clone() allows us to specify the number of vectors
182 2) Clone() returns a multivector of the same dimension
183 3) Vector norms shouldn't be negative
184 *********************************************************************/
185 {
186 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
187 MVT::MvInit(*B);
188 std::vector<MagType> norms(2*numvecs);
189 bool ResizeWarning = false;
190 if ( MVT::GetNumberVecs(*B) != numvecs ) {
191 om->stream(Warnings)
192 << "*** ERROR *** MultiVecTraits::Clone()." << endl
193 << "Did not allocate requested number of vectors." << endl;
194 return false;
195 }
196 if ( MVT::GetGlobalLength(*B) != MVT::GetGlobalLength(*A) ) {
197 om->stream(Warnings)
198 << "*** ERROR *** MultiVecTraits::Clone()." << endl
199 << "Did not allocate requested number of vectors." << endl;
200 return false;
201 }
202 MVT::MvNorm(*B, norms);
203 if ( (int)norms.size() != 2*numvecs && (ResizeWarning == false) ) {
204 om->stream(Warnings)
205 << "*** WARNING *** MultiVecTraits::MvNorm()." << endl
206 << "Method resized the output vector." << endl;
207 ResizeWarning = true;
208 }
209 for (int i=0; i<numvecs; i++) {
210 if ( norms[i] < zero_mag ) {
211 om->stream(Warnings)
212 << "*** ERROR *** MultiVecTraits::Clone()." << endl
213 << "Vector had negative norm." << endl;
214 return false;
215 }
216 }
217 }
218
219
220 /*********** MvRandom() and MvNorm() and MvInit() ********************
221 Verify:
222 1) Set vectors to zero
223 2) Check that norm is zero
224 3) Perform MvRandom.
225 4) Verify that vectors aren't zero anymore
226 5) Perform MvRandom again.
227 6) Verify that vector norms are different than before
228
229 Without knowing something about the random distribution,
230 this is about the best that we can do, to make sure that MvRandom
231 did at least *something*.
232
233 Also, make sure vector norms aren't negative.
234 *********************************************************************/
235 {
236 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
237 std::vector<MagType> norms(numvecs), norms2(numvecs);
238
239 MVT::MvInit(*B);
240 MVT::MvNorm(*B, norms);
241 for (int i=0; i<numvecs; i++) {
242 if ( norms[i] != zero_mag ) {
243 om->stream(Warnings)
244 << "*** ERROR *** MultiVecTraits::MvInit() "
245 << "and MultiVecTraits::MvNorm()" << endl
246 << "Supposedly zero vector has non-zero norm." << endl;
247 return false;
248 }
249 }
250 MVT::MvRandom(*B);
251 MVT::MvNorm(*B, norms);
252 MVT::MvRandom(*B);
253 MVT::MvNorm(*B, norms2);
254 for (int i=0; i<numvecs; i++) {
255 if ( norms[i] == zero_mag || norms2[i] == zero_mag ) {
256 om->stream(Warnings)
257 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl
258 << "Random vector was empty (very unlikely)." << endl;
259 return false;
260 }
261 else if ( norms[i] < zero_mag || norms2[i] < zero_mag ) {
262 om->stream(Warnings)
263 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl
264 << "Vector had negative norm." << endl;
265 return false;
266 }
267 else if ( norms[i] == norms2[i] ) {
268 om->stream(Warnings)
269 << "*** ERROR *** MutliVecTraits::MvRandom()." << endl
270 << "Vectors not random enough." << endl;
271 return false;
272 }
273 }
274 }
275
276
277 /*********** MvRandom() and MvNorm() and MvScale() *******************
278 Verify:
279 1) Perform MvRandom.
280 2) Verify that vectors aren't zero
281 3) Set vectors to zero via MvScale
282 4) Check that norm is zero
283 *********************************************************************/
284 {
285 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
286 std::vector<MagType> norms(numvecs);
287
288 MVT::MvRandom(*B);
289 MVT::MvScale(*B,SCT::zero());
290 MVT::MvNorm(*B, norms);
291 for (int i=0; i<numvecs; i++) {
292 if ( norms[i] != zero_mag ) {
293 om->stream(Warnings)
294 << "*** ERROR *** MultiVecTraits::MvScale(alpha) "
295 << "Supposedly zero vector has non-zero norm." << endl;
296 return false;
297 }
298 }
299
300 MVT::MvRandom(*B);
301 std::vector<ScalarType> zeros(numvecs,SCT::zero());
302 MVT::MvScale(*B,zeros);
303 MVT::MvNorm(*B, norms);
304 for (int i=0; i<numvecs; i++) {
305 if ( norms[i] != zero_mag ) {
306 om->stream(Warnings)
307 << "*** ERROR *** MultiVecTraits::MvScale(alphas) "
308 << "Supposedly zero vector has non-zero norm." << endl;
309 return false;
310 }
311 }
312 }
313
314
315 /*********** MvInit() and MvNorm() ***********************************
316 A vector of ones of dimension n should have norm sqrt(n)
317 1) Init vectors to all ones
318 2) Verify that norm is sqrt(n)
319 3) Verify that norms aren't negative
320
321 Note: I'm not sure that we can expect this to hold in practice.
322 Maybe something like abs(norm-sqrt(n)) < SCT::eps() ???
323 The sum of 1^2==1 should be n, but what about sqrt(n)?
324 They may be using a different square root than ScalartTraits
325 On my iBook G4 and on jeter, this test works.
326 Right now, this has been demoted to a warning.
327 *********************************************************************/
328 {
329 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
330 std::vector<MagType> norms(numvecs);
331
332 MVT::MvInit(*B,one);
333 MVT::MvNorm(*B, norms);
334 bool BadNormWarning = false;
335 for (int i=0; i<numvecs; i++) {
336 if ( norms[i] < zero_mag ) {
337 om->stream(Warnings)
338 << "*** ERROR *** MultiVecTraits::MvRandom()." << endl
339 << "Vector had negative norm." << endl;
340 return false;
341 }
342 else if ( norms[i] != SCT::squareroot(MVT::GetGlobalLength(*B)) && !BadNormWarning ) {
343 om->stream(Warnings)
344 << endl
345 << "Warning testing MultiVecTraits::MvInit()." << endl
346 << "Ones vector should have norm sqrt(dim)." << endl
347 << "norms[i]: " << norms[i] << "\tdim: " << MVT::GetGlobalLength(*B) << endl << endl;
348 BadNormWarning = true;
349 }
350 }
351 }
352
353
354 /*********** MvInit() and MvNorm() ***********************************
355 A vector of zeros of dimension n should have norm 0
356 1) Verify that norms aren't negative
357 2) Verify that norms are zero
358
359 We must know this works before the next tests.
360 *********************************************************************/
361 {
362 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs);
363 std::vector<MagType> norms(numvecs);
364 MVT::MvInit(*B, zero_mag);
365 MVT::MvNorm(*B, norms);
366 for (int i=0; i<numvecs; i++) {
367 if ( norms[i] < zero_mag ) {
368 om->stream(Warnings)
369 << "*** ERROR *** MultiVecTraits::MvInit()." << endl
370 << "Vector had negative norm." << endl;
371 return false;
372 }
373 else if ( norms[i] != zero_mag ) {
374 om->stream(Warnings)
375 << "*** ERROR *** MultiVecTraits::MvInit()." << endl
376 << "Zero vector should have norm zero." << endl;
377 return false;
378 }
379 }
380 }
381
382
383 /*********** CloneCopy(MV,vector<int>) and MvNorm ********************
384 1) Check quantity/length of vectors
385 2) Check vector norms for agreement
386 3) Zero out B and make sure that C norms are not affected
387 *********************************************************************/
388 {
389 Teuchos::RCP<MV> B, C;
390 std::vector<MagType> norms(numvecs), norms2(ind.size());
391
392 B = MVT::Clone(*A,numvecs);
393 MVT::MvRandom(*B);
394 MVT::MvNorm(*B, norms);
395 C = MVT::CloneCopy(*B,ind);
396 MVT::MvNorm(*C, norms2);
397 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
398 om->stream(Warnings)
399 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
400 << "Wrong number of vectors." << endl;
401 return false;
402 }
403 if ( MVT::GetGlobalLength(*C) != MVT::GetGlobalLength(*B) ) {
404 om->stream(Warnings)
405 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
406 << "Vector lengths don't match." << endl;
407 return false;
408 }
409 for (int i=0; i<numvecs_2; i++) {
410 if ( SCT::magnitude( norms2[i] - norms[ind[i]] ) > tol ) {
411 om->stream(Warnings)
412 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
413 << "Copied vectors do not agree:"
414 << norms2[i] << " != " << norms[ind[i]] << endl
415 << "Difference " << SCT::magnitude (norms2[i] - norms[ind[i]])
416 << " exceeds the tolerance 100*eps = " << tol << endl;
417
418 return false;
419 }
420 }
421 MVT::MvInit(*B,zero);
422 MVT::MvNorm(*C, norms2);
423 for (int i=0; i<numvecs_2; i++) {
424 if ( SCT::magnitude( norms2[i] - norms2[i] ) > tol ) {
425 om->stream(Warnings)
426 << "*** ERROR *** MultiVecTraits::CloneCopy(ind)." << endl
427 << "Copied vectors were not independent." << endl
428 << "Difference " << SCT::magnitude (norms2[i] - norms[i])
429 << " exceeds the tolerance 100*eps = " << tol << endl;
430 return false;
431 }
432 }
433 }
434
435
436 /*********** CloneCopy(MV) and MvNorm ********************************
437 1) Check quantity
438 2) Check value of norms
439 3) Zero out B and make sure that C is still okay
440 *********************************************************************/
441 {
442 Teuchos::RCP<MV> B, C;
443 std::vector<MagType> norms(numvecs), norms2(numvecs);
444
445 B = MVT::Clone(*A,numvecs);
446 MVT::MvRandom(*B);
447 MVT::MvNorm(*B, norms);
448 C = MVT::CloneCopy(*B);
449 MVT::MvNorm(*C, norms2);
450 if ( MVT::GetNumberVecs(*C) != numvecs ) {
451 om->stream(Warnings)
452 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
453 << "Wrong number of vectors." << endl;
454 return false;
455 }
456 for (int i=0; i<numvecs; i++) {
457 if ( SCT::magnitude( norms2[i] - norms[i] ) > tol ) {
458 om->stream(Warnings)
459 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
460 << "Copied vectors do not agree." << endl
461 << "Difference " << SCT::magnitude (norms2[i] - norms[i])
462 << " exceeds the tolerance 100*eps = " << tol << endl;
463 return false;
464 }
465 }
466 MVT::MvInit(*B,zero);
467 MVT::MvNorm(*C, norms);
468 for (int i=0; i<numvecs; i++) {
469 if ( SCT::magnitude( norms2[i] - norms[i] ) > tol ) {
470 om->stream(Warnings)
471 << "*** ERROR *** MultiVecTraits::CloneCopy()." << endl
472 << "Copied vectors were not independent." << endl
473 << "Difference " << SCT::magnitude (norms2[i] - norms[i])
474 << " exceeds the tolerance 100*eps = " << tol << endl;
475 return false;
476 }
477 }
478 }
479
480
481 /*********** CloneView(MV,vector<int>) and MvNorm ********************
482 Check that we have a view of the selected vectors
483 1) Check quantity
484 2) Check value of norms
485 3) Zero out B and make sure that C is zero as well
486 *********************************************************************/
487 {
488 Teuchos::RCP<MV> B, C;
489 std::vector<MagType> norms(numvecs), norms2(ind.size());
490
491 B = MVT::Clone(*A,numvecs);
492 MVT::MvRandom(*B);
493 MVT::MvNorm(*B, norms);
494 C = MVT::CloneViewNonConst(*B,ind);
495 MVT::MvNorm(*C, norms2);
496 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
497 om->stream(Warnings)
498 << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
499 << "Wrong number of vectors." << endl;
500 return false;
501 }
502 for (int i=0; i<numvecs_2; i++) {
503 if ( SCT::magnitude( norms2[i] - norms[ind[i]] ) > tol ) {
504 om->stream(Warnings)
505 << "*** ERROR *** MultiVecTraits::CloneView(ind)." << endl
506 << "Viewed vectors do not agree." << endl;
507 return false;
508 }
509 }
510 }
511
512
513 /*********** const CloneView(MV,vector<int>) and MvNorm() ************
514 Check that we have a view of the selected vectors.
515 1) Check quantity
516 2) Check value of norms for agreement
517 3) Zero out B and make sure that C is zerod as well
518 *********************************************************************/
519 {
520 Teuchos::RCP<MV> B;
521 Teuchos::RCP<const MV> constB, C;
522 std::vector<MagType> normsB(numvecs), normsC(ind.size());
523 std::vector<int> allind(numvecs);
524 for (int i=0; i<numvecs; i++) {
525 allind[i] = i;
526 }
527
528 B = MVT::Clone(*A,numvecs);
529 MVT::MvRandom( *B );
530 MVT::MvNorm(*B, normsB);
531 C = MVT::CloneView(*B,ind);
532 MVT::MvNorm(*C, normsC);
533 if ( MVT::GetNumberVecs(*C) != numvecs_2 ) {
534 om->stream(Warnings)
535 << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
536 << "Wrong number of vectors." << endl;
537 return false;
538 }
539 for (int i=0; i<numvecs_2; i++) {
540 if ( SCT::magnitude( normsC[i] - normsB[ind[i]] ) > tol ) {
541 om->stream(Warnings)
542 << "*** ERROR *** const MultiVecTraits::CloneView(ind)." << endl
543 << "Viewed vectors do not agree." << endl;
544 return false;
545 }
546 }
547 }
548
549
550 /*********** SetBlock() and MvNorm() *********************************
551 SetBlock() will copy the vectors from C into B
552 1) Verify that the specified vectors were copied
553 2) Verify that the other vectors were not modified
554 3) Verify that C was not modified
555 4) Change C and then check B to make sure it was not modified
556
557 Use a different index set than has been used so far (distinct entries).
558 This is because duplicate entries will cause the vector to be
559 overwritten, making it more difficult to test.
560 *********************************************************************/
561 {
562 Teuchos::RCP<MV> B, C;
563 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
564 normsC1(numvecs_2), normsC2(numvecs_2);
565
566 B = MVT::Clone(*A,numvecs);
567 C = MVT::Clone(*A,numvecs_2);
568 // Just do every other one, interleaving the vectors of C into B
569 ind.resize(numvecs_2);
570 for (int i=0; i<numvecs_2; i++) {
571 ind[i] = 2*i;
572 }
573 MVT::MvRandom(*B);
574 MVT::MvRandom(*C);
575
576 MVT::MvNorm(*B,normsB1);
577 MVT::MvNorm(*C,normsC1);
578 MVT::SetBlock(*C,ind,*B);
579 MVT::MvNorm(*B,normsB2);
580 MVT::MvNorm(*C,normsC2);
581
582 // check that C was not changed by SetBlock
583 for (int i=0; i<numvecs_2; i++) {
584 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
585 om->stream(Warnings)
586 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
587 << "Operation modified source vectors." << endl;
588 return false;
589 }
590 }
591 // check that the correct vectors of B were modified
592 // and the others were not
593 for (int i=0; i<numvecs; i++) {
594 if (i % 2 == 0) {
595 // should be a vector from C
596 if ( SCT::magnitude(normsB2[i]-normsC1[i/2]) > tol ) {
597 om->stream(Warnings)
598 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
599 << "Copied vectors do not agree." << endl
600 << "Difference " << SCT::magnitude (normsB2[i] - normsC1[i/2])
601 << " exceeds the tolerance 100*eps = " << tol << endl;
602 return false;
603 }
604 }
605 else {
606 // should be an original vector
607 if ( SCT::magnitude(normsB1[i]-normsB2[i]) > tol ) {
608 om->stream(Warnings)
609 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
610 << "Incorrect vectors were modified." << endl;
611 return false;
612 }
613 }
614 }
615 MVT::MvInit(*C,zero);
616 MVT::MvNorm(*B,normsB1);
617 // verify that we copied and didn't reference
618 for (int i=0; i<numvecs; i++) {
619 if ( SCT::magnitude(normsB1[i]-normsB2[i]) > tol ) {
620 om->stream(Warnings)
621 << "*** ERROR *** MultiVecTraits::SetBlock()." << endl
622 << "Copied vectors were not independent." << endl;
623 return false;
624 }
625 }
626 }
627
628
629 /*********** MvTransMv() *********************************************
630 Performs C = alpha * A^H * B, where
631 alpha is type ScalarType
632 A,B are type MV with p and q vectors, respectively
633 C is a SerialDenseMatrix<int,ScalarType> ALREADY sized to p by q
634
635 Verify:
636 1) C is not resized by the routine
637 3) Check that zero*(A^H B) == zero
638 3) Check inner product inequality:
639 [ |a1|*|b1| ... |ap|*|b1| ]
640 [a1 ... ap]^H [b1 ... bq] <= [ ... |ai|*|bj| ... ]
641 [ |ap|*|b1| ... |ap|*|bq| ]
642 4) Zero B and check that B^H * C is zero
643 5) Zero C and check that B^H * C is zero
644
645 Note 1: Test 4 is performed with a p x q Teuchos::SDM view of
646 a (p+1) x (q+1) Teuchos::SDM that is initialized to ones.
647 This ensures the user is correctly accessing and filling the SDM.
648
649 Note 2: Should we really require that C is correctly sized already?
650 Epetra does (and crashes if it isn't.)
651 *********************************************************************/
652 {
653 const int p = 7;
654 const int q = 9;
655 Teuchos::RCP<MV> B, C;
656 std::vector<MagType> normsB(p), normsC(q);
657 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
658
659 B = MVT::Clone(*A,p);
660 C = MVT::Clone(*A,q);
661
662 // randomize the multivectors
663 MVT::MvRandom(*B);
664 MVT::MvNorm(*B,normsB);
665 MVT::MvRandom(*C);
666 MVT::MvNorm(*C,normsC);
667
668 // perform SDM = zero() * B^H * C
669 MVT::MvTransMv( zero, *B, *C, SDM );
670
671 // check the sizes: not allowed to have shrunk
672 if ( SDM.numRows() != p || SDM.numCols() != q ) {
673 om->stream(Warnings)
674 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
675 << "Routine resized SerialDenseMatrix." << endl;
676 return false;
677 }
678
679 // check that zero**A^H*B == zero
680 if ( SDM.normOne() != zero ) {
681 om->stream(Warnings)
682 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
683 << "Scalar argument processed incorrectly." << endl;
684 return false;
685 }
686
687 // perform SDM = one * B^H * C
688 MVT::MvTransMv( one, *B, *C, SDM );
689
690 // check the norms: a^H b = |a| |b| cos(theta) <= |a| |b|
691 // with equality only when a and b are colinear
692 for (int i=0; i<p; i++) {
693 for (int j=0; j<q; j++) {
694 if ( SCT::magnitude(SDM(i,j))
695 > SCT::magnitude(normsB[i]*normsC[j]) ) {
696 om->stream(Warnings)
697 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
698 << "Triangle inequality did not hold: "
699 << SCT::magnitude(SDM(i,j))
700 << " > "
701 << SCT::magnitude(normsB[i]*normsC[j])
702 << endl;
703 return false;
704 }
705 }
706 }
707 MVT::MvInit(*C);
708 MVT::MvRandom(*B);
709 MVT::MvTransMv( one, *B, *C, SDM );
710 for (int i=0; i<p; i++) {
711 for (int j=0; j<q; j++) {
712 if ( SDM(i,j) != zero ) {
713 om->stream(Warnings)
714 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
715 << "Inner products not zero for C==0." << endl;
716 return false;
717 }
718 }
719 }
720 MVT::MvInit(*B);
721 MVT::MvRandom(*C);
722 MVT::MvTransMv( one, *B, *C, SDM );
723 for (int i=0; i<p; i++) {
724 for (int j=0; j<q; j++) {
725 if ( SDM(i,j) != zero ) {
726 om->stream(Warnings)
727 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
728 << "Inner products not zero for B==0." << endl;
729 return false;
730 }
731 }
732 }
733
734 // A larger SDM is filled with ones, initially, and a smaller
735 // view is used for the MvTransMv method. If the smaller SDM
736 // is not all zeroes, then the interface is improperly writing
737 // to the matrix object.
738 // Note: Since we didn't fail above, we know that the general
739 // inner product works, but we are checking to see if it
740 // works for a view too. This is common usage in Anasazi.
741 Teuchos::SerialDenseMatrix<int, ScalarType> largeSDM(p+1,q+1);
742 Teuchos::SerialDenseMatrix<int, ScalarType> SDM2(Teuchos::View, largeSDM, p, q);
743 largeSDM.putScalar( one );
744 MVT::MvInit(*C);
745 MVT::MvRandom(*B);
746 MVT::MvTransMv( one, *B, *C, SDM2 );
747 for (int i=0; i<p; i++) {
748 for (int j=0; j<q; j++) {
749 if ( SDM2(i,j) != zero ) {
750 om->stream(Warnings)
751 << "*** ERROR *** MultiVecTraits::MvTransMv()." << endl
752 << "Inner products not zero for C==0 when using a view into Teuchos::SerialDenseMatrix<>." << endl;
753 return false;
754 }
755 }
756 }
757 }
758
759
760 /*********** MvDot() *************************************************
761 Verify:
762 1) Results vector not resized
763 2) Inner product inequalities are satisfied
764 3) Zero vectors give zero inner products
765 *********************************************************************/
766 {
767 const int p = 7;
768 const int q = 9;
769 Teuchos::RCP<MV> B, C;
770 std::vector<ScalarType> iprods(q);
771 std::vector<MagType> normsB(p), normsC(p);
772
773 B = MVT::Clone(*A,p);
774 C = MVT::Clone(*A,p);
775
776 MVT::MvRandom(*B);
777 MVT::MvRandom(*C);
778 MVT::MvNorm(*B,normsB);
779 MVT::MvNorm(*C,normsC);
780 MVT::MvDot( *B, *C, iprods );
781 if ( (int)iprods.size() != q ) {
782 om->stream(Warnings)
783 << "*** ERROR *** MultiVecTraits::MvDot." << endl
784 << "Routine resized results vector." << endl;
785 return false;
786 }
787 for (int i=0; i<p; i++) {
788 if ( SCT::magnitude(iprods[i])
789 > SCT::magnitude(normsB[i]*normsC[i]) ) {
790 om->stream(Warnings)
791 << "*** ERROR *** MultiVecTraits::MvDot()." << endl
792 << "Inner products not valid." << endl;
793 return false;
794 }
795 }
796 MVT::MvInit(*B);
797 MVT::MvRandom(*C);
798 MVT::MvDot( *B, *C, iprods );
799 for (int i=0; i<p; i++) {
800 if ( iprods[i] != zero ) {
801 om->stream(Warnings)
802 << "*** ERROR *** MultiVecTraits::MvDot()." << endl
803 << "Inner products not zero for B==0." << endl;
804 return false;
805 }
806 }
807 MVT::MvInit(*C);
808 MVT::MvRandom(*B);
809 MVT::MvDot( *B, *C, iprods );
810 for (int i=0; i<p; i++) {
811 if ( iprods[i] != zero ) {
812 om->stream(Warnings)
813 << "*** ERROR *** MultiVecTraits::MvDot()." << endl
814 << "Inner products not zero for C==0." << endl;
815 return false;
816 }
817 }
818 }
819
820
821 /*********** MvAddMv() ***********************************************
822 D = alpha*B + beta*C
823 1) Use alpha==0,beta==1 and check that D == C
824 2) Use alpha==1,beta==0 and check that D == B
825 3) Use D==0 and D!=0 and check that result is the same
826 4) Check that input arguments are not modified
827 *********************************************************************/
828 {
829 const int p = 7;
830 Teuchos::RCP<MV> B, C, D;
831 std::vector<MagType> normsB1(p), normsB2(p),
832 normsC1(p), normsC2(p),
833 normsD1(p), normsD2(p);
834 ScalarType alpha = MagType(0.5) * SCT::one();
835 ScalarType beta = MagType(0.33) * SCT::one();
836
837 B = MVT::Clone(*A,p);
838 C = MVT::Clone(*A,p);
839 D = MVT::Clone(*A,p);
840
841 MVT::MvRandom(*B);
842 MVT::MvRandom(*C);
843 MVT::MvNorm(*B,normsB1);
844 MVT::MvNorm(*C,normsC1);
845
846 // check that 0*B+1*C == C
847 MVT::MvAddMv(zero,*B,one,*C,*D);
848 MVT::MvNorm(*B,normsB2);
849 MVT::MvNorm(*C,normsC2);
850 MVT::MvNorm(*D,normsD1);
851 for (int i=0; i<p; i++) {
852 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
853 om->stream(Warnings)
854 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
855 << "Input arguments were modified." << endl;
856 return false;
857 }
858 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
859 om->stream(Warnings)
860 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
861 << "Input arguments were modified." << endl;
862 return false;
863 }
864 else if ( SCT::magnitude(normsC1[i]-normsD1[i]) > tol ) {
865 om->stream(Warnings)
866 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
867 << "Assignment did not work." << endl;
868 return false;
869 }
870 }
871
872 // check that 1*B+0*C == B
873 MVT::MvAddMv(one,*B,zero,*C,*D);
874 MVT::MvNorm(*B,normsB2);
875 MVT::MvNorm(*C,normsC2);
876 MVT::MvNorm(*D,normsD1);
877 for (int i=0; i<p; i++) {
878 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
879 om->stream(Warnings)
880 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
881 << "Input arguments were modified." << endl;
882 return false;
883 }
884 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
885 om->stream(Warnings)
886 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
887 << "Input arguments were modified." << endl;
888 return false;
889 }
890 else if ( SCT::magnitude( normsB1[i] - normsD1[i] ) > tol ) {
891 om->stream(Warnings)
892 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
893 << "Assignment did not work." << endl;
894 return false;
895 }
896 }
897
898 // check that alpha*B+beta*C -> D is invariant under initial D
899 // first, try random D
900 MVT::MvRandom(*D);
901 MVT::MvAddMv(alpha,*B,beta,*C,*D);
902 MVT::MvNorm(*B,normsB2);
903 MVT::MvNorm(*C,normsC2);
904 MVT::MvNorm(*D,normsD1);
905 // check that input args are not modified
906 for (int i=0; i<p; i++) {
907 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
908 om->stream(Warnings)
909 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
910 << "Input arguments were modified." << endl;
911 return false;
912 }
913 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
914 om->stream(Warnings)
915 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
916 << "Input arguments were modified." << endl;
917 return false;
918 }
919 }
920 // next, try zero D
921 MVT::MvInit(*D);
922 MVT::MvAddMv(alpha,*B,beta,*C,*D);
923 MVT::MvNorm(*B,normsB2);
924 MVT::MvNorm(*C,normsC2);
925 MVT::MvNorm(*D,normsD2);
926 // check that input args are not modified and that D is the same
927 // as the above test
928 for (int i=0; i<p; i++) {
929 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
930 om->stream(Warnings)
931 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
932 << "Input arguments were modified." << endl;
933 return false;
934 }
935 else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
936 om->stream(Warnings)
937 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
938 << "Input arguments were modified." << endl;
939 return false;
940 }
941 else if ( SCT::magnitude( normsD1[i] - normsD2[i] ) > tol ) {
942 om->stream(Warnings)
943 << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
944 << "Results varies depending on initial state of dest vectors." << endl;
945 return false;
946 }
947 }
948 }
949
950 /*********** MvAddMv() ***********************************************
951 Similar to above, but where B or C are potentially the same
952 object as D. This case is commonly used, for example, to affect
953 A <- alpha*A
954 via
955 MvAddMv(alpha,A,zero,A,A)
956 ** OR **
957 MvAddMv(zero,A,alpha,A,A)
958
959 The result is that the operation has to be "atomic". That is,
960 B and C are no longer reliable after D is modified, so that
961 the assignment to D must be the last thing to occur.
962
963 D = alpha*B + beta*C
964
965 1) Use alpha==0,beta==1 and check that D == C
966 2) Use alpha==1,beta==0 and check that D == B
967 *********************************************************************/
968 {
969 const int p = 7;
970 Teuchos::RCP<MV> B, D;
971 Teuchos::RCP<const MV> C;
972 std::vector<MagType> normsB(p),
973 normsD(p);
974 std::vector<int> lclindex(p);
975 for (int i=0; i<p; i++) lclindex[i] = i;
976
977 B = MVT::Clone(*A,p);
978 C = MVT::CloneView(*B,lclindex);
979 D = MVT::CloneViewNonConst(*B,lclindex);
980
981 MVT::MvRandom(*B);
982 MVT::MvNorm(*B,normsB);
983
984 // check that 0*B+1*C == C
985 MVT::MvAddMv(zero,*B,one,*C,*D);
986 MVT::MvNorm(*D,normsD);
987 for (int i=0; i<p; i++) {
988 if ( SCT::magnitude( normsB[i] - normsD[i] ) > tol ) {
989 om->stream(Warnings)
990 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
991 << "Assignment did not work." << endl;
992 return false;
993 }
994 }
995
996 // check that 1*B+0*C == B
997 MVT::MvAddMv(one,*B,zero,*C,*D);
998 MVT::MvNorm(*D,normsD);
999 for (int i=0; i<p; i++) {
1000 if ( SCT::magnitude( normsB[i] - normsD[i] ) > tol ) {
1001 om->stream(Warnings)
1002 << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1003 << "Assignment did not work." << endl;
1004 return false;
1005 }
1006 }
1007
1008 }
1009
1010
1011 /*********** MvTimesMatAddMv() 7 by 5 ********************************
1012 C = alpha*B*SDM + beta*C
1013 1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
1014 2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
1015 3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
1016 4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
1017 5) Test with non-square matrices
1018 6) Always check that input arguments are not modified
1019 *********************************************************************/
1020 {
1021 const int p = 7, q = 5;
1022 Teuchos::RCP<MV> B, C;
1023 Teuchos::RCP<MV> Vp, Vq;
1024
1025 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1026 std::vector<MagType> normsC1(q), normsC2(q),
1027 normsB1(p), normsB2(p);
1028
1029 B = MVT::Clone(*A,p);
1030 C = MVT::Clone(*A,q);
1031
1032 // Create random SDM that is synchronized across processors.
1033 Vp = MVT::Clone(*A,p);
1034 Vq = MVT::Clone(*A,q);
1035 MVT::MvRandom(*Vp);
1036 MVT::MvRandom(*Vq);
1037 MVT::MvTransMv( one, *Vp, *Vq, SDM );
1038
1039 // Test 1: alpha==0, SDM!=0, beta==1 and check that C is unchanged
1040 MVT::MvRandom(*B);
1041 MVT::MvRandom(*C);
1042 MVT::MvNorm(*B,normsB1);
1043 MVT::MvNorm(*C,normsC1);
1044 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1045 MVT::MvNorm(*B,normsB2);
1046 MVT::MvNorm(*C,normsC2);
1047 for (int i=0; i<p; i++) {
1048 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1049 om->stream(Warnings)
1050 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1051 << "Input vectors were modified." << endl;
1052 return false;
1053 }
1054 }
1055 for (int i=0; i<q; i++) {
1056 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1057 om->stream(Warnings)
1058 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1059 << "Arithmetic test 1 failed." << endl;
1060 return false;
1061 }
1062 }
1063
1064 // Test 2: alpha==0, SDM!=0, beta==0 and check that C is set to zero
1065 MVT::MvRandom(*B);
1066 MVT::MvRandom(*C);
1067 MVT::MvNorm(*B,normsB1);
1068 MVT::MvNorm(*C,normsC1);
1069 MVT::MvRandom(*Vp);
1070 MVT::MvRandom(*Vq);
1071 MVT::MvTransMv( one, *Vp, *Vq, SDM );
1072 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1073 MVT::MvNorm(*B,normsB2);
1074 MVT::MvNorm(*C,normsC2);
1075 for (int i=0; i<p; i++) {
1076 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1077 om->stream(Warnings)
1078 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1079 << "Input vectors were modified." << endl;
1080 return false;
1081 }
1082 }
1083 for (int i=0; i<q; i++) {
1084 if ( normsC2[i] != zero ) {
1085 om->stream(Warnings)
1086 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1087 << "Arithmetic test 2 failed: "
1088 << normsC2[i]
1089 << " != "
1090 << zero
1091 << endl;
1092 return false;
1093 }
1094 }
1095
1096 // Test 3: alpha==1, SDM==|I|, beta==0 and check that C is set to B
1097 // |0|
1098 MVT::MvRandom(*B);
1099 MVT::MvRandom(*C);
1100 MVT::MvNorm(*B,normsB1);
1101 MVT::MvNorm(*C,normsC1);
1102 SDM.scale(zero);
1103 for (int i=0; i<q; i++) {
1104 SDM(i,i) = one;
1105 }
1106 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1107 MVT::MvNorm(*B,normsB2);
1108 MVT::MvNorm(*C,normsC2);
1109 for (int i=0; i<p; i++) {
1110 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1111 om->stream(Warnings)
1112 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1113 << "Input vectors were modified." << endl;
1114 return false;
1115 }
1116 }
1117 for (int i=0; i<q; i++) {
1118 if ( SCT::magnitude( normsB1[i] - normsC2[i] ) > tol ) {
1119 om->stream(Warnings)
1120 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1121 << "Arithmetic test 3 failed: "
1122 << normsB1[i]
1123 << " != "
1124 << normsC2[i]
1125 << endl;
1126 return false;
1127 }
1128 }
1129
1130 // Test 4: alpha==1, SDM==0, beta==1 and check that C is unchanged
1131 MVT::MvRandom(*B);
1132 MVT::MvRandom(*C);
1133 MVT::MvNorm(*B,normsB1);
1134 MVT::MvNorm(*C,normsC1);
1135 SDM.scale(zero);
1136 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1137 MVT::MvNorm(*B,normsB2);
1138 MVT::MvNorm(*C,normsC2);
1139 for (int i=0; i<p; i++) {
1140 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1141 om->stream(Warnings)
1142 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1143 << "Input vectors were modified." << endl;
1144 return false;
1145 }
1146 }
1147 for (int i=0; i<q; i++) {
1148 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1149 om->stream(Warnings)
1150 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1151 << "Arithmetic test 4 failed." << endl;
1152 return false;
1153 }
1154 }
1155 }
1156
1157 /*********** MvTimesMatAddMv() 5 by 7 ********************************
1158 C = alpha*B*SDM + beta*C
1159 1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
1160 2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
1161 3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
1162 4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
1163 5) Test with non-square matrices
1164 6) Always check that input arguments are not modified
1165 *********************************************************************/
1166 {
1167 const int p = 5, q = 7;
1168 Teuchos::RCP<MV> B, C;
1169 Teuchos::RCP<MV> Vp, Vq;
1170 Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1171 std::vector<MagType> normsC1(q), normsC2(q),
1172 normsB1(p), normsB2(p);
1173
1174 B = MVT::Clone(*A,p);
1175 C = MVT::Clone(*A,q);
1176
1177 // Create random SDM that is synchronized across processors.
1178 Vp = MVT::Clone(*A,p);
1179 Vq = MVT::Clone(*A,q);
1180 MVT::MvRandom(*Vp);
1181 MVT::MvRandom(*Vq);
1182 MVT::MvTransMv( one, *Vp, *Vq, SDM );
1183
1184 // Test 5: alpha==0, SDM!=0, beta==1 and check that C is unchanged
1185 MVT::MvRandom(*B);
1186 MVT::MvRandom(*C);
1187 MVT::MvNorm(*B,normsB1);
1188 MVT::MvNorm(*C,normsC1);
1189 MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1190 MVT::MvNorm(*B,normsB2);
1191 MVT::MvNorm(*C,normsC2);
1192 for (int i=0; i<p; i++) {
1193 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1194 om->stream(Warnings)
1195 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1196 << "Input vectors were modified." << endl;
1197 return false;
1198 }
1199 }
1200 for (int i=0; i<q; i++) {
1201 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1202 om->stream(Warnings)
1203 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1204 << "Arithmetic test 5 failed." << endl;
1205 return false;
1206 }
1207 }
1208
1209 // Test 6: alpha==0, SDM!=0, beta==0 and check that C is set to zero
1210 MVT::MvRandom(*B);
1211 MVT::MvRandom(*C);
1212 MVT::MvNorm(*B,normsB1);
1213 MVT::MvNorm(*C,normsC1);
1214 MVT::MvRandom(*Vp);
1215 MVT::MvRandom(*Vq);
1216 MVT::MvTransMv( one, *Vp, *Vq, SDM );
1217 MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1218 MVT::MvNorm(*B,normsB2);
1219 MVT::MvNorm(*C,normsC2);
1220 for (int i=0; i<p; i++) {
1221 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1222 om->stream(Warnings)
1223 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1224 << "Input vectors were modified." << endl;
1225 return false;
1226 }
1227 }
1228 for (int i=0; i<q; i++) {
1229 if ( normsC2[i] != zero ) {
1230 om->stream(Warnings)
1231 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1232 << "Arithmetic test 6 failed: "
1233 << normsC2[i]
1234 << " != "
1235 << zero
1236 << endl;
1237 return false;
1238 }
1239 }
1240
1241 // Test 7: alpha==1, SDM==[I 0], beta==0 and check that C is set to B
1242 MVT::MvRandom(*B);
1243 MVT::MvRandom(*C);
1244 MVT::MvNorm(*B,normsB1);
1245 MVT::MvNorm(*C,normsC1);
1246 SDM.scale(zero);
1247 for (int i=0; i<p; i++) {
1248 SDM(i,i) = one;
1249 }
1250 MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1251 MVT::MvNorm(*B,normsB2);
1252 MVT::MvNorm(*C,normsC2);
1253 for (int i=0; i<p; i++) {
1254 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1255 om->stream(Warnings)
1256 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1257 << "Input vectors were modified." << endl;
1258 return false;
1259 }
1260 }
1261 for (int i=0; i<p; i++) {
1262 if ( SCT::magnitude( normsB1[i] - normsC2[i] ) > tol ) {
1263 om->stream(Warnings)
1264 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1265 << "Arithmetic test 7 failed." << endl;
1266 return false;
1267 }
1268 }
1269 for (int i=p; i<q; i++) {
1270 if ( normsC2[i] != zero ) {
1271 om->stream(Warnings)
1272 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1273 << "Arithmetic test 7 failed." << endl;
1274 return false;
1275 }
1276 }
1277
1278 // Test 8: alpha==1, SDM==0, beta==1 and check that C is unchanged
1279 MVT::MvRandom(*B);
1280 MVT::MvRandom(*C);
1281 MVT::MvNorm(*B,normsB1);
1282 MVT::MvNorm(*C,normsC1);
1283 SDM.scale(zero);
1284 MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1285 MVT::MvNorm(*B,normsB2);
1286 MVT::MvNorm(*C,normsC2);
1287 for (int i=0; i<p; i++) {
1288 if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1289 om->stream(Warnings)
1290 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1291 << "Input vectors were modified." << endl;
1292 return false;
1293 }
1294 }
1295 for (int i=0; i<q; i++) {
1296 if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1297 om->stream(Warnings)
1298 << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1299 << "Arithmetic test 8 failed." << endl;
1300 return false;
1301 }
1302 }
1303 }
1304
1305 return true;
1306
1307 }
1308
1309
1310
1316 template< class ScalarType, class MV, class OP>
1317 bool TestOperatorTraits(
1318 const Teuchos::RCP<OutputManager<ScalarType> > &om,
1319 const Teuchos::RCP<const MV> &A,
1320 const Teuchos::RCP<const OP> &M) {
1321
1322 using std::endl;
1323
1324 /* OPT Contract:
1325 Apply()
1326 MV: OP*zero == zero
1327 Warn if OP is not deterministic (OP*A != OP*A)
1328 Does not modify input arguments
1329 *********************************************************************/
1330
1331 typedef MultiVecTraits<ScalarType, MV> MVT;
1332 typedef Teuchos::ScalarTraits<ScalarType> SCT;
1333 typedef OperatorTraits<ScalarType, MV, OP> OPT;
1334 typedef typename SCT::magnitudeType MagType;
1335
1336 const MagType tol = SCT::eps()*100;
1337
1338 const int numvecs = 10;
1339
1340 Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs),
1341 C = MVT::Clone(*A,numvecs);
1342
1343 std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
1344 normsC1(numvecs), normsC2(numvecs);
1345
1346 /*********** Apply() *************************************************
1347 Verify:
1348 1) OP*B == OP*B; OP is deterministic (just warn on this)
1349 2) OP*zero == 0
1350 3) OP*B doesn't modify B
1351 4) OP*B is invariant under initial state of destination vectors
1352 *********************************************************************/
1353 MVT::MvInit(*B);
1354 MVT::MvRandom(*C);
1355 MVT::MvNorm(*B,normsB1);
1356 OPT::Apply(*M,*B,*C);
1357 MVT::MvNorm(*B,normsB2);
1358 MVT::MvNorm(*C,normsC2);
1359 for (int i=0; i<numvecs; i++) {
1360 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1361 om->stream(Warnings)
1362 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
1363 << "Apply() modified the input vectors." << endl;
1364 return false;
1365 }
1366 if (normsC2[i] != SCT::zero()) {
1367 om->stream(Warnings)
1368 << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
1369 << "Operator applied to zero did not return zero." << endl;
1370 return false;
1371 }
1372 }
1373
1374 // If we send in a random matrix, we should not get a zero return
1375 MVT::MvRandom(*B);
1376 MVT::MvNorm(*B,normsB1);
1377 OPT::Apply(*M,*B,*C);
1378 MVT::MvNorm(*B,normsB2);
1379 MVT::MvNorm(*C,normsC2);
1380 bool ZeroWarning = false;
1381 for (int i=0; i<numvecs; i++) {
1382 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1383 om->stream(Warnings)
1384 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
1385 << "Apply() modified the input vectors." << endl;
1386 return false;
1387 }
1388 if (normsC2[i] == SCT::zero() && ZeroWarning==false ) {
1389 om->stream(Warnings)
1390 << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
1391 << "Operator applied to random vectors returned zero." << endl;
1392 ZeroWarning = true;
1393 }
1394 }
1395
1396 // Apply operator with C init'd to zero
1397 MVT::MvRandom(*B);
1398 MVT::MvNorm(*B,normsB1);
1399 MVT::MvInit(*C);
1400 OPT::Apply(*M,*B,*C);
1401 MVT::MvNorm(*B,normsB2);
1402 MVT::MvNorm(*C,normsC1);
1403 for (int i=0; i<numvecs; i++) {
1404 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1405 om->stream(Warnings)
1406 << "*** ERROR *** OperatorTraits::Apply() [3]" << endl
1407 << "Apply() modified the input vectors." << endl;
1408 return false;
1409 }
1410 }
1411
1412 // Apply operator with C init'd to random
1413 // Check that result is the same as before; warn if not.
1414 // This could be a result of a bug, or a stochastic
1415 // operator. We do not want to prejudice against a
1416 // stochastic operator.
1417 MVT::MvRandom(*C);
1418 OPT::Apply(*M,*B,*C);
1419 MVT::MvNorm(*B,normsB2);
1420 MVT::MvNorm(*C,normsC2);
1421 for (int i=0; i<numvecs; i++) {
1422 if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1423 om->stream(Warnings)
1424 << "*** ERROR *** OperatorTraits::Apply() [4]" << endl
1425 << "Apply() modified the input vectors." << endl;
1426 return false;
1427 }
1428 if ( SCT::magnitude( normsC1[i] - normsC2[i]) > tol ) {
1429 om->stream(Warnings)
1430 << endl
1431 << "*** WARNING *** OperatorTraits::Apply() [4]" << endl
1432 << "Apply() returned two different results." << endl << endl;
1433 }
1434 }
1435
1436 return true;
1437
1438 }
1439
1440}
1441
1442#endif
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
AnasaziMultiVecTraits.hpp
Declaration of basic traits for the multivector type.
AnasaziOperatorTraits.hpp
Virtual base class which defines basic traits for the operator type.
AnasaziOutputManager.hpp
Abstract class definition for Anasazi Output Managers.
AnasaziTypes.hpp
Types and exceptions used within Anasazi solvers and interfaces.
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition AnasaziMultiVecTraits.hpp:95
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition AnasaziOperatorTraits.hpp:53
Anasazi::OutputManager
Output managers remove the need for the eigensolver to know any information about the required output...
Definition AnasaziOutputManager.hpp:36
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Anasazi::TestOperatorTraits
bool TestOperatorTraits(const Teuchos::RCP< OutputManager< ScalarType > > &om, const Teuchos::RCP< const MV > &A, const Teuchos::RCP< const OP > &M)
This function tests the correctness of an operator implementation with respect to an OperatorTraits s...
Definition AnasaziMVOPTester.hpp:1317
Anasazi::TestMultiVecTraits
bool TestMultiVecTraits(const Teuchos::RCP< OutputManager< ScalarType > > &om, const Teuchos::RCP< const MV > &A)
This is a function to test the correctness of a MultiVecTraits specialization and multivector impleme...
Definition AnasaziMVOPTester.hpp:46
Generated by doxygen 1.9.8