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