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