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