Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosOrthoManagerTest.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
13 
14 #include <BelosConfigDefs.hpp>
15 #include <BelosMultiVecTraits.hpp>
16 #include <BelosOutputManager.hpp>
17 #include <BelosOrthoManagerFactory.hpp>
18 #include <Teuchos_StandardCatchMacros.hpp>
19 #include <Teuchos_TimeMonitor.hpp>
20 #include <Teuchos_SerialDenseHelpers.hpp>
21 #include <iostream>
22 #include <stdexcept>
23 
24 using std::endl;
25 
26 namespace Belos {
27  namespace Test {
28 
33  template<class Scalar, class MV>
34  class OrthoManagerBenchmarker {
35  private:
36  typedef Scalar scalar_type;
37  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
38  typedef MultiVecTraits<Scalar, MV> MVT;
39  typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
40 
41  public:
51  static void
52  baseline (const Teuchos::RCP<const MV>& X,
53  const int numCols,
54  const int numBlocks,
55  const int numTrials)
56  {
57  using Teuchos::Array;
58  using Teuchos::RCP;
59  using Teuchos::rcp;
60  using Teuchos::Time;
61  using Teuchos::TimeMonitor;
62 
63  // Make some blocks to "orthogonalize." Fill with random
64  // data. We only need X so that we can make clones (it knows
65  // its data distribution).
66  Array<RCP<MV> > V (numBlocks);
67  for (int k = 0; k < numBlocks; ++k) {
68  V[k] = MVT::Clone (*X, numCols);
69  MVT::MvRandom (*V[k]);
70  }
71 
72  // Make timers with informative labels
73  RCP<Time> timer = TimeMonitor::getNewCounter ("Baseline for OrthoManager benchmark");
74 
75  // Baseline benchmark just copies data. It's sort of a lower
76  // bound proxy for the volume of data movement done by a real
77  // OrthoManager.
78  {
79  TimeMonitor monitor (*timer);
80  for (int trial = 0; trial < numTrials; ++trial) {
81  for (int k = 0; k < numBlocks; ++k) {
82  for (int j = 0; j < k; ++j)
83  MVT::Assign (*V[j], *V[k]);
84  MVT::Assign (*X, *V[k]);
85  }
86  }
87  }
88  }
89 
121  static void
122  benchmark (const Teuchos::RCP<OrthoManager<Scalar, MV> >& orthoMan,
123  const std::string& orthoManName,
124  const std::string& normalization,
125  const Teuchos::RCP<const MV>& X,
126  const int numCols,
127  const int numBlocks,
128  const int numTrials,
129  const Teuchos::RCP<OutputManager<Scalar> >& outMan,
130  std::ostream& resultStream,
131  const bool displayResultsCompactly=false)
132  {
133  using Teuchos::Array;
134  using Teuchos::ArrayView;
135  using Teuchos::RCP;
136  using Teuchos::rcp;
137  using Teuchos::Time;
138  using Teuchos::TimeMonitor;
139  using std::endl;
140 
141  TEUCHOS_TEST_FOR_EXCEPTION(orthoMan.is_null(), std::invalid_argument,
142  "orthoMan is null");
143  TEUCHOS_TEST_FOR_EXCEPTION(X.is_null(), std::invalid_argument,
144  "X is null");
145  TEUCHOS_TEST_FOR_EXCEPTION(numCols < 1, std::invalid_argument,
146  "numCols = " << numCols << " < 1");
147  TEUCHOS_TEST_FOR_EXCEPTION(numBlocks < 1, std::invalid_argument,
148  "numBlocks = " << numBlocks << " < 1");
149  TEUCHOS_TEST_FOR_EXCEPTION(numTrials < 1, std::invalid_argument,
150  "numTrials = " << numTrials << " < 1");
151  // Debug output stream
152  std::ostream& debugOut = outMan->stream(Debug);
153 
154  // If you like, you can add the "baseline" as an approximate
155  // lower bound for orthogonalization performance. It may be
156  // useful as a sanity check to make sure that your
157  // orthogonalizations are really computing something, though
158  // testing accuracy can help with that too.
159  //
160  //baseline (X, numCols, numBlocks, numTrials);
161 
162  // Make space to put the projection and normalization
163  // coefficients.
164  Array<RCP<mat_type> > C (numBlocks);
165  for (int k = 0; k < numBlocks; ++k) {
166  C[k] = rcp (new mat_type (numCols, numCols));
167  }
168  RCP<mat_type> B (new mat_type (numCols, numCols));
169 
170  // Make some blocks to orthogonalize. Fill with random data.
171  // We won't be orthogonalizing X, or even modifying X. We
172  // only need X so that we can make clones (since X knows its
173  // data distribution).
174  Array<RCP<MV> > V (numBlocks);
175  for (int k = 0; k < numBlocks; ++k) {
176  V[k] = MVT::Clone (*X, numCols);
177  MVT::MvRandom (*V[k]);
178  }
179 
180  // Make timers with informative labels. We time an additional
181  // first run to measure the startup costs, if any, of the
182  // OrthoManager instance.
183  RCP<Time> firstRunTimer;
184  {
185  std::ostringstream os;
186  os << "OrthoManager: " << orthoManName << " first run";
187  firstRunTimer = TimeMonitor::getNewCounter (os.str());
188  }
189  RCP<Time> timer;
190  {
191  std::ostringstream os;
192  os << "OrthoManager: " << orthoManName << " total over "
193  << numTrials << " trials (excluding first run above)";
194  timer = TimeMonitor::getNewCounter (os.str());
195  }
196  // The first run lets us measure the startup costs, if any, of
197  // the OrthoManager instance, without these costs influencing
198  // the following timing runs.
199  {
200  TimeMonitor monitor (*firstRunTimer);
201  {
202  (void) orthoMan->normalize (*V[0], B);
203  for (int k = 1; k < numBlocks; ++k) {
204  // k is the number of elements in the ArrayView. We
205  // have to assign first to an ArrayView-of-RCP-of-MV,
206  // rather than to an ArrayView-of-RCP-of-const-MV, since
207  // the latter requires a reinterpret cast. Don't you
208  // love C++ type inference?
209  ArrayView<RCP<MV> > V_0k_nonconst = V.view (0, k);
210  ArrayView<RCP<const MV> > V_0k =
211  Teuchos::av_reinterpret_cast<RCP<const MV> > (V_0k_nonconst);
212  (void) orthoMan->projectAndNormalize (*V[k], C, B, V_0k);
213  }
214  }
215  // "Test" that the trial run actually orthogonalized
216  // correctly. Results are printed to the OutputManager's
217  // Belos::Debug output stream, so depending on the
218  // OutputManager's chosen verbosity level, you may or may
219  // not see the results of the test.
220  //
221  // NOTE (mfh 22 Jan 2011) For now, these results have to be
222  // inspected visually. We should add a simple automatic
223  // test.
224  debugOut << "Orthogonality of V[0:" << (numBlocks-1)
225  << "]:" << endl;
226  for (int k = 0; k < numBlocks; ++k) {
227  // Orthogonality of each block
228  debugOut << "For block V[" << k << "]:" << endl;
229  debugOut << " ||<V[" << k << "], V[" << k << "]> - I|| = "
230  << orthoMan->orthonormError(*V[k]) << endl;
231  // Relative orthogonality with the previous blocks
232  for (int j = 0; j < k; ++j) {
233  debugOut << " ||< V[" << j << "], V[" << k << "] >|| = "
234  << orthoMan->orthogError(*V[j], *V[k]) << endl;
235  }
236  }
237  }
238 
239  // Run the benchmark for numTrials trials. Time all trials as
240  // a single run.
241  {
242  TimeMonitor monitor (*timer);
243 
244  for (int trial = 0; trial < numTrials; ++trial) {
245  (void) orthoMan->normalize (*V[0], B);
246  for (int k = 1; k < numBlocks; ++k) {
247  ArrayView<RCP<MV> > V_0k_nonconst = V.view (0, k);
248  ArrayView<RCP<const MV> > V_0k =
249  Teuchos::av_reinterpret_cast<RCP<const MV> > (V_0k_nonconst);
250  (void) orthoMan->projectAndNormalize (*V[k], C, B, V_0k);
251  }
252  }
253  }
254 
255  // Report timing results.
256  if (displayResultsCompactly)
257  {
258  // The "compact" format is suitable for automatic parsing,
259  // using a CSV (Comma-Delimited Values) parser. The first
260  // "comment" line may be parsed to extract column
261  // ("field") labels; the second line contains the actual
262  // data, in ASCII comma-delimited format.
263  using std::endl;
264  resultStream << "#orthoManName"
265  << ",normalization"
266  << ",numRows"
267  << ",numCols"
268  << ",numBlocks"
269  << ",firstRunTimeInSeconds"
270  << ",timeInSeconds"
271  << ",numTrials"
272  << endl;
273  resultStream << orthoManName
274  << "," << (orthoManName=="Simple" ? normalization : "N/A")
275  << "," << MVT::GetGlobalLength(*X)
276  << "," << numCols
277  << "," << numBlocks
278  << "," << firstRunTimer->totalElapsedTime()
279  << "," << timer->totalElapsedTime()
280  << "," << numTrials
281  << endl;
282  }
283  else {
284  TimeMonitor::summarize (resultStream);
285  }
286  }
287  };
288 
292  template< class Scalar, class MV >
293  class OrthoManagerTester {
294  private:
295  typedef typename Teuchos::Array<Teuchos::RCP<MV> >::size_type size_type;
296 
297  public:
298  typedef Scalar scalar_type;
299  typedef Teuchos::ScalarTraits<scalar_type> SCT;
300  typedef typename SCT::magnitudeType magnitude_type;
301  typedef Teuchos::ScalarTraits<magnitude_type> SMT;
302  typedef MultiVecTraits<scalar_type, MV> MVT;
303  typedef Teuchos::SerialDenseMatrix<int, scalar_type> mat_type;
304 
321  static int
322  runTests (const Teuchos::RCP<OrthoManager<Scalar, MV> >& OM,
323  const bool isRankRevealing,
324  const Teuchos::RCP<MV>& S,
325  const int sizeX1,
326  const int sizeX2,
327  const Teuchos::RCP<OutputManager<Scalar> >& MyOM)
328  {
329  using Teuchos::Array;
330  using Teuchos::null;
331  using Teuchos::RCP;
332  using Teuchos::rcp;
333  using Teuchos::rcp_dynamic_cast;
334  using Teuchos::tuple;
335 
336  // Number of tests that have failed thus far.
337  int numFailed = 0;
338 
339  // Relative tolerance against which all tests are performed.
340  const magnitude_type TOL = 1.0e-12;
341  // Absolute tolerance constant
342  //const magnitude_type ATOL = 10;
343 
344  const scalar_type ZERO = SCT::zero();
345  const scalar_type ONE = SCT::one();
346 
347  // Debug output stream
348  std::ostream& debugOut = MyOM->stream(Debug);
349 
350  // Number of columns in the input "prototype" multivector S.
351  const int sizeS = MVT::GetNumberVecs (*S);
352 
353  // Create multivectors X1 and X2, using the same map as multivector
354  // S. Then, test orthogonalizing X2 against X1. After doing so, X1
355  // and X2 should each be M-orthonormal, and should be mutually
356  // M-orthogonal.
357  debugOut << "Generating X1,X2 for testing... ";
358  RCP< MV > X1 = MVT::Clone (*S, sizeX1);
359  RCP< MV > X2 = MVT::Clone (*S, sizeX2);
360  debugOut << "done." << endl;
361  {
362  magnitude_type err;
363 
364  //
365  // Fill X1 with random values, and test the normalization error.
366  //
367  debugOut << "Filling X1 with random values... ";
368  MVT::MvRandom(*X1);
369  debugOut << "done." << endl
370  << "Calling normalize() on X1... ";
371  // The Anasazi and Belos OrthoManager interfaces differ.
372  // For example, Anasazi's normalize() method accepts either
373  // one or two arguments, whereas Belos' normalize() requires
374  // two arguments.
375  const int initialX1Rank = OM->normalize(*X1, Teuchos::null);
376  TEUCHOS_TEST_FOR_EXCEPTION(initialX1Rank != sizeX1,
377  std::runtime_error,
378  "normalize(X1) returned rank "
379  << initialX1Rank << " from " << sizeX1
380  << " vectors. Cannot continue.");
381  debugOut << "done." << endl
382  << "Calling orthonormError() on X1... ";
383  err = OM->orthonormError(*X1);
384  TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
385  "After normalize(X1), orthonormError(X1) = "
386  << err << " > TOL = " << TOL);
387  debugOut << "done: ||<X1,X1> - I|| = " << err << endl;
388 
389  //
390  // Fill X2 with random values, project against X1 and normalize,
391  // and test the orthogonalization error.
392  //
393  debugOut << "Filling X2 with random values... ";
394  MVT::MvRandom(*X2);
395  debugOut << "done." << endl
396  << "Calling projectAndNormalize(X2, C, B, tuple(X1))... "
397  << std::flush;
398  // The projectAndNormalize() interface also differs between
399  // Anasazi and Belos. Anasazi's projectAndNormalize() puts
400  // the multivector and the array of multivectors first, and
401  // the (array of) SerialDenseMatrix arguments (which are
402  // optional) afterwards. Belos puts the (array of)
403  // SerialDenseMatrix arguments in the middle, and they are
404  // not optional.
405  int initialX2Rank;
406  {
407  Array<RCP<mat_type> > C (1);
408  RCP<mat_type> B = Teuchos::null;
409  initialX2Rank =
410  OM->projectAndNormalize (*X2, C, B, tuple<RCP<const MV> >(X1));
411  }
412  TEUCHOS_TEST_FOR_EXCEPTION(initialX2Rank != sizeX2,
413  std::runtime_error,
414  "projectAndNormalize(X2,X1) returned rank "
415  << initialX2Rank << " from " << sizeX2
416  << " vectors. Cannot continue.");
417  debugOut << "done." << endl
418  << "Calling orthonormError() on X2... ";
419  err = OM->orthonormError (*X2);
420  TEUCHOS_TEST_FOR_EXCEPTION(err > TOL,
421  std::runtime_error,
422  "projectAndNormalize(X2,X1) did not meet tolerance: "
423  "orthonormError(X2) = " << err << " > TOL = " << TOL);
424  debugOut << "done: || <X2,X2> - I || = " << err << endl
425  << "Calling orthogError(X2, X1)... ";
426  err = OM->orthogError (*X2, *X1);
427  TEUCHOS_TEST_FOR_EXCEPTION(err > TOL,
428  std::runtime_error,
429  "projectAndNormalize(X2,X1) did not meet tolerance: "
430  "orthogError(X2,X1) = " << err << " > TOL = " << TOL);
431  debugOut << "done: || <X2,X1> || = " << err << endl;
432  }
433 
434 #ifdef HAVE_BELOS_TSQR
435  //
436  // If OM is an OutOfPlaceNormalizerMixin, exercise the
437  // out-of-place normalization routines.
438  //
439  typedef Belos::OutOfPlaceNormalizerMixin<Scalar, MV> mixin_type;
440  RCP<mixin_type> tsqr = rcp_dynamic_cast<mixin_type>(OM);
441  if (! tsqr.is_null())
442  {
443  magnitude_type err;
444  debugOut << endl
445  << "=== OutOfPlaceNormalizerMixin tests ==="
446  << endl << endl;
447  //
448  // Fill X1_in with random values, and test the normalization
449  // error with normalizeOutOfPlace().
450  //
451  // Don't overwrite X1, else you'll mess up the tests that
452  // follow!
453  //
454  RCP<MV> X1_in = MVT::CloneCopy (*X1);
455  debugOut << "Filling X1_in with random values... ";
456  MVT::MvRandom(*X1_in);
457  debugOut << "done." << endl;
458  debugOut << "Filling X1_out with different random values...";
459  RCP<MV> X1_out = MVT::Clone(*X1_in, MVT::GetNumberVecs(*X1_in));
460  MVT::MvRandom(*X1_out);
461  debugOut << "done." << endl
462  << "Calling normalizeOutOfPlace(*X1_in, *X1_out, null)... ";
463  const int initialX1Rank =
464  tsqr->normalizeOutOfPlace(*X1_in, *X1_out, Teuchos::null);
465  TEUCHOS_TEST_FOR_EXCEPTION(initialX1Rank != sizeX1, std::runtime_error,
466  "normalizeOutOfPlace(*X1_in, *X1_out, null) "
467  "returned rank " << initialX1Rank << " from "
468  << sizeX1 << " vectors. Cannot continue.");
469  debugOut << "done." << endl
470  << "Calling orthonormError() on X1_out... ";
471  err = OM->orthonormError(*X1_out);
472  TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
473  "After calling normalizeOutOfPlace(*X1_in, "
474  "*X1_out, null), orthonormError(X1) = "
475  << err << " > TOL = " << TOL);
476  debugOut << "done: ||<X1_out,X1_out> - I|| = " << err << endl;
477 
478  //
479  // Fill X2_in with random values, project against X1_out
480  // and normalize via projectAndNormalizeOutOfPlace(), and
481  // test the orthogonalization error.
482  //
483  // Don't overwrite X2, else you'll mess up the tests that
484  // follow!
485  //
486  RCP<MV> X2_in = MVT::CloneCopy (*X2);
487  debugOut << "Filling X2_in with random values... ";
488  MVT::MvRandom(*X2_in);
489  debugOut << "done." << endl
490  << "Filling X2_out with different random values...";
491  RCP<MV> X2_out = MVT::Clone(*X2_in, MVT::GetNumberVecs(*X2_in));
492  MVT::MvRandom(*X2_out);
493  debugOut << "done." << endl
494  << "Calling projectAndNormalizeOutOfPlace(X2_in, X2_out, "
495  << "C, B, X1_out)...";
496  int initialX2Rank;
497  {
498  Array<RCP<mat_type> > C (1);
499  RCP<mat_type> B = Teuchos::null;
500  initialX2Rank =
501  tsqr->projectAndNormalizeOutOfPlace (*X2_in, *X2_out, C, B,
502  tuple<RCP<const MV> >(X1_out));
503  }
504  TEUCHOS_TEST_FOR_EXCEPTION(initialX2Rank != sizeX2,
505  std::runtime_error,
506  "projectAndNormalizeOutOfPlace(*X2_in, "
507  "*X2_out, C, B, tuple(X1_out)) returned rank "
508  << initialX2Rank << " from " << sizeX2
509  << " vectors. Cannot continue.");
510  debugOut << "done." << endl
511  << "Calling orthonormError() on X2_out... ";
512  err = OM->orthonormError (*X2_out);
513  TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
514  "projectAndNormalizeOutOfPlace(*X2_in, *X2_out, "
515  "C, B, tuple(X1_out)) did not meet tolerance: "
516  "orthonormError(X2_out) = "
517  << err << " > TOL = " << TOL);
518  debugOut << "done: || <X2_out,X2_out> - I || = " << err << endl
519  << "Calling orthogError(X2_out, X1_out)... ";
520  err = OM->orthogError (*X2_out, *X1_out);
521  TEUCHOS_TEST_FOR_EXCEPTION(err > TOL, std::runtime_error,
522  "projectAndNormalizeOutOfPlace(*X2_in, *X2_out, "
523  "C, B, tuple(X1_out)) did not meet tolerance: "
524  "orthogError(X2_out, X1_out) = "
525  << err << " > TOL = " << TOL);
526  debugOut << "done: || <X2_out,X1_out> || = " << err << endl;
527  debugOut << endl
528  << "=== Done with OutOfPlaceNormalizerMixin tests ==="
529  << endl << endl;
530  }
531 #endif // HAVE_BELOS_TSQR
532 
533  {
534  //
535  // Test project() on a random multivector S, by projecting S
536  // against various combinations of X1 and X2.
537  //
538  MVT::MvRandom(*S);
539 
540  debugOut << "Testing project() by projecting a random multivector S "
541  "against various combinations of X1 and X2 " << endl;
542  const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
543  numFailed += thisNumFailed;
544  if (thisNumFailed > 0)
545  debugOut << " *** " << thisNumFailed
546  << (thisNumFailed > 1 ? " tests" : " test")
547  << " failed." << endl;
548  }
549 
550  {
551  //
552  // Test normalize() for various deficient cases
553  //
554  debugOut << "Testing normalize() on bad multivectors " << endl;
555  const int thisNumFailed = testNormalize(OM,S,MyOM);
556  numFailed += thisNumFailed;
557  }
558 
559  if (isRankRevealing)
560  {
561  // run a X1,Y2 range multivector against P_{X1,X1} P_{Y2,Y2}
562  // note, this is allowed under the restrictions on project(),
563  // because <X1,Y2> = 0
564  // also, <Y2,Y2> = I, but <X1,X1> != I, so biOrtho must be set to false
565  // it should require randomization, as
566  // P_{X1,X1} P_{Y2,Y2} (X1*C1 + Y2*C2) = P_{X1,X1} X1*C1 = 0
567  mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
568  Teuchos::randomSyncedMatrix(C1);
569  Teuchos::randomSyncedMatrix(C2);
570  // S := X1*C1
571  MVT::MvTimesMatAddMv(ONE,*X1,C1,ZERO,*S);
572  // S := S + X2*C2
573  MVT::MvTimesMatAddMv(ONE,*X2,C2,ONE,*S);
574 
575  debugOut << "Testing project() by projecting [X1 X2]-range multivector "
576  "against P_X1 P_X2 " << endl;
577  const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
578  numFailed += thisNumFailed;
579  if (thisNumFailed > 0)
580  debugOut << " *** " << thisNumFailed
581  << (thisNumFailed > 1 ? " tests" : " test")
582  << " failed." << endl;
583  }
584 
585  // This test is only distinct from the rank-1 multivector test
586  // (below) if S has at least 3 columns.
587  if (isRankRevealing && sizeS > 2)
588  {
589  MVT::MvRandom(*S);
590  RCP<MV> mid = MVT::Clone(*S,1);
591  mat_type c(sizeS,1);
592  MVT::MvTimesMatAddMv(ONE,*S,c,ZERO,*mid);
593  std::vector<int> ind(1);
594  ind[0] = sizeS-1;
595  MVT::SetBlock(*mid,ind,*S);
596 
597  debugOut << "Testing normalize() on a rank-deficient multivector " << endl;
598  const int thisNumFailed = testNormalizeRankReveal(OM,S,MyOM);
599  numFailed += thisNumFailed;
600  if (thisNumFailed > 0)
601  debugOut << " *** " << thisNumFailed
602  << (thisNumFailed > 1 ? " tests" : " test")
603  << " failed." << endl;
604  }
605 
606  // This test will only exercise rank deficiency if S has at least 2
607  // columns.
608  if (isRankRevealing && sizeS > 1)
609  {
610  // rank-1
611  RCP<MV> one = MVT::Clone(*S,1);
612  MVT::MvRandom(*one);
613  mat_type scaleS(sizeS,1);
614  Teuchos::randomSyncedMatrix(scaleS);
615  // put multiple of column 0 in columns 0:sizeS-1
616  for (int i=0; i<sizeS; i++)
617  {
618  std::vector<int> ind(1);
619  ind[0] = i;
620  RCP<MV> Si = MVT::CloneViewNonConst(*S,ind);
621  MVT::MvAddMv(scaleS(i,0),*one,ZERO,*one,*Si);
622  }
623  debugOut << "Testing normalize() on a rank-1 multivector " << endl;
624  const int thisNumFailed = testNormalizeRankReveal(OM,S,MyOM);
625  numFailed += thisNumFailed;
626  if (thisNumFailed > 0)
627  debugOut << " *** " << thisNumFailed
628  << (thisNumFailed > 1 ? " tests" : " test")
629  << " failed." << endl;
630  }
631 
632  {
633  std::vector<int> ind(1);
634  MVT::MvRandom(*S);
635 
636  debugOut << "Testing projectAndNormalize() on a random multivector " << endl;
637  const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
638  numFailed += thisNumFailed;
639  if (thisNumFailed > 0)
640  debugOut << " *** " << thisNumFailed
641  << (thisNumFailed > 1 ? " tests" : " test")
642  << " failed." << endl;
643  }
644 
645  if (isRankRevealing)
646  {
647  // run a X1,X2 range multivector against P_X1 P_X2
648  // this is allowed as <X1,X2> == 0
649  // it should require randomization, as
650  // P_X1 P_X2 (X1*C1 + X2*C2) = P_X1 X1*C1 = 0
651  // and
652  // P_X2 P_X1 (X2*C2 + X1*C1) = P_X2 X2*C2 = 0
653  mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
654  Teuchos::randomSyncedMatrix(C1);
655  Teuchos::randomSyncedMatrix(C2);
656  MVT::MvTimesMatAddMv(ONE,*X1,C1,ZERO,*S);
657  MVT::MvTimesMatAddMv(ONE,*X2,C2,ONE,*S);
658 
659  debugOut << "Testing projectAndNormalize() by projecting [X1 X2]-range "
660  "multivector against P_X1 P_X2 " << endl;
661  const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
662  numFailed += thisNumFailed;
663  if (thisNumFailed > 0)
664  debugOut << " *** " << thisNumFailed
665  << (thisNumFailed > 1 ? " tests" : " test")
666  << " failed." << endl;
667  }
668 
669  // This test is only distinct from the rank-1 multivector test
670  // (below) if S has at least 3 columns.
671  if (isRankRevealing && sizeS > 2)
672  {
673  MVT::MvRandom(*S);
674  RCP<MV> mid = MVT::Clone(*S,1);
675  mat_type c(sizeS,1);
676  MVT::MvTimesMatAddMv(ONE,*S,c,ZERO,*mid);
677  std::vector<int> ind(1);
678  ind[0] = sizeS-1;
679  MVT::SetBlock(*mid,ind,*S);
680 
681  debugOut << "Testing projectAndNormalize() on a rank-deficient "
682  "multivector " << endl;
683  const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
684  numFailed += thisNumFailed;
685  if (thisNumFailed > 0)
686  debugOut << " *** " << thisNumFailed
687  << (thisNumFailed > 1 ? " tests" : " test")
688  << " failed." << endl;
689  }
690 
691  // This test will only exercise rank deficiency if S has at least 2
692  // columns.
693  if (isRankRevealing && sizeS > 1)
694  {
695  // rank-1
696  RCP<MV> one = MVT::Clone(*S,1);
697  MVT::MvRandom(*one);
698  mat_type scaleS(sizeS,1);
699  Teuchos::randomSyncedMatrix(scaleS);
700  // Put a multiple of column 0 in columns 0:sizeS-1.
701  for (int i=0; i<sizeS; i++)
702  {
703  std::vector<int> ind(1);
704  ind[0] = i;
705  RCP<MV> Si = MVT::CloneViewNonConst(*S,ind);
706  MVT::MvAddMv(scaleS(i,0),*one,ZERO,*one,*Si);
707  }
708  debugOut << "Testing projectAndNormalize() on a rank-1 multivector " << endl;
709  bool constantStride = true;
710  if (! MVT::HasConstantStride(*S)) {
711  debugOut << "-- S does not have constant stride" << endl;
712  constantStride = false;
713  }
714  if (! MVT::HasConstantStride(*X1)) {
715  debugOut << "-- X1 does not have constant stride" << endl;
716  constantStride = false;
717  }
718  if (! MVT::HasConstantStride(*X2)) {
719  debugOut << "-- X2 does not have constant stride" << endl;
720  constantStride = false;
721  }
722  if (! constantStride) {
723  debugOut << "-- Skipping this test, since TSQR does not work on "
724  "multivectors with nonconstant stride" << endl;
725  }
726  else {
727  const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
728  numFailed += thisNumFailed;
729  if (thisNumFailed > 0) {
730  debugOut << " *** " << thisNumFailed
731  << (thisNumFailed > 1 ? " tests" : " test")
732  << " failed." << endl;
733  }
734  }
735  }
736 
737  if (numFailed != 0) {
738  MyOM->stream(Errors) << numFailed << " total test failures." << endl;
739  }
740  return numFailed;
741  }
742 
743  private:
744 
749  static magnitude_type
750  MVDiff (const MV& X, const MV& Y)
751  {
752  using Teuchos::RCP;
753 
754  const scalar_type ONE = SCT::one();
755  const int numCols = MVT::GetNumberVecs(X);
756  TEUCHOS_TEST_FOR_EXCEPTION( (MVT::GetNumberVecs(Y) != numCols),
757  std::logic_error,
758  "MVDiff: X and Y should have the same number of columns."
759  " X has " << numCols << " column(s) and Y has "
760  << MVT::GetNumberVecs(Y) << " columns." );
761  // Resid := X
762  RCP< MV > Resid = MVT::CloneCopy(X);
763  // Resid := Resid - Y
764  MVT::MvAddMv (-ONE, Y, ONE, *Resid, *Resid);
765 
766  return frobeniusNorm (*Resid);
767  }
768 
769 
773  static magnitude_type
774  frobeniusNorm (const MV& X)
775  {
776  const scalar_type ONE = SCT::one();
777  const int numCols = MVT::GetNumberVecs(X);
778  mat_type C (numCols, numCols);
779 
780  // $C := X^* X$
781  MVT::MvTransMv (ONE, X, X, C);
782 
783  magnitude_type err (0);
784  for (int i = 0; i < numCols; ++i)
785  err += SCT::magnitude (C(i,i));
786 
787  return SCT::magnitude (SCT::squareroot (err));
788  }
789 
790 
791  static int
792  testProjectAndNormalize (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
793  const Teuchos::RCP< const MV >& S,
794  const Teuchos::RCP< const MV >& X1,
795  const Teuchos::RCP< const MV >& X2,
796  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
797  {
798  return testProjectAndNormalizeNew (OM, S, X1, X2, MyOM);
799  }
800 
805  static int
806  testProjectAndNormalizeOld (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > >& OM,
807  const Teuchos::RCP< const MV >& S,
808  const Teuchos::RCP< const MV >& X1,
809  const Teuchos::RCP< const MV >& X2,
810  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
811  {
812  using Teuchos::Array;
813  using Teuchos::null;
814  using Teuchos::RCP;
815  using Teuchos::rcp;
816  using Teuchos::tuple;
817 
818  const scalar_type ONE = SCT::one();
819  const magnitude_type ZERO = SCT::magnitude(SCT::zero());
820 
821  // Relative tolerance against which all tests are performed.
822  const magnitude_type TOL = 1.0e-12;
823  // Absolute tolerance constant
824  const magnitude_type ATOL = 10;
825 
826  const int sizeS = MVT::GetNumberVecs(*S);
827  const int sizeX1 = MVT::GetNumberVecs(*X1);
828  const int sizeX2 = MVT::GetNumberVecs(*X2);
829  int numerr = 0;
830  std::ostringstream sout;
831 
832  //
833  // output tests:
834  // <S_out,S_out> = I
835  // <S_out,X1> = 0
836  // <S_out,X2> = 0
837  // S_in = S_out B + X1 C1 + X2 C2
838  //
839  // we will loop over an integer specifying the test combinations
840  // the bit pattern for the different tests is listed in parenthesis
841  //
842  // for the projectors, test the following combinations:
843  // none (00)
844  // P_X1 (01)
845  // P_X2 (10)
846  // P_X1 P_X2 (11)
847  // P_X2 P_X1 (11)
848  // the latter two should be tested to give the same answer
849  //
850  // for each of these, we should test with C1, C2 and B
851  //
852  // if hasM:
853  // with and without MX1 (1--)
854  // with and without MX2 (1---)
855  // with and without MS (1----)
856  //
857  // as hasM controls the upper level bits, we need only run test cases 0-3 if hasM==false
858  // otherwise, we run test cases 0-31
859  //
860 
861  int numtests = 4;
862 
863  // test ortho error before orthonormalizing
864  if (X1 != null) {
865  magnitude_type err = OM->orthogError(*S,*X1);
866  sout << " || <S,X1> || before : " << err << endl;
867  }
868  if (X2 != null) {
869  magnitude_type err = OM->orthogError(*S,*X2);
870  sout << " || <S,X2> || before : " << err << endl;
871  }
872 
873  for (int t=0; t<numtests; t++) {
874 
875  Array< RCP< const MV > > theX;
876  RCP<mat_type > B = rcp( new mat_type(sizeS,sizeS) );
877  Array<RCP<mat_type > > C;
878  if ( (t % 3) == 0 ) {
879  // neither <X1,Y1> nor <X2,Y2>
880  // C, theX and theY are already empty
881  }
882  else if ( (t % 3) == 1 ) {
883  // X1
884  theX = tuple(X1);
885  C = tuple( rcp(new mat_type(sizeX1,sizeS)) );
886  }
887  else if ( (t % 3) == 2 ) {
888  // X2
889  theX = tuple(X2);
890  C = tuple( rcp(new mat_type(sizeX2,sizeS)) );
891  }
892  else {
893  // X1 and X2, and the reverse.
894  theX = tuple(X1,X2);
895  C = tuple( rcp(new mat_type(sizeX1,sizeS)),
896  rcp(new mat_type(sizeX2,sizeS)) );
897  }
898 
899  // We wrap up all the OrthoManager calls in a try-catch
900  // block, in order to check whether any of the methods throw
901  // an exception. For the tests we perform, every thrown
902  // exception is a failure.
903  try {
904  // call routine
905  // if (t && 3) == 3, {
906  // call with reversed input: X2 X1
907  // }
908  // test all outputs for correctness
909  // test all outputs for equivalence
910 
911  // here is where the outputs go
912  Array<RCP<MV> > S_outs;
913  Array<Array<RCP<mat_type > > > C_outs;
914  Array<RCP<mat_type > > B_outs;
915  RCP<MV> Scopy;
916  Array<int> ret_out;
917 
918  // copies of S,MS
919  Scopy = MVT::CloneCopy(*S);
920  // randomize this data, it should be overwritten
921  Teuchos::randomSyncedMatrix(*B);
922  for (size_type i=0; i<C.size(); i++) {
923  Teuchos::randomSyncedMatrix(*C[i]);
924  }
925  // Run test. Since S was specified by the caller and
926  // Scopy is a copy of S, we don't know what rank to expect
927  // here -- though we do require that S have rank at least
928  // one.
929  //
930  // Note that Anasazi and Belos differ, among other places,
931  // in the order of arguments to projectAndNormalize().
932  int ret = OM->projectAndNormalize(*Scopy,C,B,theX);
933  sout << "projectAndNormalize() returned rank " << ret << endl;
934  if (ret == 0) {
935  sout << " *** Error: returned rank is zero, cannot continue tests" << endl;
936  numerr++;
937  break;
938  }
939  ret_out.push_back(ret);
940  // projectAndNormalize() is only required to return a
941  // basis of rank "ret"
942  // this is what we will test:
943  // the first "ret" columns in Scopy
944  // the first "ret" rows in B
945  // save just the parts that we want
946  // we allocate S and MS for each test, so we can save these as views
947  // however, save copies of the C and B
948  if (ret < sizeS) {
949  std::vector<int> ind(ret);
950  for (int i=0; i<ret; i++) {
951  ind[i] = i;
952  }
953  S_outs.push_back( MVT::CloneViewNonConst(*Scopy,ind) );
954  B_outs.push_back( rcp( new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
955  }
956  else {
957  S_outs.push_back( Scopy );
958  B_outs.push_back( rcp( new mat_type(*B) ) );
959  }
960  C_outs.push_back( Array<RCP<mat_type > >(0) );
961  if (C.size() > 0) {
962  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
963  }
964  if (C.size() > 1) {
965  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
966  }
967 
968  // do we run the reversed input?
969  if ( (t % 3) == 3 ) {
970  // copies of S,MS
971  Scopy = MVT::CloneCopy(*S);
972 
973  // Fill the B and C[i] matrices with random data. The
974  // data will be overwritten by projectAndNormalize().
975  // Filling these matrices here is only to catch some
976  // bugs in projectAndNormalize().
977  Teuchos::randomSyncedMatrix(*B);
978  for (size_type i=0; i<C.size(); i++) {
979  Teuchos::randomSyncedMatrix(*C[i]);
980  }
981  // flip the inputs
982  theX = tuple( theX[1], theX[0] );
983  // Run test.
984  // Note that Anasazi and Belos differ, among other places,
985  // in the order of arguments to projectAndNormalize().
986  ret = OM->projectAndNormalize(*Scopy,C,B,theX);
987  sout << "projectAndNormalize() returned rank " << ret << endl;
988  if (ret == 0) {
989  sout << " *** Error: returned rank is zero, cannot continue tests" << endl;
990  numerr++;
991  break;
992  }
993  ret_out.push_back(ret);
994  // projectAndNormalize() is only required to return a
995  // basis of rank "ret"
996  // this is what we will test:
997  // the first "ret" columns in Scopy
998  // the first "ret" rows in B
999  // save just the parts that we want
1000  // we allocate S and MS for each test, so we can save these as views
1001  // however, save copies of the C and B
1002  if (ret < sizeS) {
1003  std::vector<int> ind(ret);
1004  for (int i=0; i<ret; i++) {
1005  ind[i] = i;
1006  }
1007  S_outs.push_back( MVT::CloneViewNonConst(*Scopy,ind) );
1008  B_outs.push_back( rcp( new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
1009  }
1010  else {
1011  S_outs.push_back( Scopy );
1012  B_outs.push_back( rcp( new mat_type(*B) ) );
1013  }
1014  C_outs.push_back( Array<RCP<mat_type > >() );
1015  // reverse the Cs to compensate for the reverse projectors
1016  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
1017  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
1018  // flip the inputs back
1019  theX = tuple( theX[1], theX[0] );
1020  }
1021 
1022 
1023  // test all outputs for correctness
1024  for (size_type o=0; o<S_outs.size(); o++) {
1025  // S^T M S == I
1026  {
1027  magnitude_type err = OM->orthonormError(*S_outs[o]);
1028  if (err > TOL) {
1029  sout << endl
1030  << " *** Test (number " << (t+1) << " of " << numtests
1031  << " total tests) failed: Tolerance exceeded! Error = "
1032  << err << " > TOL = " << TOL << "."
1033  << endl << endl;
1034  numerr++;
1035  }
1036  sout << " || <S,S> - I || after : " << err << endl;
1037  }
1038  // S_in = X1*C1 + C2*C2 + S_out*B
1039  {
1040  RCP<MV> tmp = MVT::Clone(*S,sizeS);
1041  MVT::MvTimesMatAddMv(ONE,*S_outs[o],*B_outs[o],ZERO,*tmp);
1042  if (C_outs[o].size() > 0) {
1043  MVT::MvTimesMatAddMv(ONE,*X1,*C_outs[o][0],ONE,*tmp);
1044  if (C_outs[o].size() > 1) {
1045  MVT::MvTimesMatAddMv(ONE,*X2,*C_outs[o][1],ONE,*tmp);
1046  }
1047  }
1048  magnitude_type err = MVDiff(*tmp,*S);
1049  if (err > ATOL*TOL) {
1050  sout << endl
1051  << " *** Test (number " << (t+1) << " of " << numtests
1052  << " total tests) failed: Tolerance exceeded! Error = "
1053  << err << " > ATOL*TOL = " << (ATOL*TOL) << "."
1054  << endl << endl;
1055  numerr++;
1056  }
1057  sout << " " << t << "|| S_in - X1*C1 - X2*C2 - S_out*B || : " << err << endl;
1058  }
1059  // <X1,S> == 0
1060  if (theX.size() > 0 && theX[0] != null) {
1061  magnitude_type err = OM->orthogError(*theX[0],*S_outs[o]);
1062  if (err > TOL) {
1063  sout << endl
1064  << " *** Test (number " << (t+1) << " of " << numtests
1065  << " total tests) failed: Tolerance exceeded! Error = "
1066  << err << " > TOL = " << TOL << "."
1067  << endl << endl;
1068  numerr++;
1069  }
1070  sout << " " << t << "|| <X[0],S> || after : " << err << endl;
1071  }
1072  // <X2,S> == 0
1073  if (theX.size() > 1 && theX[1] != null) {
1074  magnitude_type err = OM->orthogError(*theX[1],*S_outs[o]);
1075  if (err > TOL) {
1076  sout << endl
1077  << " *** Test (number " << (t+1) << " of " << numtests
1078  << " total tests) failed: Tolerance exceeded! Error = "
1079  << err << " > TOL = " << TOL << "."
1080  << endl << endl;
1081  numerr++;
1082  }
1083  sout << " " << t << "|| <X[1],S> || after : " << err << endl;
1084  }
1085  }
1086  }
1087  catch (Belos::OrthoError& e) {
1088  sout << " *** Error: OrthoManager threw exception: " << e.what() << endl;
1089  numerr++;
1090  }
1091 
1092  } // test for
1093 
1094  // NOTE (mfh 05 Nov 2010) Since Belos::MsgType is an enum,
1095  // doing bitwise logical computations on Belos::MsgType values
1096  // (such as "Debug | Errors") and passing the result into
1097  // MyOM->stream() confuses the compiler. As a result, we have
1098  // to do some type casts to make it work.
1099  const int msgType = (numerr > 0) ?
1100  (static_cast<int>(Debug) | static_cast<int>(Errors)) :
1101  static_cast<int>(Debug);
1102 
1103  // We report debug-level messages always. We also report
1104  // errors if at least one test failed.
1105  MyOM->stream(static_cast< MsgType >(msgType)) << sout.str() << endl;
1106  return numerr;
1107  }
1108 
1113  static int
1114  testNormalize (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > >& OM,
1115  const Teuchos::RCP< const MV >& S,
1116  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
1117  {
1118  using Teuchos::RCP;
1119 
1120  int numFailures = 0;
1121  const scalar_type ZERO = SCT::zero();
1122 
1123  const int msgType = (static_cast<int>(Debug) | static_cast<int>(Errors));
1124 
1125  // Check that the orthogonalization gracefully handles zero vectors.
1126  RCP<MV> zeroVec = MVT::Clone(*S,1);
1127  RCP< mat_type > bZero (new mat_type (1, 1));
1128  std::vector< magnitude_type > zeroNorm( 1 );
1129 
1130  MVT::MvInit( *zeroVec, ZERO );
1131  OM->normalize( *zeroVec, bZero );
1132  MVT::MvNorm( *zeroVec, zeroNorm );
1133  // Check if the number is a NaN, this orthogonalization fails if it is.
1134  if ( zeroNorm[0] != ZERO )
1135  {
1136  MyOM->stream(static_cast< MsgType >(msgType)) << " --> Normalization of zero vector FAILED!" << std::endl;
1137  numFailures++;
1138  }
1139 
1140  return numFailures;
1141  }
1142 
1147  static int
1148  testNormalizeRankReveal (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > >& OM,
1149  const Teuchos::RCP< const MV >& S,
1150  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
1151  {
1152  using Teuchos::Array;
1153  using Teuchos::RCP;
1154  using Teuchos::rcp;
1155  using Teuchos::tuple;
1156 
1157  const scalar_type ONE = SCT::one();
1158  std::ostringstream sout;
1159  // Total number of failed tests in this call of this routine.
1160  int numerr = 0;
1161 
1162  // Relative tolerance against which all tests are performed.
1163  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1164  // The following bounds hold for all $m \times n$ matrices $A$:
1165  // \[
1166  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1167  // \]
1168  // where $r$ is the (column) rank of $A$. We bound this above
1169  // by the number of columns in $A$.
1170  //
1171  // An accurate normalization in the Euclidean norm of a matrix
1172  // $A$ with at least as many rows m as columns n, should
1173  // produce orthogonality $\|Q^* Q - I\|_2$ less than a factor
1174  // of machine precision times a low-order polynomial in m and
1175  // n, and residual $\|A - Q B\|_2$ (where $A = Q B$ is the
1176  // computed normalization) less than that bound times the norm
1177  // of $A$.
1178  //
1179  // Since we are measuring both of these quantitites in the
1180  // Frobenius norm instead, we should scale this bound by
1181  // $\sqrt{n}$.
1182 
1183  const int numRows = MVT::GetGlobalLength(*S);
1184  const int numCols = MVT::GetNumberVecs(*S);
1185  const int sizeS = MVT::GetNumberVecs(*S);
1186 
1187  // A good heuristic is to scale the bound by the square root
1188  // of the number of floating-point operations. One could
1189  // perhaps support this theoretically, since we are using
1190  // uniform random test problems.
1191  const magnitude_type fudgeFactor =
1192  SMT::squareroot(magnitude_type(numRows) *
1193  magnitude_type(numCols) *
1194  magnitude_type(numCols));
1195  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1196  SMT::squareroot(magnitude_type(numCols));
1197 
1198  // Absolute tolerance scaling: the Frobenius norm of the test
1199  // matrix S. TOL*ATOL is the absolute tolerance for the
1200  // residual $\|A - Q*B\|_F$.
1201  const magnitude_type ATOL = frobeniusNorm (*S);
1202 
1203  sout << "The test matrix S has Frobenius norm " << ATOL
1204  << ", and the relative error tolerance is TOL = "
1205  << TOL << "." << endl;
1206 
1207  const int numtests = 1;
1208  for (int t = 0; t < numtests; ++t) {
1209 
1210  try {
1211  // call routine
1212  // test all outputs for correctness
1213 
1214  // S_copy gets a copy of S; we normalize in place, so we
1215  // need a copy to check whether the normalization
1216  // succeeded.
1217  RCP< MV > S_copy = MVT::CloneCopy (*S);
1218 
1219  // Matrix of coefficients from the normalization.
1220  RCP< mat_type > B (new mat_type (sizeS, sizeS));
1221  // The contents of B will be overwritten, but fill with
1222  // random data just to make sure that the normalization
1223  // operated on all the elements of B on which it should
1224  // operate.
1225  Teuchos::randomSyncedMatrix(*B);
1226 
1227  const int reportedRank = OM->normalize (*S_copy, B);
1228  sout << "normalize() returned rank " << reportedRank << endl;
1229  if (reportedRank == 0) {
1230  sout << " *** Error: Cannot continue, since normalize() "
1231  "reports that S has rank 0" << endl;
1232  numerr++;
1233  break;
1234  }
1235  //
1236  // We don't know in this routine whether the input
1237  // multivector S has full rank; it is only required to
1238  // have nonzero rank. Thus, we extract the first
1239  // reportedRank columns of S_copy and the first
1240  // reportedRank rows of B, and perform tests on them.
1241  //
1242 
1243  // Construct S_view, a view of the first reportedRank
1244  // columns of S_copy.
1245  std::vector<int> indices (reportedRank);
1246  for (int j = 0; j < reportedRank; ++j)
1247  indices[j] = j;
1248  RCP< MV > S_view = MVT::CloneViewNonConst (*S_copy, indices);
1249  // Construct B_top, a copy of the first reportedRank rows
1250  // of B.
1251  //
1252  // NOTE: We create this as a copy and not a view, because
1253  // otherwise it would not be safe with respect to RCPs.
1254  // This is because mat_type uses raw pointers
1255  // inside, so that a view would become invalid when B
1256  // would fall out of scope.
1257  RCP< mat_type > B_top (new mat_type (Teuchos::Copy, *B, reportedRank, sizeS));
1258 
1259  // Check ||<S_view,S_view> - I||
1260  {
1261  const magnitude_type err = OM->orthonormError(*S_view);
1262  if (err > TOL) {
1263  sout << " *** Error: Tolerance exceeded: err = "
1264  << err << " > TOL = " << TOL << endl;
1265  numerr++;
1266  }
1267  sout << " || <S,S> - I || after : " << err << endl;
1268  }
1269  // Check the residual ||Residual|| = ||S_view * B_top -
1270  // S_orig||, where S_orig is a view of the first
1271  // reportedRank columns of S.
1272  {
1273  // Residual is allocated with reportedRank columns. It
1274  // will contain the result of testing the residual error
1275  // of the normalization (i.e., $\|S - S_in*B\|$). It
1276  // should have the dimensions of S. Its initial value
1277  // is a copy of the first reportedRank columns of S.
1278  RCP< MV > Residual = MVT::CloneCopy (*S);
1279 
1280  // Residual := Residual - S_view * B_view
1281  MVT::MvTimesMatAddMv (-ONE, *S_view, *B_top, ONE, *Residual);
1282 
1283  // Compute ||Residual||
1284  const magnitude_type err = frobeniusNorm (*Residual);
1285  if (err > ATOL*TOL) {
1286  sout << " *** Error: Tolerance exceeded: err = "
1287  << err << " > ATOL*TOL = " << (ATOL*TOL) << endl;
1288  numerr++;
1289  }
1290  sout << " " << t << "|| S - Q*B || : " << err << endl;
1291  }
1292  }
1293  catch (Belos::OrthoError& e) {
1294  sout << " *** Error: the OrthoManager's normalize() method "
1295  "threw an exception: " << e.what() << endl;
1296  numerr++;
1297  }
1298 
1299  } // test for
1300 
1301  const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
1302  MyOM->stream(type) << sout.str();
1303  MyOM->stream(type) << endl;
1304 
1305  return numerr;
1306  }
1307 
1312  static int
1313  testProjectAndNormalizeNew (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1314  const Teuchos::RCP< const MV >& S,
1315  const Teuchos::RCP< const MV >& X1,
1316  const Teuchos::RCP< const MV >& X2,
1317  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
1318  {
1319  using Teuchos::Array;
1320  using Teuchos::null;
1321  using Teuchos::RCP;
1322  using Teuchos::rcp;
1323  using Teuchos::tuple;
1324 
1325  // We collect all the output in this string wrapper, and print
1326  // it at the end.
1327  std::ostringstream sout;
1328  // Total number of failed tests in this call of this routine.
1329  int numerr = 0;
1330 
1331  const int numRows = MVT::GetGlobalLength(*S);
1332  const int numCols = MVT::GetNumberVecs(*S);
1333  const int sizeS = MVT::GetNumberVecs(*S);
1334 
1335  // Relative tolerance against which all tests are performed.
1336  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1337  // The following bounds hold for all $m \times n$ matrices $A$:
1338  // \[
1339  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1340  // \]
1341  // where $r$ is the (column) rank of $A$. We bound this above
1342  // by the number of columns in $A$.
1343  //
1344  // Since we are measuring both of these quantitites in the
1345  // Frobenius norm instead, we scale all error tests by
1346  // $\sqrt{n}$.
1347  //
1348  // A good heuristic is to scale the bound by the square root
1349  // of the number of floating-point operations. One could
1350  // perhaps support this theoretically, since we are using
1351  // uniform random test problems.
1352  const magnitude_type fudgeFactor =
1353  SMT::squareroot(magnitude_type(numRows) *
1354  magnitude_type(numCols) *
1355  magnitude_type(numCols));
1356  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1357  SMT::squareroot(magnitude_type(numCols));
1358 
1359  // Absolute tolerance scaling: the Frobenius norm of the test
1360  // matrix S. TOL*ATOL is the absolute tolerance for the
1361  // residual $\|A - Q*B\|_F$.
1362  const magnitude_type ATOL = frobeniusNorm (*S);
1363 
1364  sout << "-- The test matrix S has Frobenius norm " << ATOL
1365  << ", and the relative error tolerance is TOL = "
1366  << TOL << "." << endl;
1367 
1368  // Q will contain the result of projectAndNormalize() on S.
1369  RCP< MV > Q = MVT::CloneCopy(*S);
1370  // We use this for collecting the residual error components
1371  RCP< MV > Residual = MVT::CloneCopy(*S);
1372  // Number of elements in the X array of blocks against which
1373  // to project S.
1374  const int num_X = 2;
1375  Array< RCP< const MV > > X (num_X);
1376  X[0] = MVT::CloneCopy(*X1);
1377  X[1] = MVT::CloneCopy(*X2);
1378 
1379  // Coefficients for the normalization
1380  RCP< mat_type > B (new mat_type (sizeS, sizeS));
1381 
1382  // Array of coefficients matrices from the projection.
1383  // For our first test, we allocate each of these matrices
1384  // with the proper dimensions.
1385  Array< RCP< mat_type > > C (num_X);
1386  for (int k = 0; k < num_X; ++k)
1387  {
1388  C[k] = rcp (new mat_type (MVT::GetNumberVecs(*X[k]), sizeS));
1389  Teuchos::randomSyncedMatrix(*C[k]); // will be overwritten
1390  }
1391  try {
1392  // Q*B := (I - X X^*) S
1393  const int reportedRank = OM->projectAndNormalize (*Q, C, B, X);
1394 
1395  // Pick out the first reportedRank columns of Q.
1396  std::vector<int> indices (reportedRank);
1397  for (int j = 0; j < reportedRank; ++j)
1398  indices[j] = j;
1399  RCP< const MV > Q_left = MVT::CloneView (*Q, indices);
1400 
1401  // Test whether the first reportedRank columns of Q are
1402  // orthogonal.
1403  {
1404  const magnitude_type orthoError = OM->orthonormError (*Q_left);
1405  sout << "-- ||Q(1:" << reportedRank << ")^* Q(1:" << reportedRank
1406  << ") - I||_F = " << orthoError << endl;
1407  if (orthoError > TOL)
1408  {
1409  sout << " *** Error: ||Q(1:" << reportedRank << ")^* Q(1:"
1410  << reportedRank << ") - I||_F = " << orthoError
1411  << " > TOL = " << TOL << "." << endl;
1412  numerr++;
1413  }
1414  }
1415 
1416  // Compute the residual: if successful, S = Q*B +
1417  // X (X^* S =: C) in exact arithmetic. So, the residual is
1418  // S - Q*B - X1 C1 - X2 C2.
1419  //
1420  // Residual := S
1421  MVT::MvAddMv (SCT::one(), *S, SCT::zero(), *Residual, *Residual);
1422  {
1423  // Pick out the first reportedRank rows of B. Make a deep
1424  // copy, since mat_type is not safe with respect
1425  // to RCP-based memory management (it uses raw pointers
1426  // inside).
1427  RCP< const mat_type > B_top (new mat_type (Teuchos::Copy, *B, reportedRank, B->numCols()));
1428  // Residual := Residual - Q(:, 1:reportedRank) * B(1:reportedRank, :)
1429  MVT::MvTimesMatAddMv (-SCT::one(), *Q_left, *B_top, SCT::one(), *Residual);
1430  }
1431  // Residual := Residual - X[k]*C[k]
1432  for (int k = 0; k < num_X; ++k)
1433  MVT::MvTimesMatAddMv (-SCT::one(), *X[k], *C[k], SCT::one(), *Residual);
1434  const magnitude_type residErr = frobeniusNorm (*Residual);
1435  sout << "-- ||S - Q(:, 1:" << reportedRank << ")*B(1:"
1436  << reportedRank << ", :) - X1*C1 - X2*C2||_F = "
1437  << residErr << endl;
1438  if (residErr > ATOL * TOL)
1439  {
1440  sout << " *** Error: ||S - Q(:, 1:" << reportedRank
1441  << ")*B(1:" << reportedRank << ", :) "
1442  << "- X1*C1 - X2*C2||_F = " << residErr
1443  << " > ATOL*TOL = " << (ATOL*TOL) << "." << endl;
1444  numerr++;
1445  }
1446  // Verify that Q(1:reportedRank) is orthogonal to X[k], for
1447  // all k. This test only makes sense if reportedRank > 0.
1448  if (reportedRank == 0)
1449  {
1450  sout << "-- Reported rank of Q is zero: skipping Q, X[k] "
1451  "orthogonality test." << endl;
1452  }
1453  else
1454  {
1455  for (int k = 0; k < num_X; ++k)
1456  {
1457  // Q should be orthogonal to X[k], for all k.
1458  const magnitude_type projErr = OM->orthogError(*X[k], *Q_left);
1459  sout << "-- ||<Q(1:" << reportedRank << "), X[" << k
1460  << "]>||_F = " << projErr << endl;
1461  if (projErr > ATOL*TOL)
1462  {
1463  sout << " *** Error: ||<Q(1:" << reportedRank << "), X["
1464  << k << "]>||_F = " << projErr << " > ATOL*TOL = "
1465  << (ATOL*TOL) << "." << endl;
1466  numerr++;
1467  }
1468  }
1469  }
1470  } catch (Belos::OrthoError& e) {
1471  sout << " *** Error: The OrthoManager subclass instance threw "
1472  "an exception: " << e.what() << endl;
1473  numerr++;
1474  }
1475 
1476  // Print out the collected diagnostic messages, which possibly
1477  // include error messages.
1478  const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
1479  MyOM->stream(type) << sout.str();
1480  MyOM->stream(type) << endl;
1481 
1482  return numerr;
1483  }
1484 
1485 
1489  static int
1490  testProjectNew (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1491  const Teuchos::RCP< const MV >& S,
1492  const Teuchos::RCP< const MV >& X1,
1493  const Teuchos::RCP< const MV >& X2,
1494  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
1495  {
1496  using Teuchos::Array;
1497  using Teuchos::null;
1498  using Teuchos::RCP;
1499  using Teuchos::rcp;
1500  using Teuchos::tuple;
1501 
1502  // We collect all the output in this string wrapper, and print
1503  // it at the end.
1504  std::ostringstream sout;
1505  // Total number of failed tests in this call of this routine.
1506  int numerr = 0;
1507 
1508  const int numRows = MVT::GetGlobalLength(*S);
1509  const int numCols = MVT::GetNumberVecs(*S);
1510  const int sizeS = MVT::GetNumberVecs(*S);
1511 
1512  // Relative tolerance against which all tests are performed.
1513  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1514  // The following bounds hold for all $m \times n$ matrices $A$:
1515  // \[
1516  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1517  // \]
1518  // where $r$ is the (column) rank of $A$. We bound this above
1519  // by the number of columns in $A$.
1520  //
1521  // Since we are measuring both of these quantitites in the
1522  // Frobenius norm instead, we scale all error tests by
1523  // $\sqrt{n}$.
1524  //
1525  // A good heuristic is to scale the bound by the square root
1526  // of the number of floating-point operations. One could
1527  // perhaps support this theoretically, since we are using
1528  // uniform random test problems.
1529  const magnitude_type fudgeFactor =
1530  SMT::squareroot(magnitude_type(numRows) *
1531  magnitude_type(numCols) *
1532  magnitude_type(numCols));
1533  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1534  SMT::squareroot(magnitude_type(numCols));
1535 
1536  // Absolute tolerance scaling: the Frobenius norm of the test
1537  // matrix S. TOL*ATOL is the absolute tolerance for the
1538  // residual $\|A - Q*B\|_F$.
1539  const magnitude_type ATOL = frobeniusNorm (*S);
1540 
1541  sout << "The test matrix S has Frobenius norm " << ATOL
1542  << ", and the relative error tolerance is TOL = "
1543  << TOL << "." << endl;
1544 
1545  // Make some copies of S, X1, and X2. The OrthoManager's
1546  // project() method shouldn't modify X1 or X2, but this is a a
1547  // test and we don't know that it doesn't!
1548  RCP< MV > S_copy = MVT::CloneCopy(*S);
1549  RCP< MV > Residual = MVT::CloneCopy(*S);
1550  const int num_X = 2;
1551  Array< RCP< const MV > > X (num_X);
1552  X[0] = MVT::CloneCopy(*X1);
1553  X[1] = MVT::CloneCopy(*X2);
1554 
1555  // Array of coefficients matrices from the projection.
1556  // For our first test, we allocate each of these matrices
1557  // with the proper dimensions.
1558  Array< RCP< mat_type > > C (num_X);
1559  for (int k = 0; k < num_X; ++k)
1560  {
1561  C[k] = rcp (new mat_type (MVT::GetNumberVecs(*X[k]), sizeS));
1562  Teuchos::randomSyncedMatrix(*C[k]); // will be overwritten
1563  }
1564  try {
1565  // Compute the projection: S_copy := (I - X X^*) S
1566  OM->project(*S_copy, C, X);
1567 
1568  // Compute the residual: if successful, S = S_copy + X (X^*
1569  // S =: C) in exact arithmetic. So, the residual is
1570  // S - S_copy - X1 C1 - X2 C2.
1571  //
1572  // Residual := S - S_copy
1573  MVT::MvAddMv (SCT::one(), *S, -SCT::one(), *S_copy, *Residual);
1574  // Residual := Residual - X[k]*C[k]
1575  for (int k = 0; k < num_X; ++k)
1576  MVT::MvTimesMatAddMv (-SCT::one(), *X[k], *C[k], SCT::one(), *Residual);
1577  magnitude_type residErr = frobeniusNorm (*Residual);
1578  sout << " ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr;
1579  if (residErr > ATOL * TOL)
1580  {
1581  sout << " *** Error: ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr
1582  << " > ATOL*TOL = " << (ATOL*TOL) << ".";
1583  numerr++;
1584  }
1585  for (int k = 0; k < num_X; ++k)
1586  {
1587  // S_copy should be orthogonal to X[k] now.
1588  const magnitude_type projErr = OM->orthogError(*X[k], *S_copy);
1589  if (projErr > TOL)
1590  {
1591  sout << " *** Error: S is not orthogonal to X[" << k
1592  << "] by a factor of " << projErr << " > TOL = "
1593  << TOL << ".";
1594  numerr++;
1595  }
1596  }
1597  } catch (Belos::OrthoError& e) {
1598  sout << " *** Error: The OrthoManager subclass instance threw "
1599  "an exception: " << e.what() << endl;
1600  numerr++;
1601  }
1602 
1603  // Print out the collected diagnostic messages, which possibly
1604  // include error messages.
1605  const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
1606  MyOM->stream(type) << sout.str();
1607  MyOM->stream(type) << endl;
1608 
1609  return numerr;
1610  }
1611 
1612  static int
1613  testProject (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1614  const Teuchos::RCP< const MV >& S,
1615  const Teuchos::RCP< const MV >& X1,
1616  const Teuchos::RCP< const MV >& X2,
1617  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
1618  {
1619  return testProjectNew (OM, S, X1, X2, MyOM);
1620  }
1621 
1625  static int
1626  testProjectOld (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1627  const Teuchos::RCP< const MV >& S,
1628  const Teuchos::RCP< const MV >& X1,
1629  const Teuchos::RCP< const MV >& X2,
1630  const Teuchos::RCP< Belos::OutputManager< Scalar > >& MyOM)
1631  {
1632  using Teuchos::Array;
1633  using Teuchos::null;
1634  using Teuchos::RCP;
1635  using Teuchos::rcp;
1636  using Teuchos::tuple;
1637 
1638  const scalar_type ONE = SCT::one();
1639  // We collect all the output in this string wrapper, and print
1640  // it at the end.
1641  std::ostringstream sout;
1642  // Total number of failed tests in this call of this routine.
1643  int numerr = 0;
1644 
1645  const int numRows = MVT::GetGlobalLength(*S);
1646  const int numCols = MVT::GetNumberVecs(*S);
1647  const int sizeS = MVT::GetNumberVecs(*S);
1648  const int sizeX1 = MVT::GetNumberVecs(*X1);
1649  const int sizeX2 = MVT::GetNumberVecs(*X2);
1650 
1651  // Relative tolerance against which all tests are performed.
1652  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1653  // The following bounds hold for all $m \times n$ matrices $A$:
1654  // \[
1655  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1656  // \]
1657  // where $r$ is the (column) rank of $A$. We bound this above
1658  // by the number of columns in $A$.
1659  //
1660  // Since we are measuring both of these quantitites in the
1661  // Frobenius norm instead, we scale all error tests by
1662  // $\sqrt{n}$.
1663  //
1664  // A good heuristic is to scale the bound by the square root
1665  // of the number of floating-point operations. One could
1666  // perhaps support this theoretically, since we are using
1667  // uniform random test problems.
1668  const magnitude_type fudgeFactor =
1669  SMT::squareroot(magnitude_type(numRows) *
1670  magnitude_type(numCols) *
1671  magnitude_type(numCols));
1672  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1673  SMT::squareroot(magnitude_type(numCols));
1674 
1675  // Absolute tolerance scaling: the Frobenius norm of the test
1676  // matrix S. TOL*ATOL is the absolute tolerance for the
1677  // residual $\|A - Q*B\|_F$.
1678  const magnitude_type ATOL = frobeniusNorm (*S);
1679 
1680  sout << "The test matrix S has Frobenius norm " << ATOL
1681  << ", and the relative error tolerance is TOL = "
1682  << TOL << "." << endl;
1683 
1684 
1685  //
1686  // Output tests:
1687  // <S_out,X1> = 0
1688  // <S_out,X2> = 0
1689  // S_in = S_out + X1 C1 + X2 C2
1690  //
1691  // We will loop over an integer specifying the test combinations.
1692  // The bit pattern for the different tests is listed in parentheses.
1693  //
1694  // For the projectors, test the following combinations:
1695  // none (00)
1696  // P_X1 (01)
1697  // P_X2 (10)
1698  // P_X1 P_X2 (11)
1699  // P_X2 P_X1 (11)
1700  // The latter two should be tested to give the same result.
1701  //
1702  // For each of these, we should test with C1 and C2:
1703  //
1704  // if hasM:
1705  // with and without MX1 (1--)
1706  // with and without MX2 (1---)
1707  // with and without MS (1----)
1708  //
1709  // As hasM controls the upper level bits, we need only run test
1710  // cases 0-3 if hasM==false. Otherwise, we run test cases 0-31.
1711  //
1712 
1713  int numtests = 8;
1714 
1715  // test ortho error before orthonormalizing
1716  if (X1 != null) {
1717  magnitude_type err = OM->orthogError(*S,*X1);
1718  sout << " || <S,X1> || before : " << err << endl;
1719  }
1720  if (X2 != null) {
1721  magnitude_type err = OM->orthogError(*S,*X2);
1722  sout << " || <S,X2> || before : " << err << endl;
1723  }
1724 
1725  for (int t = 0; t < numtests; ++t)
1726  {
1727  Array< RCP< const MV > > theX;
1728  Array< RCP< mat_type > > C;
1729  if ( (t % 3) == 0 ) {
1730  // neither X1 nor X2
1731  // C and theX are already empty
1732  }
1733  else if ( (t % 3) == 1 ) {
1734  // X1
1735  theX = tuple(X1);
1736  C = tuple( rcp(new mat_type(sizeX1,sizeS)) );
1737  }
1738  else if ( (t % 3) == 2 ) {
1739  // X2
1740  theX = tuple(X2);
1741  C = tuple( rcp(new mat_type(sizeX2,sizeS)) );
1742  }
1743  else {
1744  // X1 and X2, and the reverse.
1745  theX = tuple(X1,X2);
1746  C = tuple( rcp(new mat_type(sizeX1,sizeS)),
1747  rcp(new mat_type(sizeX2,sizeS)) );
1748  }
1749 
1750  try {
1751  // call routine
1752  // if (t && 3) == 3, {
1753  // call with reversed input: X2 X1
1754  // }
1755  // test all outputs for correctness
1756  // test all outputs for equivalence
1757 
1758  // here is where the outputs go
1759  Array< RCP< MV > > S_outs;
1760  Array< Array< RCP< mat_type > > > C_outs;
1761  RCP< MV > Scopy;
1762 
1763  // copies of S,MS
1764  Scopy = MVT::CloneCopy(*S);
1765  // randomize this data, it should be overwritten
1766  for (size_type i = 0; i < C.size(); ++i) {
1767  Teuchos::randomSyncedMatrix(*C[i]);
1768  }
1769  // Run test.
1770  // Note that Anasazi and Belos differ, among other places,
1771  // in the order of arguments to project().
1772  OM->project(*Scopy,C,theX);
1773  // we allocate S and MS for each test, so we can save these as views
1774  // however, save copies of the C
1775  S_outs.push_back( Scopy );
1776  C_outs.push_back( Array< RCP< mat_type > >(0) );
1777  if (C.size() > 0) {
1778  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
1779  }
1780  if (C.size() > 1) {
1781  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
1782  }
1783 
1784  // do we run the reversed input?
1785  if ( (t % 3) == 3 ) {
1786  // copies of S,MS
1787  Scopy = MVT::CloneCopy(*S);
1788  // randomize this data, it should be overwritten
1789  for (size_type i = 0; i < C.size(); ++i) {
1790  Teuchos::randomSyncedMatrix(*C[i]);
1791  }
1792  // flip the inputs
1793  theX = tuple( theX[1], theX[0] );
1794  // Run test.
1795  // Note that Anasazi and Belos differ, among other places,
1796  // in the order of arguments to project().
1797  OM->project(*Scopy,C,theX);
1798  // we allocate S and MS for each test, so we can save these as views
1799  // however, save copies of the C
1800  S_outs.push_back( Scopy );
1801  // we are in a special case: P_X1 and P_X2, so we know we applied
1802  // two projectors, and therefore have two C[i]
1803  C_outs.push_back( Array<RCP<mat_type > >() );
1804  // reverse the Cs to compensate for the reverse projectors
1805  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
1806  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
1807  // flip the inputs back
1808  theX = tuple( theX[1], theX[0] );
1809  }
1810 
1811  // test all outputs for correctness
1812  for (size_type o = 0; o < S_outs.size(); ++o) {
1813  // S_in = X1*C1 + C2*C2 + S_out
1814  {
1815  RCP<MV> tmp = MVT::CloneCopy(*S_outs[o]);
1816  if (C_outs[o].size() > 0) {
1817  MVT::MvTimesMatAddMv(ONE,*X1,*C_outs[o][0],ONE,*tmp);
1818  if (C_outs[o].size() > 1) {
1819  MVT::MvTimesMatAddMv(ONE,*X2,*C_outs[o][1],ONE,*tmp);
1820  }
1821  }
1822  magnitude_type err = MVDiff(*tmp,*S);
1823  if (err > ATOL*TOL) {
1824  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1825  numerr++;
1826  }
1827  sout << " " << t << "|| S_in - X1*C1 - X2*C2 - S_out || : " << err << endl;
1828  }
1829  // <X1,S> == 0
1830  if (theX.size() > 0 && theX[0] != null) {
1831  magnitude_type err = OM->orthogError(*theX[0],*S_outs[o]);
1832  if (err > TOL) {
1833  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1834  numerr++;
1835  }
1836  sout << " " << t << "|| <X[0],S> || after : " << err << endl;
1837  }
1838  // <X2,S> == 0
1839  if (theX.size() > 1 && theX[1] != null) {
1840  magnitude_type err = OM->orthogError(*theX[1],*S_outs[o]);
1841  if (err > TOL) {
1842  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1843  numerr++;
1844  }
1845  sout << " " << t << "|| <X[1],S> || after : " << err << endl;
1846  }
1847  }
1848 
1849  // test all outputs for equivalence
1850  // check all combinations:
1851  // output 0 == output 1
1852  // output 0 == output 2
1853  // output 1 == output 2
1854  for (size_type o1=0; o1<S_outs.size(); o1++) {
1855  for (size_type o2=o1+1; o2<S_outs.size(); o2++) {
1856  // don't need to check MS_outs because we check
1857  // S_outs and MS_outs = M*S_outs
1858  // don't need to check C_outs either
1859  //
1860  // check that S_outs[o1] == S_outs[o2]
1861  magnitude_type err = MVDiff(*S_outs[o1],*S_outs[o2]);
1862  if (err > TOL) {
1863  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1864  numerr++;
1865  }
1866  }
1867  }
1868 
1869  }
1870  catch (Belos::OrthoError& e) {
1871  sout << " ------------------------------------------- project() threw exception" << endl;
1872  sout << " Error: " << e.what() << endl;
1873  numerr++;
1874  }
1875  } // test for
1876 
1877  MsgType type = Debug;
1878  if (numerr>0) type = Errors;
1879  MyOM->stream(type) << sout.str();
1880  MyOM->stream(type) << endl;
1881 
1882  return numerr;
1883  }
1884 
1885 
1886  };
1887 
1888 
1889 
1890  } // namespace Test
1891 } // namespace Belos
1892 
1893 
Belos::OutputManager
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
Definition: BelosIteration.hpp:32
Teuchos_SerialDenseHelpers.hpp
BelosOutputManager.hpp
Class which manages the output and verbosity of the Belos solvers.
Belos::Test::OrthoManagerBenchmarker::baseline
static void baseline(const Teuchos::RCP< const MV > &X, const int numCols, const int numBlocks, const int numTrials)
Establish baseline run time for OrthoManager benchmark.
Definition: BelosOrthoManagerTest.hpp:52
Teuchos::ScalarTraits::eps
static magnitudeType eps()
Belos::MultiVecTraits::MvRandom
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
Definition: BelosMultiVecTraits.hpp:300
Belos::OutOfPlaceNormalizerMixin
Mixin for out-of-place orthogonalization.
Definition: BelosTsqrOrthoManager.hpp:47
Belos::MsgType
MsgType
Available message types recognized by the linear solvers.
Definition: BelosTypes.hpp:222
Teuchos::ScalarTraits::squareroot
static T squareroot(T x)
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Belos::MultiVecTraits::CloneView
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
Definition: BelosMultiVecTraits.hpp:160
BelosMultiVecTraits.hpp
Declaration of basic traits for the multivector type.
Belos::Test::OrthoManagerTester::mat_type
Teuchos::SerialDenseMatrix< int, scalar_type > mat_type
Definition: BelosOrthoManagerTest.hpp:303
Belos::Test::OrthoManagerTester::MVT
MultiVecTraits< scalar_type, MV > MVT
Definition: BelosOrthoManagerTest.hpp:302
Belos::MultiVecTraits::MvInit
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
Definition: BelosMultiVecTraits.hpp:305
Belos::MultiVecTraits::MvTransMv
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
Definition: BelosMultiVecTraits.hpp:243
Belos::Test::OrthoManagerTester::SCT
Teuchos::ScalarTraits< scalar_type > SCT
Definition: BelosOrthoManagerTest.hpp:299
Belos::MultiVecTraits::GetNumberVecs
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Definition: BelosMultiVecTraits.hpp:186
Belos::MultiVecTraits::Assign
static void Assign(const MV &A, MV &mv)
mv := A
Definition: BelosMultiVecTraits.hpp:295
Belos::Test::OrthoManagerBenchmarker::benchmark
static void benchmark(const Teuchos::RCP< OrthoManager< Scalar, MV > > &orthoMan, const std::string &orthoManName, const std::string &normalization, const Teuchos::RCP< const MV > &X, const int numCols, const int numBlocks, const int numTrials, const Teuchos::RCP< OutputManager< Scalar > > &outMan, std::ostream &resultStream, const bool displayResultsCompactly=false)
Benchmark the given orthogonalization manager.
Definition: BelosOrthoManagerTest.hpp:122
Belos::MultiVecTraits::MvAddMv
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Definition: BelosMultiVecTraits.hpp:228
Belos::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: BelosMultiVecTraits.hpp:97
Belos::MultiVecTraits::SetBlock
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
Definition: BelosMultiVecTraits.hpp:274
Belos::MultiVecTraits::CloneViewNonConst
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
Definition: BelosMultiVecTraits.hpp:141
Belos::MultiVecTraits::Clone
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Definition: BelosMultiVecTraits.hpp:106
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Belos::OrthoError
Exception thrown to signal error in an orthogonalization manager method.
Definition: BelosOrthoManager.hpp:51
Belos::Test::OrthoManagerTester::runTests
static int runTests(const Teuchos::RCP< OrthoManager< Scalar, MV > > &OM, const bool isRankRevealing, const Teuchos::RCP< MV > &S, const int sizeX1, const int sizeX2, const Teuchos::RCP< OutputManager< Scalar > > &MyOM)
Run all the tests.
Definition: BelosOrthoManagerTest.hpp:322
Belos::MultiVecTraits::CloneCopy
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Definition: BelosMultiVecTraits.hpp:113
Belos::MultiVecTraits::MvTimesMatAddMv
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Definition: BelosMultiVecTraits.hpp:221
Belos::MultiVecTraits::HasConstantStride
static bool HasConstantStride(const MV &mv)
Whether the given multivector mv has constant stride.
Definition: BelosMultiVecTraits.hpp:211
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Belos::Test::OrthoManagerTester::SMT
Teuchos::ScalarTraits< magnitude_type > SMT
Definition: BelosOrthoManagerTest.hpp:301
Belos::MultiVecTraits::GetGlobalLength
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Definition: BelosMultiVecTraits.hpp:182
Belos::MultiVecTraits::MvNorm
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
Definition: BelosMultiVecTraits.hpp:262
Teuchos_TimeMonitor.hpp
BelosConfigDefs.hpp
Belos header file which uses auto-configuration information to include necessary C++ headers...
Belos::Test::OrthoManagerTester
Wrapper around OrthoManager test functionality.
Definition: BelosOrthoManagerTest.hpp:293
Teuchos::RCP::is_null
bool is_null() const
Generated on Fri Dec 20 2024 09:24:49 for Belos by doxygen 1.8.5