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