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