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