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