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

Generated on Thu Mar 28 2024 09:27:41 for Belos by doxygen 1.8.5