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