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