Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_MueLu_QR_Interface_def.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 #ifndef STOKHOS_MUELU_QR_INTERFACE_DEF_HPP
42 #define STOKHOS_MUELU_QR_INTERFACE_DEF_HPP
43 
44 #include "MueLu_QR_Interface_decl.hpp"
45 #include "MueLu_Exceptions.hpp"
46 
47 namespace MueLu {
48 
49 #if defined(HAVE_STOKHOS_MUELU)
50  //Specialization for polynomial chaos expansion (PCE) scalar types.
51  template <class Scalar, class Storage, class LocalOrdinal>
52  QR_Interface< Sacado::PCE::OrthogPoly<Scalar, Storage>, LocalOrdinal>::QR_Interface(const size_t NSDim) : workSize_(NSDim), NSDim_(NSDim), info_(0) {
53  tau_ = ArrayRCP<Scalar>(NSDim);
54  work_ = ArrayRCP<Scalar>(NSDim);
55  }
56 
57  template <class Scalar, class Storage, class LocalOrdinal>
58  void QR_Interface< Sacado::PCE::OrthogPoly<Scalar, Storage>, LocalOrdinal>::Compute(LocalOrdinal const &myAggSize, ArrayRCP<Sacado::PCE::OrthogPoly<Scalar, Storage> > &localQR) {
59  int pce_sz = localQR[0].size();
60  if (localQR.size()*pce_sz > localQR_.size())
61  localQR_.resize(localQR.size()*pce_sz);
62  //convert pce to pod scalar
63  for (int i=0; i<localQR.size(); ++i) {
64  for (int j=0; j<pce_sz; ++j) {
65  localQR_[i*pce_sz+j] = (localQR[i]).coeff(j);
66  }
67  }
68  lapack_.GEQRF( myAggSize*pce_sz, Teuchos::as<int>(NSDim_), localQR_.getRawPtr(), myAggSize*pce_sz,
69  tau_.getRawPtr(), work_.getRawPtr(), workSize_, &info_ );
70  if (info_ != 0) {
71  std::string msg = "QR_Interface: dgeqrf (LAPACK QR routine) returned error code " + Teuchos::toString(info_);
72  throw(Exceptions::RuntimeError(msg));
73  }
74  //promote POD scalar back to pce
75  for (int i=0; i<localQR.size(); ++i) {
76  localQR[i].copyForWrite();
77  if (localQR[i].size() < pce_sz)
78  localQR[i].reset(localQR[0].expansion());
79  for (int j=0; j<pce_sz; ++j)
80  localQR[i].fastAccessCoeff(j) = localQR_[i*pce_sz+j];
81  }
82  // LAPACK may have determined a better length for the work array. Returns it in work[0],
83  // so we cast to avoid compiler warnings. Taking a look at the NETLIB reference implementation
84  // CGEQRF (complex), work[0] is assigned an integer, so it's safe to take the magnitude.
85  // Scalar type might be complex, so take absolute value.
86  if ( std::abs(work_[0]) > workSize_) {
87  workSize_ = (int) std::abs(work_[0]);
88  work_ = ArrayRCP<Scalar>(workSize_);
89  }
90  } //Compute
91 
92  template <class Scalar, class Storage, class LocalOrdinal>
93  void QR_Interface< Sacado::PCE::OrthogPoly<Scalar, Storage>, LocalOrdinal>::ExtractQ(LocalOrdinal const &myAggSize, ArrayRCP<Sacado::PCE::OrthogPoly<Scalar, Storage> > &localQR) {
94  //call nonmember function (perhaps specialized)
95  //Note: localQR_ already contains the proper data because of prior call to Compute, so there is no need to resize or copy.
96  // If Compute is called twice in a row, all bets are off.
97  int pce_sz = localQR[0].size();
98  LapackQR( lapack_, myAggSize*pce_sz, Teuchos::as<int>(NSDim_), localQR_, tau_, work_, workSize_, info_ );
99  if (info_ != 0) {
100  std::string msg = "QR_Interface: dorgqr (LAPACK auxiliary QR routine) returned error code " + Teuchos::toString(info_);
101  throw(Exceptions::RuntimeError(msg));
102  }
103  //promote POD scalar back to pce
104  for (int i=0; i<localQR.size(); ++i) {
105  localQR[i].copyForWrite();
106  if (localQR[i].size() < pce_sz)
107  localQR[i].reset(localQR[0].expansion());
108  for (int j=0; j<pce_sz; ++j)
109  localQR[i].fastAccessCoeff(j) = localQR_[i*pce_sz+j];
110  }
111 
112  // LAPACK may have determined a better length for the work array. Returns it in work[0],
113  // so we cast to avoid compiler warnings.
114  // Scalar type might be complex, so take absolute value.
115  if ( std::abs(work_[0]) > workSize_) {
116  workSize_ = (int) std::abs(work_[0]);
117  work_ = ArrayRCP<Scalar>(workSize_);
118  }
119  } //ExtractQ
120 #endif //if defined(HAVE_STOKHOS_MUELU)
121 
122 
123 } //namespace MueLu
124 
125 #endif // STOKHOS_MUELU_QR_INTERFACE_DEF_HPP
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr expr expr expr j