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