1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
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 Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41
42 #include "Teuchos_SerialDenseMatrix.hpp"
43 #include "Teuchos_SerialBandDenseMatrix.hpp"
44 #include "Teuchos_SerialDenseVector.hpp"
45 #include "Teuchos_SerialDenseHelpers.hpp"
46 #include "Teuchos_SerialBandDenseSolver.hpp"
47 #include "Teuchos_ScalarTraits.hpp"
48 #include "Teuchos_RCP.hpp"
49 #include "Teuchos_Version.hpp"
50
51 using Teuchos::ScalarTraits;
52 using Teuchos::SerialDenseMatrix;
53 using Teuchos::SerialBandDenseMatrix;
54 using Teuchos::SerialDenseVector;
55
56 #define OTYPE int
57 #ifdef HAVE_TEUCHOS_COMPLEX
58 #define STYPE std::complex<double>
59 #else
60 #define STYPE double
61 #endif
62
63 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and vector elements and scalars:
64 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
65 #ifdef HAVE_TEUCHOS_COMPLEX
66 #define SCALARMAX STYPE(10,0)
67 #else
68 #define SCALARMAX STYPE(10)
69 #endif
70
71 template<typename TYPE>
72 int PrintTestResults(std::string, TYPE, TYPE, bool);
73
74 int ReturnCodeCheck(std::string, int, int, bool);
75
76 typedef SerialDenseVector<OTYPE, STYPE> DVector;
77 typedef SerialDenseMatrix<OTYPE, STYPE> DMatrix;
78 typedef SerialBandDenseMatrix<OTYPE, STYPE> BDMatrix;
79
80 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
81 template<typename TYPE>
82 TYPE GetRandom(TYPE, TYPE);
83
84 // Returns a random integer between the two input parameters, inclusive
85 template<>
86 int GetRandom(int, int);
87
88 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
89 template<>
90 double GetRandom(double, double);
91
92 template<typename T>
93 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
94
95 // Generates random matrices and vectors using GetRandom()
96 Teuchos::RCP<DMatrix> GetRandomBandedMatrix(int m, int n, int kl, int ku);
97 Teuchos::RCP<DVector> GetRandomVector(int n);
98
99 // Compares the difference between two vectors using relative euclidean norms
100 // Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
101 int CompareVectors(const SerialDenseVector<OTYPE,STYPE>& Vector1,
102 const SerialDenseVector<OTYPE,STYPE>& Vector2,
103 ScalarTraits<STYPE>::magnitudeType Tolerance );
104
105 // Compares the difference between two matrices using relative euclidean norms
106 // Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
107 int CompareMatrices(const SerialDenseMatrix<OTYPE,STYPE>& Matrix1,
108 const SerialDenseMatrix<OTYPE,STYPE>& Matrix2,
109 ScalarTraits<STYPE>::magnitudeType Tolerance );
110
main(int argc,char * argv[])111 int main(int argc, char* argv[])
112
113 {
114
115 typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
116
117 int n=10, kl=2, ku=3;
118 MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
119
120 bool verbose = 0;
121 if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
122
123 if (verbose)
124 std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
125
126 int numberFailedTests = 0;
127 int returnCode = 0;
128 std::string testName = "", testType = "";
129
130 #ifdef HAVE_TEUCHOS_COMPLEX
131 testType = "COMPLEX";
132 #else
133 testType = "REAL";
134 #endif
135
136 if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL BAND DENSE SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
137
138 // Create dense matrix and vector.
139 Teuchos::RCP<DMatrix> A1 = GetRandomBandedMatrix(n,n,kl,ku);
140 Teuchos::RCP<DVector> x1 = GetRandomVector(n);
141 DVector xhat(n), b(n), bt(n);
142
143 // Create a serial dense solver.
144 Teuchos::SerialBandDenseSolver<OTYPE, STYPE> solver1;
145
146 Teuchos::RCP<BDMatrix> AB1;
147 Teuchos::RCP<DMatrix> C1;
148
149 // Convert the dense matrix to a matrix in LAPACK banded storage
150 AB1 = Teuchos::generalToBanded( A1, kl, ku, true );
151
152 // Convert the matrix in LAPACK banded storage back to a normal serial dense matrix
153 C1 = Teuchos::bandedToGeneral( AB1 );
154 testName = "Converting matrix formats: generalToBanded and bandedToGeneral random A:";
155 numberFailedTests += CompareMatrices( *A1, *C1, tol );
156 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
157
158 // Compute the right-hand side vector using multiplication.
159 returnCode = b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, *x1, ScalarTraits<STYPE>::zero());
160 testName = "Generating right-hand side vector using A*x, where x is a random vector:";
161 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
162 returnCode = bt.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, *x1, ScalarTraits<STYPE>::zero());
163 testName = "Generating right-hand side vector using A^T*x, where x is a random vector:";
164 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
165
166 #ifdef HAVE_TEUCHOS_COMPLEX
167 DVector bct(n);
168 returnCode = bct.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A1, *x1, ScalarTraits<STYPE>::zero());
169 testName = "Generating right-hand side vector using A^H*x, where x is a random vector:";
170 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
171 #endif
172
173 // Fill the solution vector with zeros.
174 xhat.putScalar( ScalarTraits<STYPE>::zero() );
175
176 // Pass in matrix and vectors
177 solver1.setMatrix( AB1 );
178 solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
179
180 // Test1: Simple factor and solve
181 returnCode = solver1.factor();
182 testName = "Simple solve: factor() random A:";
183 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
184
185 // Non-transpose solve
186 returnCode = solver1.solve();
187 testName = "Simple solve: solve() random A (NO_TRANS):";
188 numberFailedTests += CompareVectors( *x1, xhat, tol );
189 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
190
191 // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
192 xhat.putScalar( ScalarTraits<STYPE>::zero() );
193 solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
194 solver1.solveWithTransposeFlag( Teuchos::TRANS );
195 returnCode = solver1.solve();
196 testName = "Simple solve: solve() random A (TRANS):";
197 numberFailedTests += CompareVectors( *x1, xhat, tol );
198 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
199
200 #ifdef HAVE_TEUCHOS_COMPLEX
201 // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
202 xhat.putScalar( ScalarTraits<STYPE>::zero() );
203 solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
204 solver1.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
205 returnCode = solver1.solve();
206 testName = "Simple solve: solve() random A (CONJ_TRANS):";
207 numberFailedTests += CompareVectors( *x1, xhat, tol );
208 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
209 #endif
210
211 // Test2: Solve with iterative refinement.
212 #ifdef HAVE_TEUCHOSNUMERICS_EIGEN
213 // Iterative refinement not implemented in Eigen
214 #else
215 // Create random linear system
216 Teuchos::RCP<DMatrix> A2 = GetRandomBandedMatrix(n,n,kl,ku);
217 Teuchos::RCP<DVector> x2 = GetRandomVector(n);
218
219 // Create LHS through multiplication with A2
220 xhat.putScalar( ScalarTraits<STYPE>::zero() );
221 b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A2, *x2, ScalarTraits<STYPE>::zero());
222 bt.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A2, *x2, ScalarTraits<STYPE>::zero());
223 #ifdef HAVE_TEUCHOS_COMPLEX
224 bct.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A2, *x2, ScalarTraits<STYPE>::zero());
225 #endif
226
227 // Create a serial dense solver.
228 Teuchos::SerialBandDenseSolver<OTYPE, STYPE> solver2;
229 solver2.solveToRefinedSolution( true );
230
231 Teuchos::RCP<BDMatrix> AB2;
232 Teuchos::RCP<DMatrix> C2;
233 // Convert the dense matrix to a matrix in LAPACK banded storage
234 AB2 = Teuchos::generalToBanded( A2, kl, ku, true );
235
236 // Convert the matrix in LAPACK banded storage back to a normal serial dense matrix
237 C2 = Teuchos::bandedToGeneral( AB2 );
238 testName = "Converting matrix formats: generalToBanded and bandedToGeneral random A:";
239 numberFailedTests += CompareMatrices( *A2, *C2, tol );
240 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
241
242 // Pass in matrix and vectors
243 solver2.setMatrix( AB2 );
244 solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
245
246 // Factor and solve with iterative refinement.
247 returnCode = solver2.factor();
248 testName = "Solve with iterative refinement: factor() random A:";
249 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
250
251 // Non-transpose solve
252 returnCode = solver2.solve();
253 testName = "Solve with iterative refinement: solve() random A (NO_TRANS):";
254 numberFailedTests += CompareVectors( *x2, xhat, tol );
255 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
256
257 // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
258 xhat.putScalar( ScalarTraits<STYPE>::zero() );
259 solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
260 solver2.solveWithTransposeFlag( Teuchos::TRANS );
261 returnCode = solver2.solve();
262 testName = "Solve with iterative refinement: solve() random A (TRANS):";
263 numberFailedTests += CompareVectors( *x2, xhat, tol );
264 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
265
266 #ifdef HAVE_TEUCHOS_COMPLEX
267 // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
268 xhat.putScalar( ScalarTraits<STYPE>::zero() );
269 solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
270 solver2.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
271 returnCode = solver2.solve();
272 testName = "Solve with iterative refinement: solve() random A (CONJ_TRANS):";
273 numberFailedTests += CompareVectors( *x2, xhat, tol );
274 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
275 #endif
276 #endif
277
278 // Test3: Solve with matrix equilibration.
279
280 // Create random linear system
281 Teuchos::RCP<DMatrix> A3 = GetRandomBandedMatrix(n,n,kl,ku);
282 Teuchos::RCP<DVector> x3 = GetRandomVector(n);
283
284 // Create LHS through multiplication with A3
285 xhat.putScalar( ScalarTraits<STYPE>::zero() );
286 b.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A3, *x3, ScalarTraits<STYPE>::zero());
287 bt.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A3, *x3, ScalarTraits<STYPE>::zero());
288 #ifdef HAVE_TEUCHOS_COMPLEX
289 bct.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, ScalarTraits<STYPE>::one() , *A3, *x3, ScalarTraits<STYPE>::zero());
290 #endif
291
292 // Create a serial dense solver.
293 Teuchos::SerialBandDenseSolver<OTYPE, STYPE> solver3;
294 solver3.factorWithEquilibration( true );
295
296 Teuchos::RCP<BDMatrix> AB3;
297 Teuchos::RCP<DMatrix> C3;
298 // Convert the dense matrix to a matrix in LAPACK banded storage
299 AB3 = Teuchos::generalToBanded( A3, kl, ku, true );
300
301 // Convert the matrix in LAPACK banded storage back to a normal serial dense matrix
302 C3 = Teuchos::bandedToGeneral( AB3 );
303 testName = "Converting matrix formats: generalToBanded and bandedToGeneral random A:";
304 numberFailedTests += CompareMatrices( *A3, *C3, tol );
305 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
306
307 // Pass in matrix and vectors
308 solver3.setMatrix( AB3 );
309 solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
310
311 Teuchos::RCP<BDMatrix> AB3bak = Teuchos::rcp( new BDMatrix( *AB3 ) );
312 Teuchos::RCP<DVector> b3bak = Teuchos::rcp( new DVector( Teuchos::Copy, b ) );
313
314 // Factor and solve with matrix equilibration.
315 returnCode = solver3.factor();
316 testName = "Solve with matrix equilibration: factor() random A:";
317 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
318
319 // Non-transpose solve
320 returnCode = solver3.solve();
321 testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
322 numberFailedTests += CompareVectors( *x3, xhat, tol );
323 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
324
325 // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
326 xhat.putScalar( ScalarTraits<STYPE>::zero() );
327 solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bt, false ) );
328 solver3.solveWithTransposeFlag( Teuchos::TRANS );
329 returnCode = solver3.solve();
330 testName = "Solve with matrix equilibration: solve() random A (TRANS):";
331 numberFailedTests += CompareVectors( *x3, xhat, tol );
332 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
333
334 #ifdef HAVE_TEUCHOS_COMPLEX
335 // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
336 xhat.putScalar( ScalarTraits<STYPE>::zero() );
337 solver3.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &bct, false ) );
338 solver3.solveWithTransposeFlag( Teuchos::CONJ_TRANS );
339 returnCode = solver3.solve();
340 testName = "Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
341 numberFailedTests += CompareVectors( *x3, xhat, tol );
342 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
343 #endif
344
345 // Non-tranpose solve where factor is not called
346 xhat.putScalar( ScalarTraits<STYPE>::zero() );
347 solver3.setMatrix( AB3bak );
348 solver3.setVectors( Teuchos::rcp( &xhat, false ), b3bak );
349 solver3.solveWithTransposeFlag( Teuchos::NO_TRANS );
350 returnCode = solver3.solve();
351 testName = "Solve with matrix equilibration: solve() without factor() random A (NO_TRANS):";
352 numberFailedTests += CompareVectors( *x3, xhat, tol );
353 numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
354
355 //
356 // If a test failed output the number of failed tests.
357 //
358 if(numberFailedTests > 0)
359 {
360 if (verbose) {
361 std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
362 std::cout << "End Result: TEST FAILED" << std::endl;
363 return -1;
364 }
365 }
366 if(numberFailedTests == 0)
367 std::cout << "End Result: TEST PASSED!" << std::endl;
368
369 return 0;
370 }
371
372 template<typename TYPE>
PrintTestResults(std::string testName,TYPE calculatedResult,TYPE expectedResult,bool verbose)373 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
374 {
375 int result;
376 if(calculatedResult == expectedResult)
377 {
378 if(verbose) std::cout << testName << " successful." << std::endl;
379 result = 0;
380 }
381 else
382 {
383 if(verbose) std::cout << testName << " unsuccessful." << std::endl;
384 result = 1;
385 }
386 return result;
387 }
388
ReturnCodeCheck(std::string testName,int returnCode,int expectedResult,bool verbose)389 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
390 {
391 int result;
392 if(expectedResult == 0)
393 {
394 if(returnCode == 0)
395 {
396 if(verbose) std::cout << testName << " test successful." << std::endl;
397 result = 0;
398 }
399 else
400 {
401 if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
402 result = 1;
403 }
404 }
405 else
406 {
407 if(returnCode != 0)
408 {
409 if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
410 result = 0;
411 }
412 else
413 {
414 if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
415 result = 1;
416 }
417 }
418 return result;
419 }
420
421 template<typename TYPE>
GetRandom(TYPE Low,TYPE High)422 TYPE GetRandom(TYPE Low, TYPE High)
423 {
424 return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
425 }
426
427 template<typename T>
GetRandom(std::complex<T> Low,std::complex<T> High)428 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
429 {
430 T lowMag = Low.real();
431 T highMag = High.real();
432 T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
433 T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
434 return std::complex<T>( real, imag );
435 }
436
437 template<>
GetRandom(int Low,int High)438 int GetRandom(int Low, int High)
439 {
440 return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
441 }
442
443 template<>
GetRandom(double Low,double High)444 double GetRandom(double Low, double High)
445 {
446 return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
447 }
448
GetRandomBandedMatrix(int m,int n,int kl,int ku)449 Teuchos::RCP<DMatrix> GetRandomBandedMatrix(int m, int n, int kl, int ku)
450 {
451 Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
452
453 // Fill dense matrix with random entries.
454 for (int i=0; i<m; i++)
455 for (int j=0; j<n; j++)
456 if (j>= i-kl && j<=i+ku)
457 (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
458
459 return newmat;
460 }
461
GetRandomVector(int n)462 Teuchos::RCP<DVector> GetRandomVector(int n)
463 {
464 Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
465
466 // Fill dense vector with random entries.
467 for (int i=0; i<n; i++)
468 (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
469
470 return newvec;
471 }
472
473 /* Function: CompareVectors
474 Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
475 */
CompareVectors(const SerialDenseVector<OTYPE,STYPE> & Vector1,const SerialDenseVector<OTYPE,STYPE> & Vector2,ScalarTraits<STYPE>::magnitudeType Tolerance)476 int CompareVectors(const SerialDenseVector<OTYPE,STYPE>& Vector1,
477 const SerialDenseVector<OTYPE,STYPE>& Vector2,
478 ScalarTraits<STYPE>::magnitudeType Tolerance )
479 {
480 typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
481
482 SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
483 diff -= Vector2;
484
485 MagnitudeType norm_diff = diff.normFrobenius();
486 MagnitudeType norm_v2 = Vector2.normFrobenius();
487 MagnitudeType temp = norm_diff;
488 if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
489 temp /= norm_v2;
490
491 if (temp > Tolerance)
492 return 1;
493 else
494 return 0;
495 }
496
497 /* Function: CompareVectors
498 Purpose: Compares the difference between two matrices using relative euclidean-norms, i.e. ||m_1-m_2||_\inf/||m_2||_\inf
499 */
CompareMatrices(const SerialDenseMatrix<OTYPE,STYPE> & Matrix1,const SerialDenseMatrix<OTYPE,STYPE> & Matrix2,ScalarTraits<STYPE>::magnitudeType Tolerance)500 int CompareMatrices(const SerialDenseMatrix<OTYPE,STYPE>& Matrix1,
501 const SerialDenseMatrix<OTYPE,STYPE>& Matrix2,
502 ScalarTraits<STYPE>::magnitudeType Tolerance )
503 {
504 typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
505
506 SerialDenseMatrix<OTYPE,STYPE> diff( Matrix1 );
507 diff -= Matrix2;
508
509 MagnitudeType norm_diff = diff.normInf();
510 MagnitudeType norm_m2 = Matrix2.normInf();
511 MagnitudeType temp = norm_diff;
512 if (norm_m2 != ScalarTraits<MagnitudeType>::zero())
513 temp /= norm_m2;
514
515 if (temp > Tolerance)
516 return 1;
517 else
518 return 0;
519 }
520