Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Teuchos_SerialQRDenseSolver_MP_Vector.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) 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 Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef _TEUCHOS_SERIALQRDENSESOLVER_MP_VECTOR_HPP_
43#define _TEUCHOS_SERIALQRDENSESOLVER_MP_VECTOR_HPP_
48#include <new>
49#include "Teuchos_SerialQRDenseSolver.hpp"
50
51#include "Sacado_MP_Vector.hpp"
52
57namespace Teuchos {
58
59 template<typename OrdinalType, typename Storage>
60 class SerialQRDenseSolver<OrdinalType, Sacado::MP::Vector<Storage> >
61 : public CompObject,
62 public Object
63 {
64
65 public:
66
68 typedef typename ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
69
71
72
74
78
80
81
83
85 int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);
86
88
91 int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
92 const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);
94
96
97
99
101 void factorWithEquilibration(bool flag) {
102 base_QR_.factorWithEquilibration(flag);
103 }
104
106 void solveWithTranspose(bool flag) {
107 base_QR_.solveWithTranspose(flag);
108 }
109
111 void solveWithTransposeFlag(Teuchos::ETransp trans) {
112 base_QR_.solveWithTranspose(trans);
113 }
114
116
118
119
121
124 int factor() { return base_QR_.factor(); }
125
127
130 int solve() { return base_QR_.solve(); }
131
133
137 return base_QR_.computeEquilibrateScaling();
138 }
139
141
145 int equilibrateMatrix() { return base_QR_.equilibrateMatrix(); }
146
148
152 int equilibrateRHS() { return base_QR_.equilibrateRHS(); }
153
155
159 int unequilibrateLHS() { return base_QR_.unequilibrateLHS(); }
160
162
165 int formQ();
166
168
171 int formR();
172
174
177 int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);
178
180
183 int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);
185
187
188
190 bool transpose() { return base_QR_.transpose(); }
191
193 bool factored() { return base_QR_.factored(); }
194
196 bool equilibratedA() { return base_QR_.equilibratedA(); }
197
199 bool equilibratedB() { return base_QR_.equilibratedB(); }
200
202 bool shouldEquilibrate() { return base_QR_.shouldEquilibrate(); }
203
205 bool solved() { return base_QR_.solved(); }
206
208 bool formedQ() { return base_QR_.formedQ(); }
209
211 bool formedR() { return base_QR_.formedR();}
212
214
216
217
219 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}
220
222 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix() const {return(Factor_);}
223
225 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getQ() const {return(FactorQ_);}
226
228 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getR() const {return(FactorR_);}
229
231 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}
232
234 RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}
235
237 OrdinalType numRows() const {return(M_);}
238
240 OrdinalType numCols() const {return(N_);}
241
243 std::vector<ScalarType> tau() const { return base_QR_.tau(); }
244
246 MagnitudeType ANORM() const { return base_QR_.ANORM(); }
247
249
251
252
253 void Print(std::ostream& os) const;
255 protected:
256
257 typedef typename ScalarType::value_type BaseScalarType;
259 typedef SerialDenseMatrix<OrdinalType, BaseScalarType> BaseMatrixType;
260 typedef SerialDenseMatrix<OrdinalType, ScalarType> MatrixType;
261
262 void resetMatrix();
263 void resetVectors();
264 RCP<BaseMatrixType> createBaseMatrix(const RCP<MatrixType>& mat) const;
265 RCP<MatrixType> createMatrix(const RCP<BaseMatrixType>& base_mat) const;
267
268 OrdinalType M_;
269 OrdinalType N_;
270 OrdinalType SacadoSize_;
271
272 RCP<MatrixType> Matrix_;
273 RCP<MatrixType> LHS_;
274 RCP<MatrixType> RHS_;
275 RCP<MatrixType> Factor_;
276 RCP<MatrixType> FactorQ_;
277 RCP<MatrixType> FactorR_;
278
279 RCP<BaseMatrixType> Base_Matrix_;
280 RCP<BaseMatrixType> Base_LHS_;
281 RCP<BaseMatrixType> Base_RHS_;
282 RCP<BaseMatrixType> Base_Factor_;
283 RCP<BaseMatrixType> Base_FactorQ_;
284 RCP<BaseMatrixType> Base_FactorR_;
285
286 private:
287
288 // SerialQRDenseSolver copy constructor (put here because we don't want user access)
291
292 };
293
294 namespace details {
295
296 // Helper traits class for converting between arrays of
297 // Sacado::MP::Vector and its corresponding value type
298 template <typename Storage> struct MPVectorArrayHelper;
299
300 template <typename Ordinal, typename Value, typename Device>
301 struct MPVectorArrayHelper< Stokhos::DynamicStorage<Ordinal,Value,Device> >
302 {
305
306 static Value* getValueArray(Vector* v, Ordinal len) {
307 if (len == 0)
308 return 0;
309 return v[0].coeff(); // Assume data is contiguous
310 }
311
312 static Vector* getVectorArray(Value *v, Ordinal len, Ordinal vector_size){
313 if (len == 0)
314 return 0;
315 Vector* v_vec = static_cast<Vector*>(operator new(len*sizeof(Vector)));
316 for (Ordinal i=0; i<len; ++i)
317 new (v_vec+i) Vector(vector_size, v+i*vector_size, false);
318 return v_vec;
319 }
320
321 };
322
323 template <typename Ordinal, typename Value, typename Device, int Num>
324 struct MPVectorArrayHelper< Stokhos::StaticFixedStorage<Ordinal,Value,Num,Device> > {
327
328 static Value* getValueArray(Vector* v, Ordinal len)
329 {
330 return reinterpret_cast<Value*>(v);
331 }
332
333 static Vector* getVectorArray(Value *v, Ordinal len, Ordinal vector_size)
334 {
335 return reinterpret_cast<Vector*>(v);
336 }
337 };
338
339 }
340
341 // Helper traits to distinguish work arrays for real and complex-valued datatypes.
342 using namespace details;
343
344//=============================================================================
345
346template<typename OrdinalType, typename Storage>
348SerialQRDenseSolver()
349 : CompObject(),
350 Object("Teuchos::SerialQRDenseSolver"),
351 M_(0),
352 N_(0)
353{
354 resetMatrix();
355}
356
357//=============================================================================
358
359template<typename OrdinalType, typename Storage>
361resetVectors()
362{
363 LHS_ = Teuchos::null;
364 RHS_ = Teuchos::null;
365}
366//=============================================================================
367
368template<typename OrdinalType, typename Storage>
370resetMatrix()
371{
372 resetVectors();
373 M_ = 0;
374 N_ = 0;
375}
376//=============================================================================
377
378template<typename OrdinalType, typename Storage>
379RCP<SerialDenseMatrix<OrdinalType, typename Sacado::MP::Vector<Storage>::value_type> >
381createBaseMatrix(
382 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& mat) const
383{
384 BaseScalarType* base_ptr =
386 mat->values(), mat->stride()*mat->numCols());
387 RCP<BaseMatrixType> base_mat =
388 rcp(new BaseMatrixType(Teuchos::View,
389 base_ptr,
390 mat->stride()*SacadoSize_,
391 mat->numRows()*SacadoSize_,
392 mat->numCols()));
393 return base_mat;
394}
395//=============================================================================
396
397template<typename OrdinalType, typename Storage>
398RCP<SerialDenseMatrix<OrdinalType, Sacado::MP::Vector<Storage> > >
400createMatrix(
401 const RCP<SerialDenseMatrix<OrdinalType,BaseScalarType> >& base_mat) const
402{
403 ScalarType* ptr =
405 base_mat->values(), base_mat->stride()*base_mat->numCols(), SacadoSize_);
406 RCP<MatrixType> mat =
407 rcp(new MatrixType(Teuchos::View,
408 ptr,
409 base_mat->stride()/SacadoSize_,
410 base_mat->numRows()/SacadoSize_,
411 base_mat->numCols()));
412 return mat;
413}
414//=============================================================================
415
416template<typename OrdinalType, typename Storage>
418setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
419{
420 TEUCHOS_TEST_FOR_EXCEPTION(
421 A->numRows() < A->numCols(), std::invalid_argument,
422 "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
423
424 resetMatrix();
425 Matrix_ = A;
426 Factor_ = A;
427 FactorQ_ = A;
428 FactorR_ = A;
429 M_ = A->numRows();
430 N_ = A->numCols();
431 if (Storage::is_static)
432 SacadoSize_ = Storage::static_size;
433 else if (M_ > 0 && N_ > 0)
434 SacadoSize_ = (*A)(0,0).size();
435 else
436 SacadoSize_ = 0;
437
438 return base_QR_.setMatrix( createBaseMatrix(A) );
439}
440//=============================================================================
441
442template<typename OrdinalType, typename Storage>
444setVectors(
445 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
446 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
447{
448 // Check that these new vectors are consistent
449 TEUCHOS_TEST_FOR_EXCEPTION(
450 B->numCols() != X->numCols(), std::invalid_argument,
451 "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
452 TEUCHOS_TEST_FOR_EXCEPTION(
453 B->values()==0, std::invalid_argument,
454 "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
455 TEUCHOS_TEST_FOR_EXCEPTION(
456 X->values()==0, std::invalid_argument,
457 "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
458 TEUCHOS_TEST_FOR_EXCEPTION(
459 B->stride()<1, std::invalid_argument,
460 "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
461 TEUCHOS_TEST_FOR_EXCEPTION(
462 X->stride()<1, std::invalid_argument,
463 "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
464
465 resetVectors();
466 LHS_ = X;
467 RHS_ = B;
468
469 return base_QR_.setVectors( createBaseMatrix(X), createBaseMatrix(B) );
470}
471
472//=============================================================================
473
474template<typename OrdinalType, typename Storage>
476formQ() {
477 int ret = base_QR_.formQ();
478 FactorQ_ = createMatrix( base_QR_.getQ() );
479 return(ret);
480}
481
482//=============================================================================
483
484template<typename OrdinalType, typename Storage>
486formR() {
487 int ret = base_QR_.formR();
488 FactorR_ = createMatrix( base_QR_.getR() );
489 Factor_ = createMatrix( base_QR_.getFactoredMatrix() );
490 return(ret);
491}
492
493//=============================================================================
494
495template<typename OrdinalType, typename Storage>
497multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
498{
499 return base_QR_.multiplyQ(transq, createBaseMatrix(C));
500}
501
502//=============================================================================
503
504template<typename OrdinalType, typename Storage>
506solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
507{
508 return base_QR_.solveR(transr, createBaseMatrix(C));
509}
510
511//=============================================================================
512
513template<typename OrdinalType, typename Storage>
515Print(std::ostream& os) const {
516
517 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
518 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
519 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
520 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
521
522}
523
524} // namespace Teuchos
525
526#endif /* _TEUCHOS_SERIALQRDENSESOLVER_MP_VECTOR_HPP_ */
Statically allocated storage class.
SerialQRDenseSolver & operator=(const SerialQRDenseSolver &Source)
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix,...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed).
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
bool transpose()
Returns true if adjoint of this matrix has and will be used.
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
Specialization for Sacado::UQ::PCE< Storage<...> >
Top-level namespace for Stokhos classes and functions.