Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_MultiVectorStdOpsTester_def.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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 Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_MULTI_VECTOR_STD_OPS_TESTER_HPP
43#define THYRA_MULTI_VECTOR_STD_OPS_TESTER_HPP
44
45#include "Thyra_MultiVectorStdOpsTester_decl.hpp"
46#include "Thyra_MultiVectorStdOps.hpp"
47#include "Thyra_VectorStdOps.hpp"
48#include "Thyra_TestingTools.hpp"
49
50namespace Thyra {
51
52// MultiVectorStdOpsTester
53
54template <class Scalar>
56 const ScalarMag &warning_tol_in
57 ,const ScalarMag &error_tol_in
58 ,const int num_mv_cols_in
59 )
60 :warning_tol_(warning_tol_in)
61 ,error_tol_(error_tol_in)
62 ,num_mv_cols_(num_mv_cols_in)
63{}
64
65template <class Scalar>
67 const VectorSpaceBase<Scalar> &vecSpc
68 ,std::ostream *out
69 ,const bool &/* dumpAll */
70 )
71{
72 using Teuchos::as;
73 using Teuchos::ptr;
74 using Teuchos::tuple;
76
77 if(out)
78 *out << "\n*** Entering MultiVectorStdOpsTester<"<<ST::name()<<">::checkStdOps(...) ...\n"
79 << "using a \'" << vecSpc.description() << "\' object ...\n";
80
81 bool success = true;
82 if(out) *out << "\nvecSpc.dim() = " << vecSpc.dim() << std::endl;
83
84 const Ordinal n = vecSpc.dim();
85
86 const Scalar
87 two = as<Scalar>(2.0),
88 three = as<Scalar>(3.0),
89 four = as<Scalar>(4.0);
90
91 int tc = 0;
92
93 if(out) *out << "\nCreating MultiVectorBase objects V1, V2, and V3 ...\n";
95 V1 = createMembers(vecSpc,num_mv_cols()),
96 V2 = createMembers(vecSpc,num_mv_cols()),
97 V3 = createMembers(vecSpc,num_mv_cols());
98
99 if(out) *out << "\nassign(V1.ptr(),-2.0);\n";
100 assign(V1.ptr(),Scalar(-two));
101
102 Teuchos::Array<Scalar> scalars1(num_mv_cols());
103 Teuchos::Array<Scalar> scalars2(num_mv_cols());
104 Teuchos::Array<ScalarMag> mags1(num_mv_cols());
105 Teuchos::Array<ScalarMag> mags2(num_mv_cols());
106
107 // sums
108 if(out) *out << "\n"<<tc<<") sums(*V1);\n";
109 ++tc;
110 {
111 sums(*V1, scalars1());
112 for (Ordinal i = 0; i < scalars2.size(); ++i)
113 scalars2[i] = -two*as<Scalar>(vecSpc.dim());
114 if (!testRelErrors<Scalar, Scalar, ScalarMag>(
115 "sums(*V1)", scalars1(),
116 "-2.0*n", scalars2(),
117 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
118 )
119 ) success = false;
120 }
121
122 // norms_1
123 if(out) *out << "\n"<<tc<<") norms_1(*V1);\n";
124 ++tc;
125 {
126 norms_1(*V1, mags1());
127 for (Ordinal i = 0; i < mags2.size(); ++i)
128 mags2[i] = ST::magnitude(two)*as<ScalarMag>(vecSpc.dim());
129 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
130 "norms_1(*V1)", mags1(),
131 "2.0*n", mags2(),
132 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
133 )
134 ) success = false;
135 }
136
137 // norms_2
138 if(out) *out << "\n"<<tc<<") norms_2(*V1);\n";
139 ++tc;
140 {
141 norms_2(*V1, mags1());
142 for (Ordinal i = 0; i < mags2.size(); ++i)
143 mags2[i] = ST::magnitude(two * ST::squareroot(as<Scalar>(n)));
144 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
145 "norms_2(*V1)", mags1(),
146 "2.0*sqrt(n)", mags2(),
147 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
148 )
149 ) success = false;
150 }
151
152 // norms_inf
153 if(out) *out << "\n"<<tc<<") norms_inf(*V1);\n";
154 ++tc;
155 {
156 norms_inf(*V1, mags1());
157 for (Ordinal i = 0; i < mags2.size(); ++i)
158 mags2[i] = ST::magnitude(two);
159 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
160 "norms_inf(*V1)", mags1(),
161 "2.0", mags2(),
162 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
163 )
164 ) success = false;
165 }
166
167 // assign(scalar)
168 if(out) *out << "\n"<<tc<<") assign(V2.ptr(), alpha);\n";
169 ++tc;
170 {
171 assign(V2.ptr(), three);
172 norms_2(*V2, mags1());
173 for (Ordinal i = 0; i < mags2.size(); ++i)
174 mags2[i] = ST::magnitude(three * ST::squareroot(as<Scalar>(n)));
175 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
176 "norms_2(*V2)", mags1(),
177 "3.0*sqrt(n)", mags2(),
178 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
179 )
180 ) success = false;
181 }
182
183 // assign(MV)
184 if(out) *out << "\n"<<tc<<") assign(V2.ptr(), *V1);\n";
185 ++tc;
186 assign(V2.ptr(), *V1);
187 norms_2(*V1, mags1());
188 norms_2(*V2, mags2());
189 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
190 "norms_2(*V1)", mags1(),
191 "norms_2(*V2)", mags2(),
192 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
193 )
194 ) success = false;
195
196 // scale
197 if(out) *out << "\n"<<tc<<") scale(alpha,V2.ptr());\n";
198 ++tc;
199 {
200 Scalar alpha = as<Scalar>(1.2345);
201 assign(V2.ptr(), *V1);
202 scale(alpha, V2.ptr());
203 norms_2(*V1, mags1());
204 for (Ordinal i = 0; i < mags1.size(); ++i)
205 mags1[i] *= ST::magnitude(alpha);
206 norms_2(*V2, mags2());
207 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
208 "norms_2(alpha*V1)", mags1(),
209 "alpha*norms_2(V1)", mags2(),
210 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
211 )
212 ) success = false;
213 }
214
215 // scaleUpdate
216 if(out) *out << "\n"<<tc<<") scaleUpdate(a,V1,V2.ptr());\n";
217 ++tc;
218 {
219 Teuchos::RCP<VectorBase<Scalar> > a = createMember(vecSpc);
220 assign(a.ptr(), two);
221 assign(V2.ptr(), four);
222 scaleUpdate(*a, *V1, V2.ptr()); // V2(i,j) = 2.0*(-2.0) + 4.0
223 norms_2(*V2, mags1());
224 if (!testMaxErrors<Scalar>(
225 "norms_2(*V2)", mags1(),
226 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
227 )
228 ) success = false;
229 }
230
231 // update(alpha, U, V.ptr())
232 if(out) *out << "\n"<<tc<<") update(a,V1,V2.ptr());\n";
233 ++tc;
234 {
235 Scalar alpha = as<Scalar>(1.2345);
236 assign(V2.ptr(), three);
237 assign(V3.ptr(), *V1);
238 scale(alpha, V3.ptr());
239 Vp_V(V3.ptr(), *V2);
240 update(alpha, *V1, V2.ptr());
241 norms_2(*V2, mags1());
242 norms_2(*V3, mags2());
243 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
244 "norms_2(*V2)", mags1(),
245 "norms_2(*V3)", mags2(),
246 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
247 )
248 ) success = false;
249 }
250
251 // update(alpha, beta, U, V.ptr())
252 if(out) *out << "\n"<<tc<<") update(alpha,beta,*V1,V2.ptr());\n";
253 ++tc;
254 {
255 Teuchos::Array<Scalar> alpha(num_mv_cols());
256 for (Ordinal i = 0; i < alpha.size(); ++i)
257 alpha[i] = as<Scalar>(i+1);
258 Scalar beta = as<Scalar>(1.2345);
259 assign(V2.ptr(), three);
260 assign(V3.ptr(), *V1);
261 scale(beta, V3.ptr());
262 for (Ordinal i = 0; i < alpha.size(); ++i)
263 scale(alpha[i], V3->col(i).ptr());
264 Vp_V(V3.ptr(), *V2);
265 Teuchos::ArrayView<const Scalar> alphaView = alpha();
266 update(alphaView, beta, *V1, V2.ptr());
267 norms_2(*V2, mags1());
268 norms_2(*V3, mags2());
269 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
270 "norms_2(*V2)", mags1(),
271 "norms_2(*V3)", mags2(),
272 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
273 )
274 ) success = false;
275 }
276
277 // update(U, alpha, beta, V.ptr())
278 if(out) *out << "\n"<<tc<<") update(*V1,alpha,beta,V2.ptr());\n";
279 ++tc;
280 {
281 Teuchos::Array<Scalar> alpha(num_mv_cols());
282 for (Ordinal i = 0; i < alpha.size(); ++i)
283 alpha[i] = as<Scalar>(i+1);
284 Scalar beta = as<Scalar>(1.2345);
285 assign(V2.ptr(), three);
286 assign(V3.ptr(), *V2);
287 scale(beta, V3.ptr());
288 for (Ordinal i = 0; i < alpha.size(); ++i)
289 scale(alpha[i], V3->col(i).ptr());
290 Vp_V(V3.ptr(), *V1);
291 Teuchos::ArrayView<const Scalar> alphaView = alpha();
292 update(*V1, alphaView, beta, V2.ptr());
293 norms_2(*V2, mags1());
294 norms_2(*V3, mags2());
295 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
296 "norms_2(*V2)", mags1(),
297 "norms_2(*V3)", mags2(),
298 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
299 )
300 ) success = false;
301 }
302
303 // linear_combination
304 if(out) *out << "\n"<<tc<<") linear_combination({alpha,beta,gamma},{V1.ptr(),V2.ptr(),V3.ptr()},0.0,V4.ptr());\n";
305 ++tc;
306 {
307 Scalar alpha = two, beta = -three, gamma = three;
308 Teuchos::RCP<MultiVectorBase<Scalar> > V4 = createMembers(vecSpc,num_mv_cols());
309 assign(V2.ptr(), two);
310 assign(V3.ptr(), four);
311 linear_combination<Scalar>(
312 tuple<Scalar>(alpha, beta, gamma),
313 tuple<Ptr<const MultiVectorBase<Scalar> > >(V1.ptr(), V2.ptr(), V3.ptr()),
314 as<Scalar>(0.0),
315 V4.ptr()); // V4(i,j) = 2.0(-2.0) + (-3.0)(2.0) + 3.0(4.0)
316 norms_2(*V4, mags1());
317 for (Ordinal i = 0; i < mags2.size(); ++i)
318 mags2[i] = ST::magnitude(two * ST::squareroot(as<Scalar>(n)));
319 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
320 "norms_2(*V4)", mags1(),
321 "2.0*sqrt(n)", mags2(),
322 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
323 )
324 ) success = false;
325 }
326
327 if(out) *out << "\n"<<tc<<") linear_combination({alpha,beta,gamma},{V1.ptr(),V2.ptr(),V3.ptr()},0.5,V4.ptr());\n";
328 ++tc;
329 {
330 Scalar alpha = two, beta = -three, gamma = three;
331 Teuchos::RCP<MultiVectorBase<Scalar> > V4 = createMembers(vecSpc,num_mv_cols());
332 assign(V2.ptr(), two);
333 assign(V3.ptr(), four);
334 assign(V4.ptr(), -four);
335 linear_combination<Scalar>(
336 tuple<Scalar>(alpha, beta, gamma),
337 tuple<Ptr<const MultiVectorBase<Scalar> > >(V1.ptr(), V2.ptr(), V3.ptr()),
338 as<Scalar>(0.5),
339 V4.ptr()); // V4(i,j) = 0.5(-4.0) + 2.0(-2.0) + (-3.0)(2.0) + 3.0(4.0)
340 norms_2(*V4, mags1());
341 if (!testMaxErrors<Scalar>(
342 "norms_2(*V4)", mags1(),
343 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
344 )
345 ) success = false;
346 }
347
348 // Vt_S
349 if(out) *out << "\n"<<tc<<") Vt_S(V1.ptr(),alpha);\n";
350 ++tc;
351 {
352 Scalar alpha = as<Scalar>(1.2345);
353 assign(V2.ptr(), *V1);
354 Vt_S(V2.ptr(), alpha);
355 norms_2(*V1, mags1());
356 for (Ordinal i = 0; i < mags1.size(); ++i)
357 mags1[i] *= ST::magnitude(alpha);
358 norms_2(*V2, mags2());
359 if (!testRelErrors<ScalarMag, ScalarMag, ScalarMag>(
360 "norms_2(alpha*V1)", mags1(),
361 "alpha*norms_2(V1)", mags2(),
362 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
363 )
364 ) success = false;
365 }
366
367 // Vp_S
368 if(out) *out << "\n"<<tc<<") Vp_S(V2.ptr(),alpha);\n";
369 ++tc;
370 {
371 assign(V2.ptr(), *V1);
372 Vp_S(V2.ptr(), two); // V2(i,j) = -2.0 + 2.0
373 norms_2(*V2, mags1());
374 if (!testMaxErrors<Scalar>(
375 "norms_2(V2)", mags1(),
376 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
377 )
378 ) success = false;
379 }
380
381 // Vp_V
382 if(out) *out << "\n"<<tc<<") Vp_V(V2.ptr(),*V1);\n";
383 ++tc;
384 {
385 assign(V2.ptr(), two);
386 Vp_V(V2.ptr(), *V1); // V2(i,j) = -2.0 + 2.0
387 norms_2(*V2, mags1());
388 if (!testMaxErrors<Scalar>(
389 "norms_2(V2)", mags1(),
390 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
391 )
392 ) success = false;
393 }
394
395
396 // V_VpV
397 if(out) *out << "\n"<<tc<<") V_VpV(V3.ptr(),*V1,*V2);\n";
398 ++tc;
399 {
400 assign(V2.ptr(), two);
401 V_VpV(V3.ptr(), *V1, *V2); // V3(i,j) = -2.0 + 2.0
402 norms_2(*V3, mags1());
403 if (!testMaxErrors<Scalar>(
404 "norms_2(V3)", mags1(),
405 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
406 )
407 ) success = false;
408 }
409
410
411 // V_VmV
412 if(out) *out << "\n"<<tc<<") V_VmV(V3.ptr(),*V1,*V2);\n";
413 ++tc;
414 {
415 assign(V2.ptr(), -two);
416 V_VmV(V3.ptr(), *V1, *V2); // V3(i,j) = -2.0 - (-2.0)
417 norms_2(*V3, mags1());
418 if (!testMaxErrors<Scalar>(
419 "norms_2(V3)", mags1(),
420 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
421 )
422 ) success = false;
423 }
424
425
426 // V_StVpV
427 if(out) *out << "\n"<<tc<<") V_StVpV(V3.ptr(),alpha,*V1,*V2);\n";
428 ++tc;
429 {
430 Scalar alpha = as<Scalar>(1.2345);
431 assign(V2.ptr(), three);
432 V_StVpV(V3.ptr(), alpha, *V1, *V2);
433 scale(alpha, V1.ptr());
434 Vp_V(V2.ptr(), *V1);
435 V_VmV(V3.ptr(), *V2, *V3);
436 norms_2(*V3, mags1());
437 if (!testMaxErrors<Scalar>(
438 "norms_2(V3)", mags1(),
439 "error_tol", error_tol(), "warning_tol", warning_tol(), ptr(out)
440 )
441 ) success = false;
442 }
443
444
445 if(out) *out
446 << "\n*** Leaving MultiVectorStdOpsTester<"<<ST::name()<<">::checkStdOps(...) ...\n";
447
448 return success;
449
450}
451
452} // namespace Thyra
453
454#endif // THYRA_MULTI_VECTOR_STD_OPS_TESTER_HPP
size_type size() const
virtual std::string description() const
Ptr< T > ptr() const
Interface for a collection of column vectors called a multi-vector.
bool checkStdOps(const VectorSpaceBase< Scalar > &vecSpc, std::ostream *out=0, const bool &dumpAll=false)
Run the tests using a vector space.
MultiVectorStdOpsTester(const ScalarMag &warning_tol=0, const ScalarMag &error_tol=0, const int num_mv_cols=4)
Teuchos::ScalarTraits< Scalar >::magnitudeType ScalarMag
Abstract interface for objects that represent a space for vectors.
virtual Ordinal dim() const =0
Return the dimension of the vector space.
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
TypeTo as(const TypeFrom &t)