Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Teuchos_SerialQRDenseSolver_UQ_PCE.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_UQ_PCE_HPP_
43#define _TEUCHOS_SERIALQRDENSESOLVER_UQ_PCE_HPP_
48#include <new>
49#include "Teuchos_SerialQRDenseSolver.hpp"
50
51#include "Sacado_UQ_PCE.hpp"
52
57namespace Teuchos {
58
59 template<typename OrdinalType, typename Storage>
60 class SerialQRDenseSolver<OrdinalType, Sacado::UQ::PCE<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::UQ::PCE and its corresponding value type
298 template <typename Storage> struct PCEArrayHelper;
299
300 template <typename Ordinal, typename Value, typename Device>
301 struct PCEArrayHelper< Stokhos::DynamicStorage<Ordinal,Value,Device> >
302 {
305 typedef typename pce_type::cijk_type cijk_type;
306
307 static Value* getValueArray(pce_type* v, Ordinal len) {
308 if (len == 0)
309 return 0;
310 return v[0].coeff(); // Assume data is contiguous
311 }
312
313 static pce_type* getPCEArray(Value *v, Ordinal len, Ordinal vector_size){
314 if (len == 0)
315 return 0;
316 pce_type* v_pce =
317 static_cast<pce_type*>(operator new(len*sizeof(pce_type)));
318 cijk_type cijk = Kokkos::getGlobalCijkTensor<cijk_type>();
319 for (Ordinal i=0; i<len; ++i)
320 new (v_pce+i) pce_type(cijk, vector_size, v+i*vector_size, false);
321 return v_pce;
322 }
323
324 };
325
326 }
327
328 // Helper traits to distinguish work arrays for real and complex-valued datatypes.
329 using namespace details;
330
331//=============================================================================
332
333template<typename OrdinalType, typename Storage>
335SerialQRDenseSolver()
336 : CompObject(),
337 Object("Teuchos::SerialQRDenseSolver"),
338 M_(0),
339 N_(0)
340{
341 resetMatrix();
342}
343
344//=============================================================================
345
346template<typename OrdinalType, typename Storage>
348resetVectors()
349{
350 LHS_ = Teuchos::null;
351 RHS_ = Teuchos::null;
352}
353//=============================================================================
354
355template<typename OrdinalType, typename Storage>
357resetMatrix()
358{
359 resetVectors();
360 M_ = 0;
361 N_ = 0;
362}
363//=============================================================================
364
365template<typename OrdinalType, typename Storage>
366RCP<SerialDenseMatrix<OrdinalType, typename Sacado::UQ::PCE<Storage>::value_type> >
368createBaseMatrix(
369 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& mat) const
370{
371 BaseScalarType* base_ptr =
373 mat->values(), mat->stride()*mat->numCols());
374 RCP<BaseMatrixType> base_mat =
375 rcp(new BaseMatrixType(Teuchos::View,
376 base_ptr,
377 mat->stride()*SacadoSize_,
378 mat->numRows()*SacadoSize_,
379 mat->numCols()));
380 return base_mat;
381}
382//=============================================================================
383
384template<typename OrdinalType, typename Storage>
385RCP<SerialDenseMatrix<OrdinalType, Sacado::UQ::PCE<Storage> > >
387createMatrix(
388 const RCP<SerialDenseMatrix<OrdinalType,BaseScalarType> >& base_mat) const
389{
390 ScalarType* ptr =
392 base_mat->values(), base_mat->stride()*base_mat->numCols(), SacadoSize_);
393 RCP<MatrixType> mat =
394 rcp(new MatrixType(Teuchos::View,
395 ptr,
396 base_mat->stride()/SacadoSize_,
397 base_mat->numRows()/SacadoSize_,
398 base_mat->numCols()));
399 return mat;
400}
401//=============================================================================
402
403template<typename OrdinalType, typename Storage>
405setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
406{
407 TEUCHOS_TEST_FOR_EXCEPTION(
408 A->numRows() < A->numCols(), std::invalid_argument,
409 "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");
410
411 resetMatrix();
412 Matrix_ = A;
413 Factor_ = A;
414 FactorQ_ = A;
415 FactorR_ = A;
416 M_ = A->numRows();
417 N_ = A->numCols();
418 if (Storage::is_static)
419 SacadoSize_ = Storage::static_size;
420 else if (M_ > 0 && N_ > 0)
421 SacadoSize_ = (*A)(0,0).size();
422 else
423 SacadoSize_ = 0;
424
425 return base_QR_.setMatrix( createBaseMatrix(A) );
426}
427//=============================================================================
428
429template<typename OrdinalType, typename Storage>
431setVectors(
432 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
433 const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
434{
435 // Check that these new vectors are consistent
436 TEUCHOS_TEST_FOR_EXCEPTION(
437 B->numCols() != X->numCols(), std::invalid_argument,
438 "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
439 TEUCHOS_TEST_FOR_EXCEPTION(
440 B->values()==0, std::invalid_argument,
441 "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
442 TEUCHOS_TEST_FOR_EXCEPTION(
443 X->values()==0, std::invalid_argument,
444 "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
445 TEUCHOS_TEST_FOR_EXCEPTION(
446 B->stride()<1, std::invalid_argument,
447 "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
448 TEUCHOS_TEST_FOR_EXCEPTION(
449 X->stride()<1, std::invalid_argument,
450 "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");
451
452 resetVectors();
453 LHS_ = X;
454 RHS_ = B;
455
456 return base_QR_.setVectors( createBaseMatrix(X), createBaseMatrix(B) );
457}
458
459//=============================================================================
460
461template<typename OrdinalType, typename Storage>
463formQ() {
464 int ret = base_QR_.formQ();
465 FactorQ_ = createMatrix( base_QR_.getQ() );
466 return(ret);
467}
468
469//=============================================================================
470
471template<typename OrdinalType, typename Storage>
473formR() {
474 int ret = base_QR_.formR();
475 FactorR_ = createMatrix( base_QR_.getR() );
476 Factor_ = createMatrix( base_QR_.getFactoredMatrix() );
477 return(ret);
478}
479
480//=============================================================================
481
482template<typename OrdinalType, typename Storage>
484multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
485{
486 return base_QR_.multiplyQ(transq, createBaseMatrix(C));
487}
488
489//=============================================================================
490
491template<typename OrdinalType, typename Storage>
493solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
494{
495 return base_QR_.solveR(transr, createBaseMatrix(C));
496}
497
498//=============================================================================
499
500template<typename OrdinalType, typename Storage>
502Print(std::ostream& os) const {
503
504 if (Matrix_!=Teuchos::null) os << "Solver Matrix" << std::endl << *Matrix_ << std::endl;
505 if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
506 if (LHS_ !=Teuchos::null) os << "Solver LHS" << std::endl << *LHS_ << std::endl;
507 if (RHS_ !=Teuchos::null) os << "Solver RHS" << std::endl << *RHS_ << std::endl;
508
509}
510
511} // namespace Teuchos
512
513#endif /* _TEUCHOS_SERIALQRDENSESOLVER_UQ_PCE_HPP_ */
SerialQRDenseSolver & operator=(const SerialQRDenseSolver &Source)
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getFactoredMatrix() const
Returns pointer to factored matrix (assuming factorization has been performed).
bool equilibratedA()
Returns true if factor is equilibrated (factor available via getFactoredMatrix()).
int unequilibrateLHS()
Unscales the solution vectors if equilibration was used to solve the system.
bool equilibratedB()
Returns true if RHS is equilibrated (RHS available via getRHS()).
bool factored()
Returns true if matrix is factored (factor available via getFactoredMatrix()).
bool solved()
Returns true if the current set of vectors has been solved.
bool shouldEquilibrate()
Returns true if the LAPACK general rules for equilibration suggest you should equilibrate the system.
int factor()
Computes the in-place QR factorization of the matrix using the LAPACK routine _GETRF or the Eigen cla...
void factorWithEquilibration(bool flag)
Causes equilibration to be called just before the matrix factorization as part of the call to factor.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getQ() const
Returns pointer to Q (assuming factorization has been performed).
void solveWithTranspose(bool flag)
If flag is true, causes all subsequent function calls to work with the adjoint of this matrix,...
std::vector< ScalarType > tau() const
Returns pointer to pivot vector (if factorization has been computed), zero otherwise.
bool transpose()
Returns true if adjoint of this matrix has and will be used.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getLHS() const
Returns pointer to current LHS.
MagnitudeType ANORM() const
Returns the absolute value of the largest element of this matrix (returns -1 if not yet computed).
void solveWithTransposeFlag(Teuchos::ETransp trans)
All subsequent function calls will work with the transpose-type set by this method (Teuchos::NO_TRANS...
int solve()
Computes the solution X to AX = B for the this matrix and the B provided to SetVectors()....
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getRHS() const
Returns pointer to current RHS.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getMatrix() const
Returns pointer to current matrix.
RCP< SerialDenseMatrix< OrdinalType, ScalarType > > getR() const
Returns pointer to R (assuming factorization has been performed).
Specialization for Sacado::UQ::PCE< Storage<...> >
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
Top-level namespace for Stokhos classes and functions.