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