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