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