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

Generated for Belos by doxygen 1.10.0