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