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