Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
cxx_main_qr_solver.cpp
Go to the documentation of this file.
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 
46 #include "Teuchos_ScalarTraits.hpp"
47 #include "Teuchos_RCP.hpp"
48 #include "Teuchos_Version.hpp"
49 
51 using Teuchos::LAPACK;
54 
55 #define OTYPE int
56 #ifdef HAVE_TEUCHOS_COMPLEX
57 #define STYPE std::complex<double>
58 #else
59 #define STYPE double
60 #endif
61 
62 // SCALARMAX defines the maximum positive value (with a little leeway) generated for matrix and vector elements and scalars:
63 // random numbers in [-SCALARMAX, SCALARMAX] will be generated.
64 #ifdef HAVE_TEUCHOS_COMPLEX
65 #define SCALARMAX STYPE(10,0)
66 #else
67 #define SCALARMAX STYPE(10)
68 #endif
69 
70 template<typename TYPE>
71 int PrintTestResults(std::string, TYPE, TYPE, bool);
72 
73 int ReturnCodeCheck(std::string, int, int, bool);
74 
77 
78 // Returns ScalarTraits<TYPE>::random() (the input parameters are ignored)
79 template<typename TYPE>
80 TYPE GetRandom(TYPE, TYPE);
81 
82 // Returns a random integer between the two input parameters, inclusive
83 template<>
84 int GetRandom(int, int);
85 
86 // Returns a random double between the two input parameters, plus or minus a random number between 0 and 1
87 template<>
88 double GetRandom(double, double);
89 
90 template<typename T>
91 std::complex<T> GetRandom( std::complex<T>, std::complex<T> );
92 
93 // Generates random matrices and vectors using GetRandom()
96 
97 // Compares the difference between two vectors using relative euclidean norms
98 // Returns 1 if the comparison failed, the relative difference is greater than the tolerance.
100  const SerialDenseVector<OTYPE,STYPE>& Vector2,
102 
103 int main(int argc, char* argv[])
104 {
105  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
106 
108 
109  int n=8, m=12;
110  MagnitudeType tol = 1e-12*ScalarTraits<MagnitudeType>::one();
111 
112  bool verbose = 0;
113  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
114 
115  if (verbose)
116  std::cout << Teuchos::Teuchos_Version() << std::endl << std::endl;
117 
118  int numberFailedTests = 0;
119  int returnCode = 0;
120  std::string testName = "", testType = "";
121 
122 #ifdef HAVE_TEUCHOS_COMPLEX
123  testType = "COMPLEX";
124 #else
125  testType = "REAL";
126 #endif
127 
128  if (verbose) std::cout<<std::endl<<"********** CHECKING TEUCHOS SERIAL QR SOLVER - " << testType << "-VALUED **********"<<std::endl<<std::endl;
129 
130  // Create dense matrix and vector.
134  DVector xhat(n), b(m), xhatt(m);
135 
136  // Compute the right-hand side vector using multiplication.
138  testName = "Generating right-hand side vector using A*x, where x is a random vector:";
139  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
140 
141 #ifdef HAVE_TEUCHOS_COMPLEX
142  DVector bct(n);
144  testName = "Generating right-hand side vector using A^H*x, where x is a random vector:";
145  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
146 #else
147  DVector bt(n);
149  testName = "Generating right-hand side vector using A^T*x, where x is a random vector:";
150  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
151 #endif
152  DVector bp, bp2;
154  DVector tmp(n), tmp2(m);
155 
156  // Fill the solution vector with zeros.
157  xhat.putScalar( ScalarTraits<STYPE>::zero() );
159 
160  // Create a serial dense solver.
162 
163  // Pass in matrix and vectors
164  solver1.setMatrix( A1 );
165  solver1.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
166 
167  // Test1: Simple factor and solve
168  returnCode = solver1.factor();
169  testName = "Simple solve: factor() random A:";
170  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
171  returnCode = solver1.formQ();
172  testName = "Simple solve: formQ() random A:";
173  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
174 
175  // Non-transpose solve
176  returnCode = solver1.solve();
177  testName = "Simple solve: solve() random A (NO_TRANS):";
178  numberFailedTests += CompareVectors( *x1, xhat, tol );
179  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
180  testName = "Simple solve: multiplyQ(TRANS) solveR (NO_TRANS) random A (NO_TRANS):";
181  bp = b;
182  returnCode = solver1.multiplyQ( Teuchos::TRANS, bp );
183  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
184  for (OTYPE i=0; i<n; i++) {tmp(i) = bp(i);}
185  returnCode = solver1.solveR( Teuchos::NO_TRANS, tmp );
186  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
187  numberFailedTests += CompareVectors( tmp, xhat, tol );
188 
189 #ifdef HAVE_TEUCHOS_COMPLEX
190  testName = "Simple solve: formQ() solve with explicit Q implicit R random A (NO_TRANS):";
191  Q = solver1.getQ();
192  bp.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, 1.0, *Q, b, 0.0);
193  for (OTYPE i=0; i<n; i++) {tmp(i) = bp(i);}
194  returnCode = solver1.solveR( Teuchos::NO_TRANS, tmp );
195  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
196  numberFailedTests += CompareVectors( tmp, xhat, tol );
197  testName = "Simple solve: formQ() solve with explicit Q explicit R random A (NO_TRANS):";
198  Q = solver1.getQ();
199  returnCode = solver1.formR();
200  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
201  bp.multiply(Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, 1.0, *Q, b, 0.0);
202  for (OTYPE i=0; i<n; i++) {tmp(i) = bp(i);}
203  returnCode = solver1.solveR( Teuchos::NO_TRANS, tmp );
204  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
205  numberFailedTests += CompareVectors( tmp, xhat, tol );
206 #else
207  testName = "Simple solve: formQ() solve with explicit Q random A (NO_TRANS):";
208  Q = solver1.getQ();
209  bp.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *Q, b, 0.0);
210  for (OTYPE i=0; i<n; i++) {tmp(i) = bp(i);}
211  returnCode = solver1.solveR( Teuchos::NO_TRANS, tmp );
212  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
213  numberFailedTests += CompareVectors( tmp, xhat, tol );
214  testName = "Simple solve: formQ() solve with explicit Q explicit R random A (NO_TRANS):";
215  Q = solver1.getQ();
216  returnCode = solver1.formR();
217  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
218  bp.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, *Q, b, 0.0);
219  for (OTYPE i=0; i<n; i++) {tmp(i) = bp(i);}
220  returnCode = solver1.solveR( Teuchos::NO_TRANS, tmp );
221  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
222  numberFailedTests += CompareVectors( tmp, xhat, tol );
223 #endif
224 
225 #ifdef HAVE_TEUCHOS_COMPLEX
226  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
228  solver1.setVectors( Teuchos::rcp( &xhatt, false ), Teuchos::rcp( &bct, false ) );
230  returnCode = solver1.solve();
231  testName = "Simple solve: solve() random A (CONJ_TRANS):";
232  if (m == n)
233  numberFailedTests += CompareVectors( *x1t, xhatt, tol );
234  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
235  testName = "Simple solve: solveR (NO_TRANS) multiplyQ(TRANS) random A (CONJ_TRANS):";
236  bp = bct;
237  returnCode = solver1.solveR( Teuchos::TRANS, bp );
238  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
239  for (OTYPE i=0; i<n; i++) {tmp2(i) = bp(i);}
240  returnCode = solver1.multiplyQ( Teuchos::NO_TRANS, tmp2 );
241  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
242  if (m == n)
243  numberFailedTests += CompareVectors(*x1t, tmp2, tol );
244  numberFailedTests += CompareVectors( tmp2, xhatt, tol );
245  testName = "Simple solve: formQ() solve with explicit Q random A (CONJ_TRANS):";
246  bp = bct;
247  Q = solver1.getQ();
248  returnCode = solver1.solveR( Teuchos::TRANS, bp );
249  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
250  tmp2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *Q, bp, 0.0);
251  if (m == n)
252  numberFailedTests += CompareVectors(*x1t, tmp2, tol );
253  numberFailedTests += CompareVectors( tmp2, xhatt, tol );
254 
255 #else
256  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
258  solver1.setVectors( Teuchos::rcp( &xhatt, false ), Teuchos::rcp( &bt, false ) );
260  returnCode = solver1.solve();
261  testName = "Simple solve: solve() random A (TRANS):";
262  if (m == n)
263  numberFailedTests += CompareVectors( *x1t, xhatt, tol );
264  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
265  bp = bt;
266  testName = "Simple solve: solveR multiplyQ(TRANS) (NO_TRANS) random A (TRANS):";
267  returnCode = solver1.solveR( Teuchos::TRANS, bp );
268  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
269  for (OTYPE i=0; i<n; i++) {tmp2(i) = bp(i);}
270  returnCode = solver1.multiplyQ( Teuchos::NO_TRANS, tmp2 );
271  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
272  if (m == n)
273  numberFailedTests += CompareVectors(*x1t, tmp2, tol );
274  numberFailedTests += CompareVectors( tmp2, xhatt, tol );
275  testName = "Simple solve: formQ() solve with explicit Q random A (CONJ_TRANS):";
276  bp = bt;
277  Q = solver1.getQ();
278  returnCode = solver1.solveR( Teuchos::TRANS, bp );
279  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
280  tmp2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, *Q, bp, 0.0);
281  if (m == n)
282  numberFailedTests += CompareVectors(*x1t, tmp2, tol );
283  numberFailedTests += CompareVectors( tmp2, xhatt, tol );
284 
285 #endif
286 
287  // Test2: Solve with matrix equilibration.
288 
289  // Create random linear system
293 
294  // Create LHS through multiplication with A2
295  xhat.putScalar( ScalarTraits<STYPE>::zero() );
298 #ifdef HAVE_TEUCHOS_COMPLEX
300 #else
302 #endif
303 
304  // Create a serial dense solver.
306  solver2.factorWithEquilibration( true );
307 
308  // Pass in matrix and vectors
309  MagnitudeType safeMin = ScalarTraits<STYPE>::sfmin();
310  MagnitudeType prec = ScalarTraits<STYPE>::prec();
312  MagnitudeType smlnum = ScalarTraits<STYPE>::magnitude(safeMin/prec);
313  MagnitudeType bignum = ScalarTraits<STYPE>::magnitude(sOne/smlnum);
314  (void) bignum; // Silence "unused variable" compiler warning.
315  MagnitudeType smlnum2 = smlnum/2;
316  (void) smlnum2; // Silence "unused variable" compiler warning.
318  for (OTYPE j = 0; j < A2->numCols(); j++) {
319  for (OTYPE i = 0; i < A2->numRows(); i++) {
320  anorm = TEUCHOS_MAX( anorm, ScalarTraits<STYPE>::magnitude((*A2)(i,j)) );
321  }
322  }
323  OTYPE BW = 0;
324  (void) BW; // Silence "unused variable" compiler warning.
325  OTYPE info = 0;
326  (void) info; // Silence "unused variable" compiler warning.
327  // TODO: fix scaling test
328  // L.LASCL(Teuchos::ETypeChar[Teuchos::FULL], BW, BW, anorm, smlnum2, A2->numRows(), A2->numCols(), A2->values(), A2->stride(), &info);
329  solver2.setMatrix( A2 );
330  solver2.setVectors( Teuchos::rcp( &xhat, false ), Teuchos::rcp( &b, false ) );
331 
332  // Factor and solve with matrix equilibration.
333  returnCode = solver2.factor();
334  testName = "Solve with matrix equilibration: factor() random A:";
335  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
336 
337  // Non-transpose solve
338  returnCode = solver2.solve();
339  testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
340  numberFailedTests += CompareVectors( *x2, xhat, tol );
341  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
342 
343 #ifdef HAVE_TEUCHOS_COMPLEX
344  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
346  solver2.setVectors( Teuchos::rcp( &xhatt, false ), Teuchos::rcp( &bct, false ) );
348  returnCode = solver2.solve();
349  testName = "Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
350  if (m == n)
351  numberFailedTests += CompareVectors( *x2t, xhatt, tol );
352  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
353 
354 #else
355  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
357  solver2.setVectors( Teuchos::rcp( &xhatt, false ), Teuchos::rcp( &bt, false ) );
359  returnCode = solver2.solve();
360  testName = "Solve with matrix equilibration: solve() random A (TRANS):";
361  if (m == n)
362  numberFailedTests += CompareVectors( *x2t, xhatt, tol );
363  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
364 #endif
365 
366  //
367  // If a test failed output the number of failed tests.
368  //
369  if(numberFailedTests > 0)
370  {
371  if (verbose) {
372  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
373  std::cout << "End Result: TEST FAILED" << std::endl;
374  return -1;
375  }
376  }
377  if(numberFailedTests == 0)
378  std::cout << "End Result: TEST PASSED" << std::endl;
379 
380  return 0;
381 }
382 
383 template<typename TYPE>
384 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
385 {
386  int result;
387  if(calculatedResult == expectedResult)
388  {
389  if(verbose) std::cout << testName << " successful." << std::endl;
390  result = 0;
391  }
392  else
393  {
394  if(verbose) std::cout << testName << " unsuccessful." << std::endl;
395  result = 1;
396  }
397  return result;
398 }
399 
400 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
401 {
402  int result;
403  if(expectedResult == 0)
404  {
405  if(returnCode == 0)
406  {
407  if(verbose) std::cout << testName << " test successful." << std::endl;
408  result = 0;
409  }
410  else
411  {
412  if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
413  result = 1;
414  }
415  }
416  else
417  {
418  if(returnCode != 0)
419  {
420  if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
421  result = 0;
422  }
423  else
424  {
425  if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
426  result = 1;
427  }
428  }
429  return result;
430 }
431 
432 template<typename TYPE>
433 TYPE GetRandom(TYPE Low, TYPE High)
434 {
435  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
436 }
437 
438 template<typename T>
439 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
440 {
441  T lowMag = Low.real();
442  T highMag = High.real();
443  T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
444  T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
445  return std::complex<T>( real, imag );
446 }
447 
448 template<>
449 int GetRandom(int Low, int High)
450 {
451  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
452 }
453 
454 template<>
455 double GetRandom(double Low, double High)
456 {
457  return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
458 }
459 
461 {
462  Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
463 
464  // Fill dense matrix with random entries.
465  for (int i=0; i<m; i++)
466  for (int j=0; j<n; j++) {
467  (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
468  }
469  return newmat;
470 }
471 
473 {
474  Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
475 
476  // Fill dense vector with random entries.
477  for (int i=0; i<n; i++) {
478  (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
479  }
480 
481  return newvec;
482 }
483 
484 /* Function: CompareVectors
485  Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
486 */
488  const SerialDenseVector<OTYPE,STYPE>& Vector2,
490 {
491  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
492 
493  SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
494  diff -= Vector2;
495 
496  MagnitudeType norm_diff = diff.normFrobenius();
497  MagnitudeType norm_v2 = Vector2.normFrobenius();
498  MagnitudeType temp = norm_diff;
499  if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
500  temp /= norm_v2;
501 
502  if (temp > Tolerance)
503  return 1;
504  else
505  return 0;
506 }
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
SerialDenseVector< OTYPE, STYPE > DVector
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
Sets the pointers for coefficient matrix.
Templated serial dense matrix class.
int multiplyQ(ETransp transq, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Left multiply the input matrix by the unitary matrix Q or its adjoint.
int solveR(ETransp transr, SerialDenseMatrix< OrdinalType, ScalarType > &C)
Solve input matrix on the left with the upper triangular matrix R or its adjoint. ...
int formQ()
Explicitly forms the unitary matrix Q.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
int multiply(ETransp transa, ETransp transb, ScalarType alpha, const SerialDenseMatrix< OrdinalType, ScalarType > &A, const SerialDenseMatrix< OrdinalType, ScalarType > &B, ScalarType beta)
Multiply A * B and add them to this; this = beta * this + alpha*A*B.
TYPE GetRandom(TYPE, TYPE)
Templated class for solving dense linear problems.
#define STYPE
This class creates and provides basic support for dense vectors of templated type as a specialization...
#define SCALARMAX
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()..
This structure defines some basic traits for a scalar field type.
A class for solving dense linear problems.
SerialDenseMatrix< OTYPE, STYPE > DMatrix
Teuchos::RCP< DVector > GetRandomVector(int n)
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Deprecated.
OrdinalType numRows() const
Returns the row dimension of this matrix.
int PrintTestResults(std::string, TYPE, TYPE, bool)
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Set all values in the matrix to a constant value.
#define OTYPE
Teuchos::RCP< DMatrix > GetRandomMatrix(int m, int n)
The Templated LAPACK Wrapper Class.
#define TEUCHOS_MAX(x, y)
std::string Teuchos_Version()
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor...
Templated serial dense vector class.
Defines basic traits for the scalar field type.
int CompareVectors(const SerialDenseVector< OTYPE, STYPE > &Vector1, const SerialDenseVector< OTYPE, STYPE > &Vector2, ScalarTraits< STYPE >::magnitudeType Tolerance)
Smart reference counting pointer class for automatic garbage collection.
int formR()
Explicitly forms the upper triangular matrix R.
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
int ReturnCodeCheck(std::string, int, int, bool)
OrdinalType numCols() const
Returns the column dimension of this matrix.
Reference-counted pointer class and non-member templated function implementations.
int main(int argc, char *argv[])
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
Sets the pointers for left and right hand side vector(s).
This class creates and provides basic support for dense rectangular matrix of templated type...