Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
334 
335  // Factor and solve with matrix equilibration.
336  returnCode = solver2.factor();
337  testName = "Solve with matrix equilibration: factor() random A:";
338  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
339 
340  // Non-transpose solve
341  returnCode = solver2.solve();
342  testName = "Solve with matrix equilibration: solve() random A (NO_TRANS):";
343  numberFailedTests += CompareVectors( *x2, xhat, tol );
344  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
345 
346 #ifdef HAVE_TEUCHOS_COMPLEX
347  // Conjugate tranpose solve (can be done after factorization, since factorization doesn't depend on this)
349  solver2.setVectors( Teuchos::rcp( &xhatt, false ), Teuchos::rcp( &bct, false ) );
351  returnCode = solver2.solve();
352  testName = "Solve with matrix equilibration: solve() random A (CONJ_TRANS):";
353  if (m == n)
354  numberFailedTests += CompareVectors( *x2t, xhatt, tol );
355  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
356 
357 #else
358  // Tranpose solve (can be done after factorization, since factorization doesn't depend on this)
360  solver2.setVectors( Teuchos::rcp( &xhatt, false ), Teuchos::rcp( &bt, false ) );
362  returnCode = solver2.solve();
363  testName = "Solve with matrix equilibration: solve() random A (TRANS):";
364  if (m == n)
365  numberFailedTests += CompareVectors( *x2t, xhatt, tol );
366  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
367 #endif
368 
369  // Non-transpose solve without call to factor.
370  xhat.putScalar( ScalarTraits<STYPE>::zero() );
371  solver2.setMatrix( A2bak );
372  solver2.setVectors( Teuchos::rcp( &xhat, false ), b2bak );
374  returnCode = solver2.solve();
375  testName = "Solve with matrix equilibration: solve() without factor() random A (NO_TRANS):";
376  numberFailedTests += CompareVectors( *x2, xhat, tol );
377  numberFailedTests += ReturnCodeCheck(testName, returnCode, 0, verbose);
378 
379  //
380  // If a test failed output the number of failed tests.
381  //
382  if(numberFailedTests > 0)
383  {
384  if (verbose) {
385  std::cout << "Number of failed tests: " << numberFailedTests << std::endl;
386  std::cout << "End Result: TEST FAILED" << std::endl;
387  return -1;
388  }
389  }
390  if(numberFailedTests == 0)
391  std::cout << "End Result: TEST PASSED" << std::endl;
392 
393  return 0;
394 }
395 
396 template<typename TYPE>
397 int PrintTestResults(std::string testName, TYPE calculatedResult, TYPE expectedResult, bool verbose)
398 {
399  int result;
400  if(calculatedResult == expectedResult)
401  {
402  if(verbose) std::cout << testName << " successful." << std::endl;
403  result = 0;
404  }
405  else
406  {
407  if(verbose) std::cout << testName << " unsuccessful." << std::endl;
408  result = 1;
409  }
410  return result;
411 }
412 
413 int ReturnCodeCheck(std::string testName, int returnCode, int expectedResult, bool verbose)
414 {
415  int result;
416  if(expectedResult == 0)
417  {
418  if(returnCode == 0)
419  {
420  if(verbose) std::cout << testName << " test successful." << std::endl;
421  result = 0;
422  }
423  else
424  {
425  if(verbose) std::cout << testName << " test unsuccessful. Return code was " << returnCode << "." << std::endl;
426  result = 1;
427  }
428  }
429  else
430  {
431  if(returnCode != 0)
432  {
433  if(verbose) std::cout << testName << " test successful -- failed as expected." << std::endl;
434  result = 0;
435  }
436  else
437  {
438  if(verbose) std::cout << testName << " test unsuccessful -- did not fail as expected. Return code was " << returnCode << "." << std::endl;
439  result = 1;
440  }
441  }
442  return result;
443 }
444 
445 template<typename TYPE>
446 TYPE GetRandom(TYPE Low, TYPE High)
447 {
448  return ((TYPE)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
449 }
450 
451 template<typename T>
452 std::complex<T> GetRandom( std::complex<T> Low, std::complex<T> High)
453 {
454  T lowMag = Low.real();
455  T highMag = High.real();
456  T real = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
457  T imag = (T)(((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (highMag - lowMag + ScalarTraits<T>::one())) + lowMag;
458  return std::complex<T>( real, imag );
459 }
460 
461 template<>
462 int GetRandom(int Low, int High)
463 {
464  return ((int)((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low);
465 }
466 
467 template<>
468 double GetRandom(double Low, double High)
469 {
470  return (((double)((1.0 * ScalarTraits<int>::random()) / RAND_MAX) * (High - Low + 1)) + Low + ScalarTraits<double>::random());
471 }
472 
474 {
475  Teuchos::RCP<DMatrix> newmat = Teuchos::rcp( new DMatrix(m,n) );
476 
477  // Fill dense matrix with random entries.
478  for (int i=0; i<m; i++)
479  for (int j=0; j<n; j++) {
480  (*newmat)(i,j) = GetRandom(-SCALARMAX, SCALARMAX);
481  }
482  return newmat;
483 }
484 
486 {
487  Teuchos::RCP<DVector> newvec = Teuchos::rcp( new DVector( n ) );
488 
489  // Fill dense vector with random entries.
490  for (int i=0; i<n; i++) {
491  (*newvec)(i) = GetRandom(-SCALARMAX, SCALARMAX);
492  }
493 
494  return newvec;
495 }
496 
497 /* Function: CompareVectors
498  Purpose: Compares the difference between two vectors using relative euclidean-norms, i.e. ||v_1-v_2||_2/||v_2||_2
499 */
501  const SerialDenseVector<OTYPE,STYPE>& Vector2,
503 {
504  typedef ScalarTraits<STYPE>::magnitudeType MagnitudeType;
505 
506  SerialDenseVector<OTYPE,STYPE> diff( Vector1 );
507  diff -= Vector2;
508 
509  MagnitudeType norm_diff = diff.normFrobenius();
510  MagnitudeType norm_v2 = Vector2.normFrobenius();
511  MagnitudeType temp = norm_diff;
512  if (norm_v2 != ScalarTraits<MagnitudeType>::zero())
513  temp /= norm_v2;
514 
515  if (temp > Tolerance)
516  return 1;
517  else
518  return 0;
519 }
int PrintTestResults(std::string, TYPE, TYPE, bool)
Non-member helper functions on the templated serial, dense matrix/vector classes. ...
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.
Teuchos::SerialDenseVector< int, std::complex< Real > > DVector
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.
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()..
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Returns the Frobenius-norm of the matrix.
This structure defines some basic traits for a scalar field type.
A class for solving dense linear problems.
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.
TYPE GetRandom(TYPE, TYPE)
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.
bool CompareVectors(TYPE1 *Vector1, TYPE2 *Vector2, OType Size, TYPE2 Tolerance)
#define TEUCHOS_MAX(x, y)
Teuchos::SerialDenseMatrix< int, std::complex< Real > > DMatrix
std::string Teuchos_Version()
int ReturnCodeCheck(std::string, int, int, bool)
int main(int argc, char *argv[])
OrdinalType numCols() const
Returns the column dimension of this matrix.
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.
Teuchos::RCP< DVector > GetRandomVector(int n)
Smart reference counting pointer class for automatic garbage collection.
int formR()
Explicitly forms the upper triangular matrix R.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
Reference-counted pointer class and non-member templated function implementations.
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).
OrdinalType numRows() const
Returns the row dimension of this matrix.
This class creates and provides basic support for dense rectangular matrix of templated type...