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 
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 
564  {
565  //
566  // Test project() on a random multivector S, by projecting S
567  // against various combinations of X1 and X2.
568  //
569  MVT::MvRandom(*S);
570 
571  debugOut << "Testing project() by projecting a random multivector S "
572  "against various combinations of X1 and X2 " << endl;
573  const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
574  numFailed += thisNumFailed;
575  if (thisNumFailed > 0)
576  debugOut << " *** " << thisNumFailed
577  << (thisNumFailed > 1 ? " tests" : " test")
578  << " failed." << endl;
579  }
580 
581  if (isRankRevealing)
582  {
583  // run a X1,Y2 range multivector against P_{X1,X1} P_{Y2,Y2}
584  // note, this is allowed under the restrictions on project(),
585  // because <X1,Y2> = 0
586  // also, <Y2,Y2> = I, but <X1,X1> != I, so biOrtho must be set to false
587  // it should require randomization, as
588  // P_{X1,X1} P_{Y2,Y2} (X1*C1 + Y2*C2) = P_{X1,X1} X1*C1 = 0
589  mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
590  Teuchos::randomSyncedMatrix(C1);
591  Teuchos::randomSyncedMatrix(C2);
592  // S := X1*C1
593  MVT::MvTimesMatAddMv(ONE,*X1,C1,ZERO,*S);
594  // S := S + X2*C2
595  MVT::MvTimesMatAddMv(ONE,*X2,C2,ONE,*S);
596 
597  debugOut << "Testing project() by projecting [X1 X2]-range multivector "
598  "against P_X1 P_X2 " << endl;
599  const int thisNumFailed = testProject(OM,S,X1,X2,MyOM);
600  numFailed += thisNumFailed;
601  if (thisNumFailed > 0)
602  debugOut << " *** " << thisNumFailed
603  << (thisNumFailed > 1 ? " tests" : " test")
604  << " failed." << endl;
605  }
606 
607  // This test is only distinct from the rank-1 multivector test
608  // (below) if S has at least 3 columns.
609  if (isRankRevealing && sizeS > 2)
610  {
611  MVT::MvRandom(*S);
612  RCP<MV> mid = MVT::Clone(*S,1);
613  mat_type c(sizeS,1);
614  MVT::MvTimesMatAddMv(ONE,*S,c,ZERO,*mid);
615  std::vector<int> ind(1);
616  ind[0] = sizeS-1;
617  MVT::SetBlock(*mid,ind,*S);
618 
619  debugOut << "Testing normalize() on a rank-deficient multivector " << endl;
620  const int thisNumFailed = testNormalize(OM,S,MyOM);
621  numFailed += thisNumFailed;
622  if (thisNumFailed > 0)
623  debugOut << " *** " << thisNumFailed
624  << (thisNumFailed > 1 ? " tests" : " test")
625  << " failed." << endl;
626  }
627 
628  // This test will only exercise rank deficiency if S has at least 2
629  // columns.
630  if (isRankRevealing && sizeS > 1)
631  {
632  // rank-1
633  RCP<MV> one = MVT::Clone(*S,1);
634  MVT::MvRandom(*one);
635  mat_type scaleS(sizeS,1);
636  Teuchos::randomSyncedMatrix(scaleS);
637  // put multiple of column 0 in columns 0:sizeS-1
638  for (int i=0; i<sizeS; i++)
639  {
640  std::vector<int> ind(1);
641  ind[0] = i;
642  RCP<MV> Si = MVT::CloneViewNonConst(*S,ind);
643  MVT::MvAddMv(scaleS[i],*one,ZERO,*one,*Si);
644  }
645  debugOut << "Testing normalize() on a rank-1 multivector " << endl;
646  const int thisNumFailed = testNormalize(OM,S,MyOM);
647  numFailed += thisNumFailed;
648  if (thisNumFailed > 0)
649  debugOut << " *** " << thisNumFailed
650  << (thisNumFailed > 1 ? " tests" : " test")
651  << " failed." << endl;
652  }
653 
654  {
655  std::vector<int> ind(1);
656  MVT::MvRandom(*S);
657 
658  debugOut << "Testing projectAndNormalize() on a random multivector " << endl;
659  const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
660  numFailed += thisNumFailed;
661  if (thisNumFailed > 0)
662  debugOut << " *** " << thisNumFailed
663  << (thisNumFailed > 1 ? " tests" : " test")
664  << " failed." << endl;
665  }
666 
667  if (isRankRevealing)
668  {
669  // run a X1,X2 range multivector against P_X1 P_X2
670  // this is allowed as <X1,X2> == 0
671  // it should require randomization, as
672  // P_X1 P_X2 (X1*C1 + X2*C2) = P_X1 X1*C1 = 0
673  // and
674  // P_X2 P_X1 (X2*C2 + X1*C1) = P_X2 X2*C2 = 0
675  mat_type C1(sizeX1,sizeS), C2(sizeX2,sizeS);
676  Teuchos::randomSyncedMatrix(C1);
677  Teuchos::randomSyncedMatrix(C2);
678  MVT::MvTimesMatAddMv(ONE,*X1,C1,ZERO,*S);
679  MVT::MvTimesMatAddMv(ONE,*X2,C2,ONE,*S);
680 
681  debugOut << "Testing projectAndNormalize() by projecting [X1 X2]-range "
682  "multivector against P_X1 P_X2 " << 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 is only distinct from the rank-1 multivector test
692  // (below) if S has at least 3 columns.
693  if (isRankRevealing && sizeS > 2)
694  {
695  MVT::MvRandom(*S);
696  RCP<MV> mid = MVT::Clone(*S,1);
697  mat_type c(sizeS,1);
698  MVT::MvTimesMatAddMv(ONE,*S,c,ZERO,*mid);
699  std::vector<int> ind(1);
700  ind[0] = sizeS-1;
701  MVT::SetBlock(*mid,ind,*S);
702 
703  debugOut << "Testing projectAndNormalize() on a rank-deficient "
704  "multivector " << endl;
705  const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
706  numFailed += thisNumFailed;
707  if (thisNumFailed > 0)
708  debugOut << " *** " << thisNumFailed
709  << (thisNumFailed > 1 ? " tests" : " test")
710  << " failed." << endl;
711  }
712 
713  // This test will only exercise rank deficiency if S has at least 2
714  // columns.
715  if (isRankRevealing && sizeS > 1)
716  {
717  // rank-1
718  RCP<MV> one = MVT::Clone(*S,1);
719  MVT::MvRandom(*one);
720  mat_type scaleS(sizeS,1);
721  Teuchos::randomSyncedMatrix(scaleS);
722  // Put a multiple of column 0 in columns 0:sizeS-1.
723  for (int i=0; i<sizeS; i++)
724  {
725  std::vector<int> ind(1);
726  ind[0] = i;
727  RCP<MV> Si = MVT::CloneViewNonConst(*S,ind);
728  MVT::MvAddMv(scaleS[i],*one,ZERO,*one,*Si);
729  }
730  debugOut << "Testing projectAndNormalize() on a rank-1 multivector " << endl;
731  bool constantStride = true;
732  if (! MVT::HasConstantStride(*S)) {
733  debugOut << "-- S does not have constant stride" << endl;
734  constantStride = false;
735  }
736  if (! MVT::HasConstantStride(*X1)) {
737  debugOut << "-- X1 does not have constant stride" << endl;
738  constantStride = false;
739  }
740  if (! MVT::HasConstantStride(*X2)) {
741  debugOut << "-- X2 does not have constant stride" << endl;
742  constantStride = false;
743  }
744  if (! constantStride) {
745  debugOut << "-- Skipping this test, since TSQR does not work on "
746  "multivectors with nonconstant stride" << endl;
747  }
748  else {
749  const int thisNumFailed = testProjectAndNormalize(OM,S,X1,X2,MyOM);
750  numFailed += thisNumFailed;
751  if (thisNumFailed > 0) {
752  debugOut << " *** " << thisNumFailed
753  << (thisNumFailed > 1 ? " tests" : " test")
754  << " failed." << endl;
755  }
756  }
757  }
758 
759  if (numFailed != 0) {
760  MyOM->stream(Errors) << numFailed << " total test failures." << endl;
761  }
762  return numFailed;
763  }
764 
765  private:
766 
771  static magnitude_type
772  MVDiff (const MV& X, const MV& Y)
773  {
774  using Teuchos::RCP;
775 
776  const scalar_type ONE = SCT::one();
777  const int numCols = MVT::GetNumberVecs(X);
779  std::logic_error,
780  "MVDiff: X and Y should have the same number of columns."
781  " X has " << numCols << " column(s) and Y has "
782  << MVT::GetNumberVecs(Y) << " columns." );
783  // Resid := X
784  RCP< MV > Resid = MVT::CloneCopy(X);
785  // Resid := Resid - Y
786  MVT::MvAddMv (-ONE, Y, ONE, *Resid, *Resid);
787 
788  return frobeniusNorm (*Resid);
789  }
790 
791 
795  static magnitude_type
796  frobeniusNorm (const MV& X)
797  {
798  const scalar_type ONE = SCT::one();
799  const int numCols = MVT::GetNumberVecs(X);
800  mat_type C (numCols, numCols);
801 
802  // $C := X^* X$
803  MVT::MvTransMv (ONE, X, X, C);
804 
805  magnitude_type err (0);
806  for (int i = 0; i < numCols; ++i)
807  err += SCT::magnitude (C(i,i));
808 
809  return SCT::magnitude (SCT::squareroot (err));
810  }
811 
812 
813  static int
814  testProjectAndNormalize (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
815  const Teuchos::RCP< const MV >& S,
816  const Teuchos::RCP< const MV >& X1,
817  const Teuchos::RCP< const MV >& X2,
819  {
820  return testProjectAndNormalizeNew (OM, S, X1, X2, MyOM);
821  }
822 
827  static int
828  testProjectAndNormalizeOld (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > >& OM,
829  const Teuchos::RCP< const MV >& S,
830  const Teuchos::RCP< const MV >& X1,
831  const Teuchos::RCP< const MV >& X2,
833  {
834  using Teuchos::Array;
835  using Teuchos::null;
836  using Teuchos::RCP;
837  using Teuchos::rcp;
838  using Teuchos::tuple;
839 
840  const scalar_type ONE = SCT::one();
841  const magnitude_type ZERO = SCT::magnitude(SCT::zero());
842 
843  // Relative tolerance against which all tests are performed.
844  const magnitude_type TOL = 1.0e-12;
845  // Absolute tolerance constant
846  const magnitude_type ATOL = 10;
847 
848  const int sizeS = MVT::GetNumberVecs(*S);
849  const int sizeX1 = MVT::GetNumberVecs(*X1);
850  const int sizeX2 = MVT::GetNumberVecs(*X2);
851  int numerr = 0;
852  std::ostringstream sout;
853 
854  //
855  // output tests:
856  // <S_out,S_out> = I
857  // <S_out,X1> = 0
858  // <S_out,X2> = 0
859  // S_in = S_out B + X1 C1 + X2 C2
860  //
861  // we will loop over an integer specifying the test combinations
862  // the bit pattern for the different tests is listed in parenthesis
863  //
864  // for the projectors, test the following combinations:
865  // none (00)
866  // P_X1 (01)
867  // P_X2 (10)
868  // P_X1 P_X2 (11)
869  // P_X2 P_X1 (11)
870  // the latter two should be tested to give the same answer
871  //
872  // for each of these, we should test with C1, C2 and B
873  //
874  // if hasM:
875  // with and without MX1 (1--)
876  // with and without MX2 (1---)
877  // with and without MS (1----)
878  //
879  // as hasM controls the upper level bits, we need only run test cases 0-3 if hasM==false
880  // otherwise, we run test cases 0-31
881  //
882 
883  int numtests = 4;
884 
885  // test ortho error before orthonormalizing
886  if (X1 != null) {
887  magnitude_type err = OM->orthogError(*S,*X1);
888  sout << " || <S,X1> || before : " << err << endl;
889  }
890  if (X2 != null) {
891  magnitude_type err = OM->orthogError(*S,*X2);
892  sout << " || <S,X2> || before : " << err << endl;
893  }
894 
895  for (int t=0; t<numtests; t++) {
896 
897  Array< RCP< const MV > > theX;
898  RCP<mat_type > B = rcp( new mat_type(sizeS,sizeS) );
899  Array<RCP<mat_type > > C;
900  if ( (t % 3) == 0 ) {
901  // neither <X1,Y1> nor <X2,Y2>
902  // C, theX and theY are already empty
903  }
904  else if ( (t % 3) == 1 ) {
905  // X1
906  theX = tuple(X1);
907  C = tuple( rcp(new mat_type(sizeX1,sizeS)) );
908  }
909  else if ( (t % 3) == 2 ) {
910  // X2
911  theX = tuple(X2);
912  C = tuple( rcp(new mat_type(sizeX2,sizeS)) );
913  }
914  else {
915  // X1 and X2, and the reverse.
916  theX = tuple(X1,X2);
917  C = tuple( rcp(new mat_type(sizeX1,sizeS)),
918  rcp(new mat_type(sizeX2,sizeS)) );
919  }
920 
921  // We wrap up all the OrthoManager calls in a try-catch
922  // block, in order to check whether any of the methods throw
923  // an exception. For the tests we perform, every thrown
924  // exception is a failure.
925  try {
926  // call routine
927  // if (t && 3) == 3, {
928  // call with reversed input: X2 X1
929  // }
930  // test all outputs for correctness
931  // test all outputs for equivalence
932 
933  // here is where the outputs go
934  Array<RCP<MV> > S_outs;
935  Array<Array<RCP<mat_type > > > C_outs;
936  Array<RCP<mat_type > > B_outs;
937  RCP<MV> Scopy;
938  Array<int> ret_out;
939 
940  // copies of S,MS
941  Scopy = MVT::CloneCopy(*S);
942  // randomize this data, it should be overwritten
943  Teuchos::randomSyncedMatrix(B);
944  for (size_type i=0; i<C.size(); i++) {
945  Teuchos::randomSyncedmatrix(*C[i]);
946  }
947  // Run test. Since S was specified by the caller and
948  // Scopy is a copy of S, we don't know what rank to expect
949  // here -- though we do require that S have rank at least
950  // one.
951  //
952  // Note that Anasazi and Belos differ, among other places,
953  // in the order of arguments to projectAndNormalize().
954  int ret = OM->projectAndNormalize(*Scopy,C,B,theX);
955  sout << "projectAndNormalize() returned rank " << ret << endl;
956  if (ret == 0) {
957  sout << " *** Error: returned rank is zero, cannot continue tests" << endl;
958  numerr++;
959  break;
960  }
961  ret_out.push_back(ret);
962  // projectAndNormalize() is only required to return a
963  // basis of rank "ret"
964  // this is what we will test:
965  // the first "ret" columns in Scopy
966  // the first "ret" rows in B
967  // save just the parts that we want
968  // we allocate S and MS for each test, so we can save these as views
969  // however, save copies of the C and B
970  if (ret < sizeS) {
971  std::vector<int> ind(ret);
972  for (int i=0; i<ret; i++) {
973  ind[i] = i;
974  }
975  S_outs.push_back( MVT::CloneViewNonConst(*Scopy,ind) );
976  B_outs.push_back( rcp( new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
977  }
978  else {
979  S_outs.push_back( Scopy );
980  B_outs.push_back( rcp( new mat_type(*B) ) );
981  }
982  C_outs.push_back( Array<RCP<mat_type > >(0) );
983  if (C.size() > 0) {
984  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
985  }
986  if (C.size() > 1) {
987  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
988  }
989 
990  // do we run the reversed input?
991  if ( (t % 3) == 3 ) {
992  // copies of S,MS
993  Scopy = MVT::CloneCopy(*S);
994 
995  // Fill the B and C[i] matrices with random data. The
996  // data will be overwritten by projectAndNormalize().
997  // Filling these matrices here is only to catch some
998  // bugs in projectAndNormalize().
999  Teuchos::randomSyncedMatrix(B);
1000  for (size_type i=0; i<C.size(); i++) {
1001  Teuchos::randomSyncedMatrix(*C[i]);
1002  }
1003  // flip the inputs
1004  theX = tuple( theX[1], theX[0] );
1005  // Run test.
1006  // Note that Anasazi and Belos differ, among other places,
1007  // in the order of arguments to projectAndNormalize().
1008  ret = OM->projectAndNormalize(*Scopy,C,B,theX);
1009  sout << "projectAndNormalize() returned rank " << ret << endl;
1010  if (ret == 0) {
1011  sout << " *** Error: returned rank is zero, cannot continue tests" << endl;
1012  numerr++;
1013  break;
1014  }
1015  ret_out.push_back(ret);
1016  // projectAndNormalize() is only required to return a
1017  // basis of rank "ret"
1018  // this is what we will test:
1019  // the first "ret" columns in Scopy
1020  // the first "ret" rows in B
1021  // save just the parts that we want
1022  // we allocate S and MS for each test, so we can save these as views
1023  // however, save copies of the C and B
1024  if (ret < sizeS) {
1025  std::vector<int> ind(ret);
1026  for (int i=0; i<ret; i++) {
1027  ind[i] = i;
1028  }
1029  S_outs.push_back( MVT::CloneViewNonConst(*Scopy,ind) );
1030  B_outs.push_back( rcp( new mat_type(Teuchos::Copy,*B,ret,sizeS) ) );
1031  }
1032  else {
1033  S_outs.push_back( Scopy );
1034  B_outs.push_back( rcp( new mat_type(*B) ) );
1035  }
1036  C_outs.push_back( Array<RCP<mat_type > >() );
1037  // reverse the Cs to compensate for the reverse projectors
1038  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
1039  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
1040  // flip the inputs back
1041  theX = tuple( theX[1], theX[0] );
1042  }
1043 
1044 
1045  // test all outputs for correctness
1046  for (size_type o=0; o<S_outs.size(); o++) {
1047  // S^T M S == I
1048  {
1049  magnitude_type err = OM->orthonormError(*S_outs[o]);
1050  if (err > TOL) {
1051  sout << endl
1052  << " *** Test (number " << (t+1) << " of " << numtests
1053  << " total tests) failed: Tolerance exceeded! Error = "
1054  << err << " > TOL = " << TOL << "."
1055  << endl << endl;
1056  numerr++;
1057  }
1058  sout << " || <S,S> - I || after : " << err << endl;
1059  }
1060  // S_in = X1*C1 + C2*C2 + S_out*B
1061  {
1062  RCP<MV> tmp = MVT::Clone(*S,sizeS);
1063  MVT::MvTimesMatAddMv(ONE,*S_outs[o],*B_outs[o],ZERO,*tmp);
1064  if (C_outs[o].size() > 0) {
1065  MVT::MvTimesMatAddMv(ONE,*X1,*C_outs[o][0],ONE,*tmp);
1066  if (C_outs[o].size() > 1) {
1067  MVT::MvTimesMatAddMv(ONE,*X2,*C_outs[o][1],ONE,*tmp);
1068  }
1069  }
1070  magnitude_type err = MVDiff(*tmp,*S);
1071  if (err > ATOL*TOL) {
1072  sout << endl
1073  << " *** Test (number " << (t+1) << " of " << numtests
1074  << " total tests) failed: Tolerance exceeded! Error = "
1075  << err << " > ATOL*TOL = " << (ATOL*TOL) << "."
1076  << endl << endl;
1077  numerr++;
1078  }
1079  sout << " " << t << "|| S_in - X1*C1 - X2*C2 - S_out*B || : " << err << endl;
1080  }
1081  // <X1,S> == 0
1082  if (theX.size() > 0 && theX[0] != null) {
1083  magnitude_type err = OM->orthogError(*theX[0],*S_outs[o]);
1084  if (err > TOL) {
1085  sout << endl
1086  << " *** Test (number " << (t+1) << " of " << numtests
1087  << " total tests) failed: Tolerance exceeded! Error = "
1088  << err << " > TOL = " << TOL << "."
1089  << endl << endl;
1090  numerr++;
1091  }
1092  sout << " " << t << "|| <X[0],S> || after : " << err << endl;
1093  }
1094  // <X2,S> == 0
1095  if (theX.size() > 1 && theX[1] != null) {
1096  magnitude_type err = OM->orthogError(*theX[1],*S_outs[o]);
1097  if (err > TOL) {
1098  sout << endl
1099  << " *** Test (number " << (t+1) << " of " << numtests
1100  << " total tests) failed: Tolerance exceeded! Error = "
1101  << err << " > TOL = " << TOL << "."
1102  << endl << endl;
1103  numerr++;
1104  }
1105  sout << " " << t << "|| <X[1],S> || after : " << err << endl;
1106  }
1107  }
1108  }
1109  catch (Belos::OrthoError& e) {
1110  sout << " *** Error: OrthoManager threw exception: " << e.what() << endl;
1111  numerr++;
1112  }
1113 
1114  } // test for
1115 
1116  // NOTE (mfh 05 Nov 2010) Since Belos::MsgType is an enum,
1117  // doing bitwise logical computations on Belos::MsgType values
1118  // (such as "Debug | Errors") and passing the result into
1119  // MyOM->stream() confuses the compiler. As a result, we have
1120  // to do some type casts to make it work.
1121  const int msgType = (numerr > 0) ?
1122  (static_cast<int>(Debug) | static_cast<int>(Errors)) :
1123  static_cast<int>(Debug);
1124 
1125  // We report debug-level messages always. We also report
1126  // errors if at least one test failed.
1127  MyOM->stream(static_cast< MsgType >(msgType)) << sout.str() << endl;
1128  return numerr;
1129  }
1130 
1135  static int
1136  testNormalize (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > >& OM,
1137  const Teuchos::RCP< const MV >& S,
1139  {
1140  using Teuchos::Array;
1141  using Teuchos::RCP;
1142  using Teuchos::rcp;
1143  using Teuchos::tuple;
1144 
1145  const scalar_type ONE = SCT::one();
1146  std::ostringstream sout;
1147  // Total number of failed tests in this call of this routine.
1148  int numerr = 0;
1149 
1150  // Relative tolerance against which all tests are performed.
1151  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1152  // The following bounds hold for all $m \times n$ matrices $A$:
1153  // \[
1154  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1155  // \]
1156  // where $r$ is the (column) rank of $A$. We bound this above
1157  // by the number of columns in $A$.
1158  //
1159  // An accurate normalization in the Euclidean norm of a matrix
1160  // $A$ with at least as many rows m as columns n, should
1161  // produce orthogonality $\|Q^* Q - I\|_2$ less than a factor
1162  // of machine precision times a low-order polynomial in m and
1163  // n, and residual $\|A - Q B\|_2$ (where $A = Q B$ is the
1164  // computed normalization) less than that bound times the norm
1165  // of $A$.
1166  //
1167  // Since we are measuring both of these quantitites in the
1168  // Frobenius norm instead, we should scale this bound by
1169  // $\sqrt{n}$.
1170 
1171  const int numRows = MVT::GetGlobalLength(*S);
1172  const int numCols = MVT::GetNumberVecs(*S);
1173  const int sizeS = MVT::GetNumberVecs(*S);
1174 
1175  // A good heuristic is to scale the bound by the square root
1176  // of the number of floating-point operations. One could
1177  // perhaps support this theoretically, since we are using
1178  // uniform random test problems.
1179  const magnitude_type fudgeFactor =
1180  SMT::squareroot(magnitude_type(numRows) *
1181  magnitude_type(numCols) *
1182  magnitude_type(numCols));
1183  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1184  SMT::squareroot(magnitude_type(numCols));
1185 
1186  // Absolute tolerance scaling: the Frobenius norm of the test
1187  // matrix S. TOL*ATOL is the absolute tolerance for the
1188  // residual $\|A - Q*B\|_F$.
1189  const magnitude_type ATOL = frobeniusNorm (*S);
1190 
1191  sout << "The test matrix S has Frobenius norm " << ATOL
1192  << ", and the relative error tolerance is TOL = "
1193  << TOL << "." << endl;
1194 
1195  const int numtests = 1;
1196  for (int t = 0; t < numtests; ++t) {
1197 
1198  try {
1199  // call routine
1200  // test all outputs for correctness
1201 
1202  // S_copy gets a copy of S; we normalize in place, so we
1203  // need a copy to check whether the normalization
1204  // succeeded.
1205  RCP< MV > S_copy = MVT::CloneCopy (*S);
1206 
1207  // Matrix of coefficients from the normalization.
1208  RCP< mat_type > B (new mat_type (sizeS, sizeS));
1209  // The contents of B will be overwritten, but fill with
1210  // random data just to make sure that the normalization
1211  // operated on all the elements of B on which it should
1212  // operate.
1213  Teuchos::randomSyncedMatrix(B);
1214 
1215  const int reportedRank = OM->normalize (*S_copy, B);
1216  sout << "normalize() returned rank " << reportedRank << endl;
1217  if (reportedRank == 0) {
1218  sout << " *** Error: Cannot continue, since normalize() "
1219  "reports that S has rank 0" << endl;
1220  numerr++;
1221  break;
1222  }
1223  //
1224  // We don't know in this routine whether the input
1225  // multivector S has full rank; it is only required to
1226  // have nonzero rank. Thus, we extract the first
1227  // reportedRank columns of S_copy and the first
1228  // reportedRank rows of B, and perform tests on them.
1229  //
1230 
1231  // Construct S_view, a view of the first reportedRank
1232  // columns of S_copy.
1233  std::vector<int> indices (reportedRank);
1234  for (int j = 0; j < reportedRank; ++j)
1235  indices[j] = j;
1236  RCP< MV > S_view = MVT::CloneViewNonConst (*S_copy, indices);
1237  // Construct B_top, a copy of the first reportedRank rows
1238  // of B.
1239  //
1240  // NOTE: We create this as a copy and not a view, because
1241  // otherwise it would not be safe with respect to RCPs.
1242  // This is because mat_type uses raw pointers
1243  // inside, so that a view would become invalid when B
1244  // would fall out of scope.
1245  RCP< mat_type > B_top (new mat_type (Teuchos::Copy, *B, reportedRank, sizeS));
1246 
1247  // Check ||<S_view,S_view> - I||
1248  {
1249  const magnitude_type err = OM->orthonormError(*S_view);
1250  if (err > TOL) {
1251  sout << " *** Error: Tolerance exceeded: err = "
1252  << err << " > TOL = " << TOL << endl;
1253  numerr++;
1254  }
1255  sout << " || <S,S> - I || after : " << err << endl;
1256  }
1257  // Check the residual ||Residual|| = ||S_view * B_top -
1258  // S_orig||, where S_orig is a view of the first
1259  // reportedRank columns of S.
1260  {
1261  // Residual is allocated with reportedRank columns. It
1262  // will contain the result of testing the residual error
1263  // of the normalization (i.e., $\|S - S_in*B\|$). It
1264  // should have the dimensions of S. Its initial value
1265  // is a copy of the first reportedRank columns of S.
1266  RCP< MV > Residual = MVT::CloneCopy (*S);
1267 
1268  // Residual := Residual - S_view * B_view
1269  MVT::MvTimesMatAddMv (-ONE, *S_view, *B_top, ONE, *Residual);
1270 
1271  // Compute ||Residual||
1272  const magnitude_type err = frobeniusNorm (*Residual);
1273  if (err > ATOL*TOL) {
1274  sout << " *** Error: Tolerance exceeded: err = "
1275  << err << " > ATOL*TOL = " << (ATOL*TOL) << endl;
1276  numerr++;
1277  }
1278  sout << " " << t << "|| S - Q*B || : " << err << endl;
1279  }
1280  }
1281  catch (Belos::OrthoError& e) {
1282  sout << " *** Error: the OrthoManager's normalize() method "
1283  "threw an exception: " << e.what() << endl;
1284  numerr++;
1285  }
1286 
1287  } // test for
1288 
1289  const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
1290  MyOM->stream(type) << sout.str();
1291  MyOM->stream(type) << endl;
1292 
1293  return numerr;
1294  }
1295 
1300  static int
1301  testProjectAndNormalizeNew (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1302  const Teuchos::RCP< const MV >& S,
1303  const Teuchos::RCP< const MV >& X1,
1304  const Teuchos::RCP< const MV >& X2,
1306  {
1307  using Teuchos::Array;
1308  using Teuchos::null;
1309  using Teuchos::RCP;
1310  using Teuchos::rcp;
1311  using Teuchos::tuple;
1312 
1313  // We collect all the output in this string wrapper, and print
1314  // it at the end.
1315  std::ostringstream sout;
1316  // Total number of failed tests in this call of this routine.
1317  int numerr = 0;
1318 
1319  const int numRows = MVT::GetGlobalLength(*S);
1320  const int numCols = MVT::GetNumberVecs(*S);
1321  const int sizeS = MVT::GetNumberVecs(*S);
1322 
1323  // Relative tolerance against which all tests are performed.
1324  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1325  // The following bounds hold for all $m \times n$ matrices $A$:
1326  // \[
1327  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1328  // \]
1329  // where $r$ is the (column) rank of $A$. We bound this above
1330  // by the number of columns in $A$.
1331  //
1332  // Since we are measuring both of these quantitites in the
1333  // Frobenius norm instead, we scale all error tests by
1334  // $\sqrt{n}$.
1335  //
1336  // A good heuristic is to scale the bound by the square root
1337  // of the number of floating-point operations. One could
1338  // perhaps support this theoretically, since we are using
1339  // uniform random test problems.
1340  const magnitude_type fudgeFactor =
1341  SMT::squareroot(magnitude_type(numRows) *
1342  magnitude_type(numCols) *
1343  magnitude_type(numCols));
1344  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1345  SMT::squareroot(magnitude_type(numCols));
1346 
1347  // Absolute tolerance scaling: the Frobenius norm of the test
1348  // matrix S. TOL*ATOL is the absolute tolerance for the
1349  // residual $\|A - Q*B\|_F$.
1350  const magnitude_type ATOL = frobeniusNorm (*S);
1351 
1352  sout << "-- The test matrix S has Frobenius norm " << ATOL
1353  << ", and the relative error tolerance is TOL = "
1354  << TOL << "." << endl;
1355 
1356  // Q will contain the result of projectAndNormalize() on S.
1357  RCP< MV > Q = MVT::CloneCopy(*S);
1358  // We use this for collecting the residual error components
1359  RCP< MV > Residual = MVT::CloneCopy(*S);
1360  // Number of elements in the X array of blocks against which
1361  // to project S.
1362  const int num_X = 2;
1363  Array< RCP< const MV > > X (num_X);
1364  X[0] = MVT::CloneCopy(*X1);
1365  X[1] = MVT::CloneCopy(*X2);
1366 
1367  // Coefficients for the normalization
1368  RCP< mat_type > B (new mat_type (sizeS, sizeS));
1369 
1370  // Array of coefficients matrices from the projection.
1371  // For our first test, we allocate each of these matrices
1372  // with the proper dimensions.
1373  Array< RCP< mat_type > > C (num_X);
1374  for (int k = 0; k < num_X; ++k)
1375  {
1376  C[k] = rcp (new mat_type (MVT::GetNumberVecs(*X[k]), sizeS));
1377  Teuchos::randomSyncedMatrix(*C[k]); // will be overwritten
1378  }
1379  try {
1380  // Q*B := (I - X X^*) S
1381  const int reportedRank = OM->projectAndNormalize (*Q, C, B, X);
1382 
1383  // Pick out the first reportedRank columns of Q.
1384  std::vector<int> indices (reportedRank);
1385  for (int j = 0; j < reportedRank; ++j)
1386  indices[j] = j;
1387  RCP< const MV > Q_left = MVT::CloneView (*Q, indices);
1388 
1389  // Test whether the first reportedRank columns of Q are
1390  // orthogonal.
1391  {
1392  const magnitude_type orthoError = OM->orthonormError (*Q_left);
1393  sout << "-- ||Q(1:" << reportedRank << ")^* Q(1:" << reportedRank
1394  << ") - I||_F = " << orthoError << endl;
1395  if (orthoError > TOL)
1396  {
1397  sout << " *** Error: ||Q(1:" << reportedRank << ")^* Q(1:"
1398  << reportedRank << ") - I||_F = " << orthoError
1399  << " > TOL = " << TOL << "." << endl;
1400  numerr++;
1401  }
1402  }
1403 
1404  // Compute the residual: if successful, S = Q*B +
1405  // X (X^* S =: C) in exact arithmetic. So, the residual is
1406  // S - Q*B - X1 C1 - X2 C2.
1407  //
1408  // Residual := S
1409  MVT::MvAddMv (SCT::one(), *S, SCT::zero(), *Residual, *Residual);
1410  {
1411  // Pick out the first reportedRank rows of B. Make a deep
1412  // copy, since mat_type is not safe with respect
1413  // to RCP-based memory management (it uses raw pointers
1414  // inside).
1415  RCP< const mat_type > B_top (new mat_type (Teuchos::Copy, *B, reportedRank, B->numCols()));
1416  // Residual := Residual - Q(:, 1:reportedRank) * B(1:reportedRank, :)
1417  MVT::MvTimesMatAddMv (-SCT::one(), *Q_left, *B_top, SCT::one(), *Residual);
1418  }
1419  // Residual := Residual - X[k]*C[k]
1420  for (int k = 0; k < num_X; ++k)
1421  MVT::MvTimesMatAddMv (-SCT::one(), *X[k], *C[k], SCT::one(), *Residual);
1422  const magnitude_type residErr = frobeniusNorm (*Residual);
1423  sout << "-- ||S - Q(:, 1:" << reportedRank << ")*B(1:"
1424  << reportedRank << ", :) - X1*C1 - X2*C2||_F = "
1425  << residErr << endl;
1426  if (residErr > ATOL * TOL)
1427  {
1428  sout << " *** Error: ||S - Q(:, 1:" << reportedRank
1429  << ")*B(1:" << reportedRank << ", :) "
1430  << "- X1*C1 - X2*C2||_F = " << residErr
1431  << " > ATOL*TOL = " << (ATOL*TOL) << "." << endl;
1432  numerr++;
1433  }
1434  // Verify that Q(1:reportedRank) is orthogonal to X[k], for
1435  // all k. This test only makes sense if reportedRank > 0.
1436  if (reportedRank == 0)
1437  {
1438  sout << "-- Reported rank of Q is zero: skipping Q, X[k] "
1439  "orthogonality test." << endl;
1440  }
1441  else
1442  {
1443  for (int k = 0; k < num_X; ++k)
1444  {
1445  // Q should be orthogonal to X[k], for all k.
1446  const magnitude_type projErr = OM->orthogError(*X[k], *Q_left);
1447  sout << "-- ||<Q(1:" << reportedRank << "), X[" << k
1448  << "]>||_F = " << projErr << endl;
1449  if (projErr > ATOL*TOL)
1450  {
1451  sout << " *** Error: ||<Q(1:" << reportedRank << "), X["
1452  << k << "]>||_F = " << projErr << " > ATOL*TOL = "
1453  << (ATOL*TOL) << "." << endl;
1454  numerr++;
1455  }
1456  }
1457  }
1458  } catch (Belos::OrthoError& e) {
1459  sout << " *** Error: The OrthoManager subclass instance threw "
1460  "an exception: " << e.what() << endl;
1461  numerr++;
1462  }
1463 
1464  // Print out the collected diagnostic messages, which possibly
1465  // include error messages.
1466  const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
1467  MyOM->stream(type) << sout.str();
1468  MyOM->stream(type) << endl;
1469 
1470  return numerr;
1471  }
1472 
1473 
1477  static int
1478  testProjectNew (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1479  const Teuchos::RCP< const MV >& S,
1480  const Teuchos::RCP< const MV >& X1,
1481  const Teuchos::RCP< const MV >& X2,
1483  {
1484  using Teuchos::Array;
1485  using Teuchos::null;
1486  using Teuchos::RCP;
1487  using Teuchos::rcp;
1488  using Teuchos::tuple;
1489 
1490  // We collect all the output in this string wrapper, and print
1491  // it at the end.
1492  std::ostringstream sout;
1493  // Total number of failed tests in this call of this routine.
1494  int numerr = 0;
1495 
1496  const int numRows = MVT::GetGlobalLength(*S);
1497  const int numCols = MVT::GetNumberVecs(*S);
1498  const int sizeS = MVT::GetNumberVecs(*S);
1499 
1500  // Relative tolerance against which all tests are performed.
1501  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1502  // The following bounds hold for all $m \times n$ matrices $A$:
1503  // \[
1504  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1505  // \]
1506  // where $r$ is the (column) rank of $A$. We bound this above
1507  // by the number of columns in $A$.
1508  //
1509  // Since we are measuring both of these quantitites in the
1510  // Frobenius norm instead, we scale all error tests by
1511  // $\sqrt{n}$.
1512  //
1513  // A good heuristic is to scale the bound by the square root
1514  // of the number of floating-point operations. One could
1515  // perhaps support this theoretically, since we are using
1516  // uniform random test problems.
1517  const magnitude_type fudgeFactor =
1518  SMT::squareroot(magnitude_type(numRows) *
1519  magnitude_type(numCols) *
1520  magnitude_type(numCols));
1521  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1522  SMT::squareroot(magnitude_type(numCols));
1523 
1524  // Absolute tolerance scaling: the Frobenius norm of the test
1525  // matrix S. TOL*ATOL is the absolute tolerance for the
1526  // residual $\|A - Q*B\|_F$.
1527  const magnitude_type ATOL = frobeniusNorm (*S);
1528 
1529  sout << "The test matrix S has Frobenius norm " << ATOL
1530  << ", and the relative error tolerance is TOL = "
1531  << TOL << "." << endl;
1532 
1533  // Make some copies of S, X1, and X2. The OrthoManager's
1534  // project() method shouldn't modify X1 or X2, but this is a a
1535  // test and we don't know that it doesn't!
1536  RCP< MV > S_copy = MVT::CloneCopy(*S);
1537  RCP< MV > Residual = MVT::CloneCopy(*S);
1538  const int num_X = 2;
1539  Array< RCP< const MV > > X (num_X);
1540  X[0] = MVT::CloneCopy(*X1);
1541  X[1] = MVT::CloneCopy(*X2);
1542 
1543  // Array of coefficients matrices from the projection.
1544  // For our first test, we allocate each of these matrices
1545  // with the proper dimensions.
1546  Array< RCP< mat_type > > C (num_X);
1547  for (int k = 0; k < num_X; ++k)
1548  {
1549  C[k] = rcp (new mat_type (MVT::GetNumberVecs(*X[k]), sizeS));
1550  Teuchos::randomSyncedMatrix(*C[k]); // will be overwritten
1551  }
1552  try {
1553  // Compute the projection: S_copy := (I - X X^*) S
1554  OM->project(*S_copy, C, X);
1555 
1556  // Compute the residual: if successful, S = S_copy + X (X^*
1557  // S =: C) in exact arithmetic. So, the residual is
1558  // S - S_copy - X1 C1 - X2 C2.
1559  //
1560  // Residual := S - S_copy
1561  MVT::MvAddMv (SCT::one(), *S, -SCT::one(), *S_copy, *Residual);
1562  // Residual := Residual - X[k]*C[k]
1563  for (int k = 0; k < num_X; ++k)
1564  MVT::MvTimesMatAddMv (-SCT::one(), *X[k], *C[k], SCT::one(), *Residual);
1565  magnitude_type residErr = frobeniusNorm (*Residual);
1566  sout << " ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr;
1567  if (residErr > ATOL * TOL)
1568  {
1569  sout << " *** Error: ||S - S_copy - X1*C1 - X2*C2||_F = " << residErr
1570  << " > ATOL*TOL = " << (ATOL*TOL) << ".";
1571  numerr++;
1572  }
1573  for (int k = 0; k < num_X; ++k)
1574  {
1575  // S_copy should be orthogonal to X[k] now.
1576  const magnitude_type projErr = OM->orthogError(*X[k], *S_copy);
1577  if (projErr > TOL)
1578  {
1579  sout << " *** Error: S is not orthogonal to X[" << k
1580  << "] by a factor of " << projErr << " > TOL = "
1581  << TOL << ".";
1582  numerr++;
1583  }
1584  }
1585  } catch (Belos::OrthoError& e) {
1586  sout << " *** Error: The OrthoManager subclass instance threw "
1587  "an exception: " << e.what() << endl;
1588  numerr++;
1589  }
1590 
1591  // Print out the collected diagnostic messages, which possibly
1592  // include error messages.
1593  const MsgType type = (numerr == 0) ? Debug : static_cast<MsgType> (static_cast<int>(Errors) | static_cast<int>(Debug));
1594  MyOM->stream(type) << sout.str();
1595  MyOM->stream(type) << endl;
1596 
1597  return numerr;
1598  }
1599 
1600  static int
1601  testProject (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1602  const Teuchos::RCP< const MV >& S,
1603  const Teuchos::RCP< const MV >& X1,
1604  const Teuchos::RCP< const MV >& X2,
1606  {
1607  return testProjectNew (OM, S, X1, X2, MyOM);
1608  }
1609 
1613  static int
1614  testProjectOld (const Teuchos::RCP< Belos::OrthoManager< Scalar, MV > > OM,
1615  const Teuchos::RCP< const MV >& S,
1616  const Teuchos::RCP< const MV >& X1,
1617  const Teuchos::RCP< const MV >& X2,
1619  {
1620  using Teuchos::Array;
1621  using Teuchos::null;
1622  using Teuchos::RCP;
1623  using Teuchos::rcp;
1624  using Teuchos::tuple;
1625 
1626  const scalar_type ONE = SCT::one();
1627  // We collect all the output in this string wrapper, and print
1628  // it at the end.
1629  std::ostringstream sout;
1630  // Total number of failed tests in this call of this routine.
1631  int numerr = 0;
1632 
1633  const int numRows = MVT::GetGlobalLength(*S);
1634  const int numCols = MVT::GetNumberVecs(*S);
1635  const int sizeS = MVT::GetNumberVecs(*S);
1636  const int sizeX1 = MVT::GetNumberVecs(*X1);
1637  const int sizeX2 = MVT::GetNumberVecs(*X2);
1638 
1639  // Relative tolerance against which all tests are performed.
1640  // We are measuring things in the Frobenius norm $\| \cdot \|_F$.
1641  // The following bounds hold for all $m \times n$ matrices $A$:
1642  // \[
1643  // \|A\|_2 \leq \|A\|_F \leq \sqrt{r} \|A\|_2,
1644  // \]
1645  // where $r$ is the (column) rank of $A$. We bound this above
1646  // by the number of columns in $A$.
1647  //
1648  // Since we are measuring both of these quantitites in the
1649  // Frobenius norm instead, we scale all error tests by
1650  // $\sqrt{n}$.
1651  //
1652  // A good heuristic is to scale the bound by the square root
1653  // of the number of floating-point operations. One could
1654  // perhaps support this theoretically, since we are using
1655  // uniform random test problems.
1656  const magnitude_type fudgeFactor =
1657  SMT::squareroot(magnitude_type(numRows) *
1658  magnitude_type(numCols) *
1659  magnitude_type(numCols));
1660  const magnitude_type TOL = SMT::eps() * fudgeFactor *
1661  SMT::squareroot(magnitude_type(numCols));
1662 
1663  // Absolute tolerance scaling: the Frobenius norm of the test
1664  // matrix S. TOL*ATOL is the absolute tolerance for the
1665  // residual $\|A - Q*B\|_F$.
1666  const magnitude_type ATOL = frobeniusNorm (*S);
1667 
1668  sout << "The test matrix S has Frobenius norm " << ATOL
1669  << ", and the relative error tolerance is TOL = "
1670  << TOL << "." << endl;
1671 
1672 
1673  //
1674  // Output tests:
1675  // <S_out,X1> = 0
1676  // <S_out,X2> = 0
1677  // S_in = S_out + X1 C1 + X2 C2
1678  //
1679  // We will loop over an integer specifying the test combinations.
1680  // The bit pattern for the different tests is listed in parentheses.
1681  //
1682  // For the projectors, test the following combinations:
1683  // none (00)
1684  // P_X1 (01)
1685  // P_X2 (10)
1686  // P_X1 P_X2 (11)
1687  // P_X2 P_X1 (11)
1688  // The latter two should be tested to give the same result.
1689  //
1690  // For each of these, we should test with C1 and C2:
1691  //
1692  // if hasM:
1693  // with and without MX1 (1--)
1694  // with and without MX2 (1---)
1695  // with and without MS (1----)
1696  //
1697  // As hasM controls the upper level bits, we need only run test
1698  // cases 0-3 if hasM==false. Otherwise, we run test cases 0-31.
1699  //
1700 
1701  int numtests = 8;
1702 
1703  // test ortho error before orthonormalizing
1704  if (X1 != null) {
1705  magnitude_type err = OM->orthogError(*S,*X1);
1706  sout << " || <S,X1> || before : " << err << endl;
1707  }
1708  if (X2 != null) {
1709  magnitude_type err = OM->orthogError(*S,*X2);
1710  sout << " || <S,X2> || before : " << err << endl;
1711  }
1712 
1713  for (int t = 0; t < numtests; ++t)
1714  {
1715  Array< RCP< const MV > > theX;
1716  Array< RCP< mat_type > > C;
1717  if ( (t % 3) == 0 ) {
1718  // neither X1 nor X2
1719  // C and theX are already empty
1720  }
1721  else if ( (t % 3) == 1 ) {
1722  // X1
1723  theX = tuple(X1);
1724  C = tuple( rcp(new mat_type(sizeX1,sizeS)) );
1725  }
1726  else if ( (t % 3) == 2 ) {
1727  // X2
1728  theX = tuple(X2);
1729  C = tuple( rcp(new mat_type(sizeX2,sizeS)) );
1730  }
1731  else {
1732  // X1 and X2, and the reverse.
1733  theX = tuple(X1,X2);
1734  C = tuple( rcp(new mat_type(sizeX1,sizeS)),
1735  rcp(new mat_type(sizeX2,sizeS)) );
1736  }
1737 
1738  try {
1739  // call routine
1740  // if (t && 3) == 3, {
1741  // call with reversed input: X2 X1
1742  // }
1743  // test all outputs for correctness
1744  // test all outputs for equivalence
1745 
1746  // here is where the outputs go
1747  Array< RCP< MV > > S_outs;
1748  Array< Array< RCP< mat_type > > > C_outs;
1749  RCP< MV > Scopy;
1750 
1751  // copies of S,MS
1752  Scopy = MVT::CloneCopy(*S);
1753  // randomize this data, it should be overwritten
1754  for (size_type i = 0; i < C.size(); ++i) {
1755  Teuchos::randomSyncedMatrix(*C[i]);
1756  }
1757  // Run test.
1758  // Note that Anasazi and Belos differ, among other places,
1759  // in the order of arguments to project().
1760  OM->project(*Scopy,C,theX);
1761  // we allocate S and MS for each test, so we can save these as views
1762  // however, save copies of the C
1763  S_outs.push_back( Scopy );
1764  C_outs.push_back( Array< RCP< mat_type > >(0) );
1765  if (C.size() > 0) {
1766  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
1767  }
1768  if (C.size() > 1) {
1769  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
1770  }
1771 
1772  // do we run the reversed input?
1773  if ( (t % 3) == 3 ) {
1774  // copies of S,MS
1775  Scopy = MVT::CloneCopy(*S);
1776  // randomize this data, it should be overwritten
1777  for (size_type i = 0; i < C.size(); ++i) {
1778  Teuchos::randomSyncedMatrix(*C[i]);
1779  }
1780  // flip the inputs
1781  theX = tuple( theX[1], theX[0] );
1782  // Run test.
1783  // Note that Anasazi and Belos differ, among other places,
1784  // in the order of arguments to project().
1785  OM->project(*Scopy,C,theX);
1786  // we allocate S and MS for each test, so we can save these as views
1787  // however, save copies of the C
1788  S_outs.push_back( Scopy );
1789  // we are in a special case: P_X1 and P_X2, so we know we applied
1790  // two projectors, and therefore have two C[i]
1791  C_outs.push_back( Array<RCP<mat_type > >() );
1792  // reverse the Cs to compensate for the reverse projectors
1793  C_outs.back().push_back( rcp( new mat_type(*C[1]) ) );
1794  C_outs.back().push_back( rcp( new mat_type(*C[0]) ) );
1795  // flip the inputs back
1796  theX = tuple( theX[1], theX[0] );
1797  }
1798 
1799  // test all outputs for correctness
1800  for (size_type o = 0; o < S_outs.size(); ++o) {
1801  // S_in = X1*C1 + C2*C2 + S_out
1802  {
1803  RCP<MV> tmp = MVT::CloneCopy(*S_outs[o]);
1804  if (C_outs[o].size() > 0) {
1805  MVT::MvTimesMatAddMv(ONE,*X1,*C_outs[o][0],ONE,*tmp);
1806  if (C_outs[o].size() > 1) {
1807  MVT::MvTimesMatAddMv(ONE,*X2,*C_outs[o][1],ONE,*tmp);
1808  }
1809  }
1810  magnitude_type err = MVDiff(*tmp,*S);
1811  if (err > ATOL*TOL) {
1812  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1813  numerr++;
1814  }
1815  sout << " " << t << "|| S_in - X1*C1 - X2*C2 - S_out || : " << err << endl;
1816  }
1817  // <X1,S> == 0
1818  if (theX.size() > 0 && theX[0] != null) {
1819  magnitude_type err = OM->orthogError(*theX[0],*S_outs[o]);
1820  if (err > TOL) {
1821  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1822  numerr++;
1823  }
1824  sout << " " << t << "|| <X[0],S> || after : " << err << endl;
1825  }
1826  // <X2,S> == 0
1827  if (theX.size() > 1 && theX[1] != null) {
1828  magnitude_type err = OM->orthogError(*theX[1],*S_outs[o]);
1829  if (err > TOL) {
1830  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1831  numerr++;
1832  }
1833  sout << " " << t << "|| <X[1],S> || after : " << err << endl;
1834  }
1835  }
1836 
1837  // test all outputs for equivalence
1838  // check all combinations:
1839  // output 0 == output 1
1840  // output 0 == output 2
1841  // output 1 == output 2
1842  for (size_type o1=0; o1<S_outs.size(); o1++) {
1843  for (size_type o2=o1+1; o2<S_outs.size(); o2++) {
1844  // don't need to check MS_outs because we check
1845  // S_outs and MS_outs = M*S_outs
1846  // don't need to check C_outs either
1847  //
1848  // check that S_outs[o1] == S_outs[o2]
1849  magnitude_type err = MVDiff(*S_outs[o1],*S_outs[o2]);
1850  if (err > TOL) {
1851  sout << " vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv tolerance exceeded! test failed!" << endl;
1852  numerr++;
1853  }
1854  }
1855  }
1856 
1857  }
1858  catch (Belos::OrthoError& e) {
1859  sout << " ------------------------------------------- project() threw exception" << endl;
1860  sout << " Error: " << e.what() << endl;
1861  numerr++;
1862  }
1863  } // test for
1864 
1865  MsgType type = Debug;
1866  if (numerr>0) type = Errors;
1867  MyOM->stream(type) << sout.str();
1868  MyOM->stream(type) << endl;
1869 
1870  return numerr;
1871  }
1872 
1873 
1874  };
1875 
1876 
1877 
1878  } // namespace Test
1879 } // namespace Belos
1880 
1881 
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:262
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 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.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Wrapper around OrthoManager test functionality.
bool is_null() const

Generated on Fri Jun 5 2020 10:20:43 for Belos by doxygen 1.8.5