Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AnasaziMVOPTester.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
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_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::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 = 0.5 * SCT::one(),
866  beta = 0.33 * SCT::one();
867 
868  B = MVT::Clone(*A,p);
869  C = MVT::Clone(*A,p);
870  D = MVT::Clone(*A,p);
871 
872  MVT::MvRandom(*B);
873  MVT::MvRandom(*C);
874  MVT::MvNorm(*B,normsB1);
875  MVT::MvNorm(*C,normsC1);
876 
877  // check that 0*B+1*C == C
878  MVT::MvAddMv(zero,*B,one,*C,*D);
879  MVT::MvNorm(*B,normsB2);
880  MVT::MvNorm(*C,normsC2);
881  MVT::MvNorm(*D,normsD1);
882  for (int i=0; i<p; i++) {
883  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
884  om->stream(Warnings)
885  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
886  << "Input arguments were modified." << endl;
887  return false;
888  }
889  else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
890  om->stream(Warnings)
891  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
892  << "Input arguments were modified." << endl;
893  return false;
894  }
895  else if ( SCT::magnitude(normsC1[i]-normsD1[i]) > tol ) {
896  om->stream(Warnings)
897  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
898  << "Assignment did not work." << endl;
899  return false;
900  }
901  }
902 
903  // check that 1*B+0*C == B
904  MVT::MvAddMv(one,*B,zero,*C,*D);
905  MVT::MvNorm(*B,normsB2);
906  MVT::MvNorm(*C,normsC2);
907  MVT::MvNorm(*D,normsD1);
908  for (int i=0; i<p; i++) {
909  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
910  om->stream(Warnings)
911  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
912  << "Input arguments were modified." << endl;
913  return false;
914  }
915  else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
916  om->stream(Warnings)
917  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
918  << "Input arguments were modified." << endl;
919  return false;
920  }
921  else if ( SCT::magnitude( normsB1[i] - normsD1[i] ) > tol ) {
922  om->stream(Warnings)
923  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
924  << "Assignment did not work." << endl;
925  return false;
926  }
927  }
928 
929  // check that alpha*B+beta*C -> D is invariant under initial D
930  // first, try random D
931  MVT::MvRandom(*D);
932  MVT::MvAddMv(alpha,*B,beta,*C,*D);
933  MVT::MvNorm(*B,normsB2);
934  MVT::MvNorm(*C,normsC2);
935  MVT::MvNorm(*D,normsD1);
936  // check that input args are not modified
937  for (int i=0; i<p; i++) {
938  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
939  om->stream(Warnings)
940  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
941  << "Input arguments were modified." << endl;
942  return false;
943  }
944  else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
945  om->stream(Warnings)
946  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
947  << "Input arguments were modified." << endl;
948  return false;
949  }
950  }
951  // next, try zero D
952  MVT::MvInit(*D);
953  MVT::MvAddMv(alpha,*B,beta,*C,*D);
954  MVT::MvNorm(*B,normsB2);
955  MVT::MvNorm(*C,normsC2);
956  MVT::MvNorm(*D,normsD2);
957  // check that input args are not modified and that D is the same
958  // as the above test
959  for (int i=0; i<p; i++) {
960  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
961  om->stream(Warnings)
962  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
963  << "Input arguments were modified." << endl;
964  return false;
965  }
966  else if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
967  om->stream(Warnings)
968  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
969  << "Input arguments were modified." << endl;
970  return false;
971  }
972  else if ( SCT::magnitude( normsD1[i] - normsD2[i] ) > tol ) {
973  om->stream(Warnings)
974  << "*** ERROR *** MultiVecTraits::MvAddMv()." << endl
975  << "Results varies depending on initial state of dest vectors." << endl;
976  return false;
977  }
978  }
979  }
980 
981  /*********** MvAddMv() ***********************************************
982  Similar to above, but where B or C are potentially the same
983  object as D. This case is commonly used, for example, to affect
984  A <- alpha*A
985  via
986  MvAddMv(alpha,A,zero,A,A)
987  ** OR **
988  MvAddMv(zero,A,alpha,A,A)
989 
990  The result is that the operation has to be "atomic". That is,
991  B and C are no longer reliable after D is modified, so that
992  the assignment to D must be the last thing to occur.
993 
994  D = alpha*B + beta*C
995 
996  1) Use alpha==0,beta==1 and check that D == C
997  2) Use alpha==1,beta==0 and check that D == B
998  *********************************************************************/
999  {
1000  const int p = 7;
1001  Teuchos::RCP<MV> B, D;
1002  Teuchos::RCP<const MV> C;
1003  std::vector<MagType> normsB(p),
1004  normsD(p);
1005  std::vector<int> lclindex(p);
1006  for (int i=0; i<p; i++) lclindex[i] = i;
1007 
1008  B = MVT::Clone(*A,p);
1009  C = MVT::CloneView(*B,lclindex);
1010  D = MVT::CloneViewNonConst(*B,lclindex);
1011 
1012  MVT::MvRandom(*B);
1013  MVT::MvNorm(*B,normsB);
1014 
1015  // check that 0*B+1*C == C
1016  MVT::MvAddMv(zero,*B,one,*C,*D);
1017  MVT::MvNorm(*D,normsD);
1018  for (int i=0; i<p; i++) {
1019  if ( SCT::magnitude( normsB[i] - normsD[i] ) > tol ) {
1020  om->stream(Warnings)
1021  << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1022  << "Assignment did not work." << endl;
1023  return false;
1024  }
1025  }
1026 
1027  // check that 1*B+0*C == B
1028  MVT::MvAddMv(one,*B,zero,*C,*D);
1029  MVT::MvNorm(*D,normsD);
1030  for (int i=0; i<p; i++) {
1031  if ( SCT::magnitude( normsB[i] - normsD[i] ) > tol ) {
1032  om->stream(Warnings)
1033  << "*** ERROR *** MultiVecTraits::MvAddMv() #2" << endl
1034  << "Assignment did not work." << endl;
1035  return false;
1036  }
1037  }
1038 
1039  }
1040 
1041 
1042  /*********** MvTimesMatAddMv() 7 by 5 ********************************
1043  C = alpha*B*SDM + beta*C
1044  1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
1045  2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
1046  3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
1047  4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
1048  5) Test with non-square matrices
1049  6) Always check that input arguments are not modified
1050  *********************************************************************/
1051  {
1052  const int p = 7, q = 5;
1053  Teuchos::RCP<MV> B, C;
1054  Teuchos::RCP<MV> Vp, Vq;
1055 
1056  Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1057  std::vector<MagType> normsC1(q), normsC2(q),
1058  normsB1(p), normsB2(p);
1059 
1060  B = MVT::Clone(*A,p);
1061  C = MVT::Clone(*A,q);
1062 
1063  // Create random SDM that is synchronized across processors.
1064  Vp = MVT::Clone(*A,p);
1065  Vq = MVT::Clone(*A,q);
1066  MVT::MvRandom(*Vp);
1067  MVT::MvRandom(*Vq);
1068  MVT::MvTransMv( one, *Vp, *Vq, SDM );
1069 
1070  // Test 1: alpha==0, SDM!=0, beta==1 and check that C is unchanged
1071  MVT::MvRandom(*B);
1072  MVT::MvRandom(*C);
1073  MVT::MvNorm(*B,normsB1);
1074  MVT::MvNorm(*C,normsC1);
1075  MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1076  MVT::MvNorm(*B,normsB2);
1077  MVT::MvNorm(*C,normsC2);
1078  for (int i=0; i<p; i++) {
1079  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1080  om->stream(Warnings)
1081  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1082  << "Input vectors were modified." << endl;
1083  return false;
1084  }
1085  }
1086  for (int i=0; i<q; i++) {
1087  if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1088  om->stream(Warnings)
1089  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1090  << "Arithmetic test 1 failed." << endl;
1091  return false;
1092  }
1093  }
1094 
1095  // Test 2: alpha==0, SDM!=0, beta==0 and check that C is set to zero
1096  MVT::MvRandom(*B);
1097  MVT::MvRandom(*C);
1098  MVT::MvNorm(*B,normsB1);
1099  MVT::MvNorm(*C,normsC1);
1100  MVT::MvRandom(*Vp);
1101  MVT::MvRandom(*Vq);
1102  MVT::MvTransMv( one, *Vp, *Vq, SDM );
1103  MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1104  MVT::MvNorm(*B,normsB2);
1105  MVT::MvNorm(*C,normsC2);
1106  for (int i=0; i<p; i++) {
1107  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1108  om->stream(Warnings)
1109  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1110  << "Input vectors were modified." << endl;
1111  return false;
1112  }
1113  }
1114  for (int i=0; i<q; i++) {
1115  if ( normsC2[i] != zero ) {
1116  om->stream(Warnings)
1117  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1118  << "Arithmetic test 2 failed: "
1119  << normsC2[i]
1120  << " != "
1121  << zero
1122  << endl;
1123  return false;
1124  }
1125  }
1126 
1127  // Test 3: alpha==1, SDM==|I|, beta==0 and check that C is set to B
1128  // |0|
1129  MVT::MvRandom(*B);
1130  MVT::MvRandom(*C);
1131  MVT::MvNorm(*B,normsB1);
1132  MVT::MvNorm(*C,normsC1);
1133  SDM.scale(zero);
1134  for (int i=0; i<q; i++) {
1135  SDM(i,i) = one;
1136  }
1137  MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1138  MVT::MvNorm(*B,normsB2);
1139  MVT::MvNorm(*C,normsC2);
1140  for (int i=0; i<p; i++) {
1141  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1142  om->stream(Warnings)
1143  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1144  << "Input vectors were modified." << endl;
1145  return false;
1146  }
1147  }
1148  for (int i=0; i<q; i++) {
1149  if ( SCT::magnitude( normsB1[i] - normsC2[i] ) > tol ) {
1150  om->stream(Warnings)
1151  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1152  << "Arithmetic test 3 failed: "
1153  << normsB1[i]
1154  << " != "
1155  << normsC2[i]
1156  << endl;
1157  return false;
1158  }
1159  }
1160 
1161  // Test 4: alpha==1, SDM==0, beta==1 and check that C is unchanged
1162  MVT::MvRandom(*B);
1163  MVT::MvRandom(*C);
1164  MVT::MvNorm(*B,normsB1);
1165  MVT::MvNorm(*C,normsC1);
1166  SDM.scale(zero);
1167  MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1168  MVT::MvNorm(*B,normsB2);
1169  MVT::MvNorm(*C,normsC2);
1170  for (int i=0; i<p; i++) {
1171  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1172  om->stream(Warnings)
1173  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1174  << "Input vectors were modified." << endl;
1175  return false;
1176  }
1177  }
1178  for (int i=0; i<q; i++) {
1179  if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1180  om->stream(Warnings)
1181  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1182  << "Arithmetic test 4 failed." << endl;
1183  return false;
1184  }
1185  }
1186  }
1187 
1188  /*********** MvTimesMatAddMv() 5 by 7 ********************************
1189  C = alpha*B*SDM + beta*C
1190  1) Use alpha==0, SDM!=0, beta==1 and check that C is unchanged
1191  2) Use alpha==0, SDM!=0, beta==0 and check that C is set to zero
1192  3) Use alpha==1, SDM==I, beta==0 and check that C is set to B
1193  4) Use alpha==1, SDM==0, beta==1 and check that C is unchanged
1194  5) Test with non-square matrices
1195  6) Always check that input arguments are not modified
1196  *********************************************************************/
1197  {
1198  const int p = 5, q = 7;
1199  Teuchos::RCP<MV> B, C;
1200  Teuchos::RCP<MV> Vp, Vq;
1201  Teuchos::SerialDenseMatrix<int,ScalarType> SDM(p,q);
1202  std::vector<MagType> normsC1(q), normsC2(q),
1203  normsB1(p), normsB2(p);
1204 
1205  B = MVT::Clone(*A,p);
1206  C = MVT::Clone(*A,q);
1207 
1208  // Create random SDM that is synchronized across processors.
1209  Vp = MVT::Clone(*A,p);
1210  Vq = MVT::Clone(*A,q);
1211  MVT::MvRandom(*Vp);
1212  MVT::MvRandom(*Vq);
1213  MVT::MvTransMv( one, *Vp, *Vq, SDM );
1214 
1215  // Test 5: alpha==0, SDM!=0, beta==1 and check that C is unchanged
1216  MVT::MvRandom(*B);
1217  MVT::MvRandom(*C);
1218  MVT::MvNorm(*B,normsB1);
1219  MVT::MvNorm(*C,normsC1);
1220  MVT::MvTimesMatAddMv(zero,*B,SDM,one,*C);
1221  MVT::MvNorm(*B,normsB2);
1222  MVT::MvNorm(*C,normsC2);
1223  for (int i=0; i<p; i++) {
1224  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1225  om->stream(Warnings)
1226  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1227  << "Input vectors were modified." << endl;
1228  return false;
1229  }
1230  }
1231  for (int i=0; i<q; i++) {
1232  if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1233  om->stream(Warnings)
1234  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1235  << "Arithmetic test 5 failed." << endl;
1236  return false;
1237  }
1238  }
1239 
1240  // Test 6: alpha==0, SDM!=0, beta==0 and check that C is set to zero
1241  MVT::MvRandom(*B);
1242  MVT::MvRandom(*C);
1243  MVT::MvNorm(*B,normsB1);
1244  MVT::MvNorm(*C,normsC1);
1245  MVT::MvRandom(*Vp);
1246  MVT::MvRandom(*Vq);
1247  MVT::MvTransMv( one, *Vp, *Vq, SDM );
1248  MVT::MvTimesMatAddMv(zero,*B,SDM,zero,*C);
1249  MVT::MvNorm(*B,normsB2);
1250  MVT::MvNorm(*C,normsC2);
1251  for (int i=0; i<p; i++) {
1252  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1253  om->stream(Warnings)
1254  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1255  << "Input vectors were modified." << endl;
1256  return false;
1257  }
1258  }
1259  for (int i=0; i<q; i++) {
1260  if ( normsC2[i] != zero ) {
1261  om->stream(Warnings)
1262  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1263  << "Arithmetic test 6 failed: "
1264  << normsC2[i]
1265  << " != "
1266  << zero
1267  << endl;
1268  return false;
1269  }
1270  }
1271 
1272  // Test 7: alpha==1, SDM==[I 0], beta==0 and check that C is set to B
1273  MVT::MvRandom(*B);
1274  MVT::MvRandom(*C);
1275  MVT::MvNorm(*B,normsB1);
1276  MVT::MvNorm(*C,normsC1);
1277  SDM.scale(zero);
1278  for (int i=0; i<p; i++) {
1279  SDM(i,i) = one;
1280  }
1281  MVT::MvTimesMatAddMv(one,*B,SDM,zero,*C);
1282  MVT::MvNorm(*B,normsB2);
1283  MVT::MvNorm(*C,normsC2);
1284  for (int i=0; i<p; i++) {
1285  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1286  om->stream(Warnings)
1287  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1288  << "Input vectors were modified." << endl;
1289  return false;
1290  }
1291  }
1292  for (int i=0; i<p; i++) {
1293  if ( SCT::magnitude( normsB1[i] - normsC2[i] ) > tol ) {
1294  om->stream(Warnings)
1295  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1296  << "Arithmetic test 7 failed." << endl;
1297  return false;
1298  }
1299  }
1300  for (int i=p; i<q; i++) {
1301  if ( normsC2[i] != zero ) {
1302  om->stream(Warnings)
1303  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1304  << "Arithmetic test 7 failed." << endl;
1305  return false;
1306  }
1307  }
1308 
1309  // Test 8: alpha==1, SDM==0, beta==1 and check that C is unchanged
1310  MVT::MvRandom(*B);
1311  MVT::MvRandom(*C);
1312  MVT::MvNorm(*B,normsB1);
1313  MVT::MvNorm(*C,normsC1);
1314  SDM.scale(zero);
1315  MVT::MvTimesMatAddMv(one,*B,SDM,one,*C);
1316  MVT::MvNorm(*B,normsB2);
1317  MVT::MvNorm(*C,normsC2);
1318  for (int i=0; i<p; i++) {
1319  if ( SCT::magnitude( normsB1[i] - normsB2[i] ) > tol ) {
1320  om->stream(Warnings)
1321  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1322  << "Input vectors were modified." << endl;
1323  return false;
1324  }
1325  }
1326  for (int i=0; i<q; i++) {
1327  if ( SCT::magnitude( normsC1[i] - normsC2[i] ) > tol ) {
1328  om->stream(Warnings)
1329  << "*** ERROR *** MultiVecTraits::MvTimesMatAddMv()." << endl
1330  << "Arithmetic test 8 failed." << endl;
1331  return false;
1332  }
1333  }
1334  }
1335 
1336  return true;
1337 
1338  }
1339 
1340 
1341 
1347  template< class ScalarType, class MV, class OP>
1348  bool TestOperatorTraits(
1349  const Teuchos::RCP<OutputManager<ScalarType> > &om,
1350  const Teuchos::RCP<const MV> &A,
1351  const Teuchos::RCP<const OP> &M) {
1352 
1353  using std::endl;
1354 
1355  /* OPT Contract:
1356  Apply()
1357  MV: OP*zero == zero
1358  Warn if OP is not deterministic (OP*A != OP*A)
1359  Does not modify input arguments
1360  *********************************************************************/
1361 
1362  typedef MultiVecTraits<ScalarType, MV> MVT;
1363  typedef Teuchos::ScalarTraits<ScalarType> SCT;
1364  typedef OperatorTraits<ScalarType, MV, OP> OPT;
1365  typedef typename SCT::magnitudeType MagType;
1366 
1367  const MagType tol = SCT::eps()*100;
1368 
1369  const int numvecs = 10;
1370 
1371  Teuchos::RCP<MV> B = MVT::Clone(*A,numvecs),
1372  C = MVT::Clone(*A,numvecs);
1373 
1374  std::vector<MagType> normsB1(numvecs), normsB2(numvecs),
1375  normsC1(numvecs), normsC2(numvecs);
1376 
1377  /*********** Apply() *************************************************
1378  Verify:
1379  1) OP*B == OP*B; OP is deterministic (just warn on this)
1380  2) OP*zero == 0
1381  3) OP*B doesn't modify B
1382  4) OP*B is invariant under initial state of destination vectors
1383  *********************************************************************/
1384  MVT::MvInit(*B);
1385  MVT::MvRandom(*C);
1386  MVT::MvNorm(*B,normsB1);
1387  OPT::Apply(*M,*B,*C);
1388  MVT::MvNorm(*B,normsB2);
1389  MVT::MvNorm(*C,normsC2);
1390  for (int i=0; i<numvecs; i++) {
1391  if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1392  om->stream(Warnings)
1393  << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
1394  << "Apply() modified the input vectors." << endl;
1395  return false;
1396  }
1397  if (normsC2[i] != SCT::zero()) {
1398  om->stream(Warnings)
1399  << "*** ERROR *** OperatorTraits::Apply() [1]" << endl
1400  << "Operator applied to zero did not return zero." << endl;
1401  return false;
1402  }
1403  }
1404 
1405  // If we send in a random matrix, we should not get a zero return
1406  MVT::MvRandom(*B);
1407  MVT::MvNorm(*B,normsB1);
1408  OPT::Apply(*M,*B,*C);
1409  MVT::MvNorm(*B,normsB2);
1410  MVT::MvNorm(*C,normsC2);
1411  bool ZeroWarning = false;
1412  for (int i=0; i<numvecs; i++) {
1413  if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1414  om->stream(Warnings)
1415  << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
1416  << "Apply() modified the input vectors." << endl;
1417  return false;
1418  }
1419  if (normsC2[i] == SCT::zero() && ZeroWarning==false ) {
1420  om->stream(Warnings)
1421  << "*** ERROR *** OperatorTraits::Apply() [2]" << endl
1422  << "Operator applied to random vectors returned zero." << endl;
1423  ZeroWarning = true;
1424  }
1425  }
1426 
1427  // Apply operator with C init'd to zero
1428  MVT::MvRandom(*B);
1429  MVT::MvNorm(*B,normsB1);
1430  MVT::MvInit(*C);
1431  OPT::Apply(*M,*B,*C);
1432  MVT::MvNorm(*B,normsB2);
1433  MVT::MvNorm(*C,normsC1);
1434  for (int i=0; i<numvecs; i++) {
1435  if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1436  om->stream(Warnings)
1437  << "*** ERROR *** OperatorTraits::Apply() [3]" << endl
1438  << "Apply() modified the input vectors." << endl;
1439  return false;
1440  }
1441  }
1442 
1443  // Apply operator with C init'd to random
1444  // Check that result is the same as before; warn if not.
1445  // This could be a result of a bug, or a stochastic
1446  // operator. We do not want to prejudice against a
1447  // stochastic operator.
1448  MVT::MvRandom(*C);
1449  OPT::Apply(*M,*B,*C);
1450  MVT::MvNorm(*B,normsB2);
1451  MVT::MvNorm(*C,normsC2);
1452  for (int i=0; i<numvecs; i++) {
1453  if ( SCT::magnitude( normsB2[i] - normsB1[i] ) > tol ) {
1454  om->stream(Warnings)
1455  << "*** ERROR *** OperatorTraits::Apply() [4]" << endl
1456  << "Apply() modified the input vectors." << endl;
1457  return false;
1458  }
1459  if ( SCT::magnitude( normsC1[i] - normsC2[i]) > tol ) {
1460  om->stream(Warnings)
1461  << endl
1462  << "*** WARNING *** OperatorTraits::Apply() [4]" << endl
1463  << "Apply() returned two different results." << endl << endl;
1464  }
1465  }
1466 
1467  return true;
1468 
1469  }
1470 
1471 }
1472 
1473 #endif
Teuchos::SerialDenseMatrix::normOne
ScalarTraits< ScalarType >::magnitudeType normOne() const
Anasazi::TestOperatorTraits
bool TestOperatorTraits(const Teuchos::RCP< OutputManager< ScalarType > > &om, const Teuchos::RCP< const MV > &A, const Teuchos::RCP< const OP > &M)
This function tests the correctness of an operator implementation with respect to an OperatorTraits s...
Definition: AnasaziMVOPTester.hpp:1348
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
Teuchos::SerialDenseMatrix::scale
int scale(const ScalarType alpha)
AnasaziOutputManager.hpp
Abstract class definition for Anasazi Output Managers.
Anasazi::OutputManager
Output managers remove the need for the eigensolver to know any information about the required output...
Definition: AnasaziOutputManager.hpp:68
Teuchos::SerialDenseMatrix::putScalar
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: AnasaziMultiVecTraits.hpp: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...
Teuchos::SerialDenseMatrix::numCols
OrdinalType numCols() const
AnasaziTypes.hpp
Types and exceptions used within Anasazi solvers and interfaces.
Teuchos_RCP.hpp
Teuchos_as.hpp
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::SerialDenseMatrix::numRows
OrdinalType numRows() const
Generated by   doxygen 1.8.5