Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Tpetra_TsqrAdaptor_UQ_PCE.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Tpetra: Templated Linear Algebra Services Package
5// Copyright (2008) Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
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// ************************************************************************
38// @HEADER
39
40#ifndef TPETRA_TSQR_ADAPTOR_UQ_PCE_HPP
41#define TPETRA_TSQR_ADAPTOR_UQ_PCE_HPP
42
43#include <Tpetra_ConfigDefs.hpp> // HAVE_TPETRA_TSQR, etc.
44
45#ifdef HAVE_TPETRA_TSQR
46
48
49# include "Tsqr_NodeTsqrFactory.hpp" // create intranode TSQR object
50# include "Tsqr.hpp" // full (internode + intranode) TSQR
51# include "Tsqr_DistTsqr.hpp" // internode TSQR
52// Subclass of TSQR::MessengerBase, implemented using Teuchos
53// communicator template helper functions
54# include "Tsqr_TeuchosMessenger.hpp"
55# include "Tpetra_MultiVector.hpp"
56# include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
57# include <stdexcept>
58
59// Base TsqrAdator template we will specialize
60# include "Tpetra_TsqrAdaptor.hpp"
61
62namespace Tpetra {
63
70 template <class Storage, class LO, class GO, class Node>
71 class TsqrAdaptor< Tpetra::MultiVector< Sacado::UQ::PCE<Storage>,
72 LO, GO, Node > > :
73 public Teuchos::ParameterListAcceptorDefaultBase {
74 public:
75 typedef Tpetra::MultiVector< Sacado::UQ::PCE<Storage>, LO, GO, Node > MV;
76 typedef typename MV::scalar_type mp_scalar_type;
77
78 // For Sacado::UQ::PCE< Storage<Ordinal,Scalar,Device> > this is Scalar
79 typedef typename mp_scalar_type::scalar_type scalar_type;
80 typedef typename mp_scalar_type::ordinal_type mp_ordinal_type;
81 typedef typename MV::local_ordinal_type ordinal_type;
82 typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
83 typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
84
85 private:
86 using node_tsqr_factory_type =
87 TSQR::NodeTsqrFactory<scalar_type, ordinal_type,
88 typename MV::device_type>;
89 using node_tsqr_type = TSQR::NodeTsqr<ordinal_type, scalar_type>;
90 using dist_tsqr_type = TSQR::DistTsqr<ordinal_type, scalar_type>;
91 using tsqr_type = TSQR::Tsqr<ordinal_type, scalar_type>;
92
93 public:
100 TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
101 nodeTsqr_ (node_tsqr_factory_type::getNodeTsqr ()),
102 distTsqr_ (new dist_tsqr_type),
103 tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
104 ready_ (false)
105 {
106 setParameterList (plist);
107 }
108
110 TsqrAdaptor () :
111 nodeTsqr_ (new node_tsqr_factory_type::getNodeTsqr ()),
112 distTsqr_ (new dist_tsqr_type),
113 tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
114 ready_ (false)
115 {
116 setParameterList (Teuchos::null);
117 }
118
119 Teuchos::RCP<const Teuchos::ParameterList>
120 getValidParameters () const
121 {
122 using Teuchos::RCP;
123 using Teuchos::rcp;
124 using Teuchos::ParameterList;
125 using Teuchos::parameterList;
126
127 if (defaultParams_.is_null()) {
128 RCP<ParameterList> params = parameterList ("TSQR implementation");
129 params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
130 params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
131 defaultParams_ = params;
132 }
133 return defaultParams_;
134 }
135
136 void
137 setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
138 {
139 using Teuchos::ParameterList;
140 using Teuchos::parameterList;
141 using Teuchos::RCP;
142 using Teuchos::sublist;
143
144 RCP<ParameterList> params = plist.is_null() ?
145 parameterList (*getValidParameters ()) : plist;
146 nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
147 distTsqr_->setParameterList (sublist (params, "DistTsqr"));
148
149 this->setMyParamList (params);
150 }
151
173 void
174 factorExplicit (MV& A,
175 MV& Q,
176 dense_matrix_type& R,
177 const bool forceNonnegativeDiagonal=false)
178 {
179 prepareTsqr (Q); // Finish initializing TSQR.
180
181 ordinal_type numRows;
182 ordinal_type numCols;
183 ordinal_type LDA;
184 ordinal_type LDQ;
185 scalar_type* A_ptr;
186 scalar_type* Q_ptr;
187
188 getNonConstView (numRows, numCols, A_ptr, LDA, A);
189 getNonConstView (numRows, numCols, Q_ptr, LDQ, Q);
190 const bool contiguousCacheBlocks = false;
191 tsqr_->factorExplicitRaw (numRows, numCols, A_ptr, LDA,
192 Q_ptr, LDQ, R.values (), R.stride (),
193 contiguousCacheBlocks,
194 forceNonnegativeDiagonal);
195 }
196
227 int
228 revealRank (MV& Q,
229 dense_matrix_type& R,
230 const magnitude_type& tol)
231 {
232 prepareTsqr (Q); // Finish initializing TSQR.
233
234 // FIXME (mfh 18 Oct 2010) Check Teuchos::Comm<int> object in Q
235 // to make sure it is the same communicator as the one we are
236 // using in our dist_tsqr_type implementation.
237
238 ordinal_type numRows;
239 ordinal_type numCols;
240 scalar_type* Q_ptr;
241 ordinal_type LDQ;
242 getNonConstView (numRows, numCols, Q_ptr, LDQ, Q);
243 const bool contiguousCacheBlocks = false;
244 return tsqr_->revealRankRaw (numRows, numCols, Q_ptr, LDQ,
245 R.values (), R.stride (), tol,
246 contiguousCacheBlocks);
247 }
248
249 private:
251 Teuchos::RCP<node_tsqr_type> nodeTsqr_;
252
254 Teuchos::RCP<dist_tsqr_type> distTsqr_;
255
257 Teuchos::RCP<tsqr_type> tsqr_;
258
260 mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
261
263 bool ready_;
264
285 void
286 prepareTsqr (const MV& mv)
287 {
288 if (! ready_) {
289 prepareDistTsqr (mv);
290 ready_ = true;
291 }
292 }
293
300 void
301 prepareDistTsqr (const MV& mv)
302 {
303 using Teuchos::RCP;
304 using Teuchos::rcp_implicit_cast;
305 typedef TSQR::TeuchosMessenger<scalar_type> mess_type;
306 typedef TSQR::MessengerBase<scalar_type> base_mess_type;
307
308 RCP<const Teuchos::Comm<int> > comm = mv.getMap()->getComm();
309 RCP<mess_type> mess (new mess_type (comm));
310 RCP<base_mess_type> messBase = rcp_implicit_cast<base_mess_type> (mess);
311 distTsqr_->init (messBase);
312 }
313
320 static void
321 getNonConstView (ordinal_type& numRows,
322 ordinal_type& numCols,
323 scalar_type*& A_ptr,
324 ordinal_type& LDA,
325 const MV& A)
326 {
327 // FIXME (mfh 25 Oct 2010) We should be able to run TSQR even if
328 // storage of A uses nonconstant stride internally. We would
329 // have to copy and pack into a matrix with constant stride, and
330 // then unpack on exit. For now we choose just to raise an
331 // exception.
332 TEUCHOS_TEST_FOR_EXCEPTION
333 (! A.isConstantStride(), std::invalid_argument,
334 "TSQR does not currently support Tpetra::MultiVector "
335 "inputs that do not have constant stride.");
336
337 // FIXME (mfh 16 Jan 2016) When I got here, I found strides[0]
338 // instead of strides[1] for the stride. I don't think this is
339 // right. However, I don't know about these Stokhos scalar
340 // types so I'll just do what was here.
341 //
342 // STOKHOS' TYPES ARE NOT TESTED WITH TSQR REGULARLY SO IT IS
343 // POSSIBLE THAT THE ORIGINAL CODE WAS WRONG.
344
345 typedef typename MV::dual_view_type view_type;
346 typedef typename view_type::t_dev::array_type flat_array_type;
347
348 // Reinterpret the data as a longer array of the base scalar
349 // type. TSQR currently forbids MultiVector input with
350 // nonconstant stride, so we need not worry about that here.
351
352 flat_array_type flat_mv = A.getLocalViewDevice();
353
354 numRows = static_cast<ordinal_type> (flat_mv.extent(0));
355 numCols = static_cast<ordinal_type> (flat_mv.extent(1));
356 A_ptr = flat_mv.data ();
357
358 ordinal_type strides[2];
359 flat_mv.stride (strides);
360 LDA = strides[0];
361 }
362 };
363} // namespace Tpetra
364
365#endif // HAVE_TPETRA_TSQR
366
367#endif // TPETRA_TSQR_ADAPTOR_UQ_PCE_HPP
KokkosClassic::DefaultNode::DefaultNodeType Node