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