Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Epetra_TsqrAdaptor.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 EPETRA_TSQRADAPTOR_HPP
41#define EPETRA_TSQRADAPTOR_HPP
42
54
55#include "Tpetra_ConfigDefs.hpp"
56
57#if defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
58
59#include "Tsqr_NodeTsqrFactory.hpp" // create intranode TSQR object
60#include "Tsqr.hpp" // full (internode + intranode) TSQR
61#include "Tsqr_DistTsqr.hpp" // internode TSQR
62#include "Epetra_Comm.h"
63// Subclass of TSQR::MessengerBase, implemented using Teuchos
64// communicator template helper functions
66#include "Epetra_MultiVector.h"
67#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
68#include <stdexcept>
69
70namespace Epetra {
71
96 class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
97 public:
98 typedef Epetra_MultiVector MV;
99
106 typedef double scalar_type;
107
114 typedef int ordinal_type;
115
121 using device_type =
122 Kokkos::Device<Kokkos::DefaultHostExecutionSpace,
123 Kokkos::HostSpace>;
124
133 using dense_matrix_type =
134 Teuchos::SerialDenseMatrix<ordinal_type, scalar_type>;
135
141 using magnitude_type = double;
142
143 private:
144 using matview_type = TSQR::MatView<ordinal_type, scalar_type>;
145 using node_tsqr_factory_type =
146 TSQR::NodeTsqrFactory<scalar_type, ordinal_type, device_type>;
147 // Don't need a "typename" here, because there are no template
148 // parameters involved in the type definition.
149 using node_tsqr_type = TSQR::NodeTsqr<ordinal_type, scalar_type>;
150 using dist_tsqr_type = TSQR::DistTsqr<ordinal_type, scalar_type>;
151 using tsqr_type = TSQR::Tsqr<ordinal_type, scalar_type>;
152
153 public:
160 TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& plist) :
161 nodeTsqr_ (node_tsqr_factory_type::getNodeTsqr ()),
162 distTsqr_ (new dist_tsqr_type),
163 tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
164 ready_ (false)
165 {
166 setParameterList (plist);
167 }
168
170 TsqrAdaptor () :
171 nodeTsqr_ (node_tsqr_factory_type::getNodeTsqr ()),
172 distTsqr_ (new dist_tsqr_type),
173 tsqr_ (new tsqr_type (nodeTsqr_, distTsqr_)),
174 ready_ (false)
175 {
176 setParameterList (Teuchos::null);
177 }
178
179 Teuchos::RCP<const Teuchos::ParameterList>
180 getValidParameters () const
181 {
182 using Teuchos::RCP;
183 using Teuchos::rcp;
184 using Teuchos::ParameterList;
185 using Teuchos::parameterList;
186
187 if (defaultParams_.is_null()) {
188 RCP<ParameterList> params = parameterList ("TSQR implementation");
189 params->set ("NodeTsqr", *(nodeTsqr_->getValidParameters ()));
190 params->set ("DistTsqr", *(distTsqr_->getValidParameters ()));
191 defaultParams_ = params;
192 }
193 return defaultParams_;
194 }
195
196 void
197 setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
198 {
199 using Teuchos::ParameterList;
200 using Teuchos::parameterList;
201 using Teuchos::RCP;
202 using Teuchos::sublist;
203
204 RCP<ParameterList> params = plist.is_null() ?
205 parameterList (*getValidParameters ()) : plist;
206 nodeTsqr_->setParameterList (sublist (params, "NodeTsqr"));
207 distTsqr_->setParameterList (sublist (params, "DistTsqr"));
208
209 this->setMyParamList (params);
210 }
211
233 void
234 factorExplicit (MV& A,
235 MV& Q,
236 dense_matrix_type& R,
237 const bool forceNonnegativeDiagonal=false)
238 {
239 prepareTsqr (Q); // Finish initializing TSQR.
240
241 scalar_type* const A_ptr = A.Values ();
242 scalar_type* const Q_ptr = Q.Values ();
243 scalar_type* const R_ptr = R.values ();
244 const ordinal_type numRows = A.MyLength ();
245 const ordinal_type numCols = A.NumVectors ();
246 const ordinal_type lda = A.Stride ();
247 const ordinal_type ldq = Q.Stride ();
248 const ordinal_type ldr = R.stride ();
249
250 const bool contiguousCacheBlocks = false;
251 tsqr_->factorExplicitRaw (numRows, numCols, A_ptr, lda,
252 Q_ptr, ldq, R_ptr, ldr,
253 contiguousCacheBlocks,
254 forceNonnegativeDiagonal);
255 }
256
287 int
288 revealRank (MV& Q,
289 dense_matrix_type& R,
290 const magnitude_type& tol)
291 {
292 TEUCHOS_TEST_FOR_EXCEPTION
293 (! Q.ConstantStride (), std::invalid_argument, "TsqrAdaptor::"
294 "revealRank: Input MultiVector Q must have constant stride.");
295 prepareTsqr (Q); // Finish initializing TSQR.
296 // FIXME (mfh 25 Oct 2010) Check Epetra_Comm object in Q to make
297 // sure it is the same communicator as the one we are using in
298 // our dist_tsqr_type implementation.
299 return tsqr_->revealRankRaw (Q.MyLength (), Q.NumVectors (),
300 Q.Values (), Q.Stride (),
301 R.values (), R.stride (), tol, false);
302 }
303
304 private:
306 Teuchos::RCP<node_tsqr_type> nodeTsqr_;
307
309 Teuchos::RCP<dist_tsqr_type> distTsqr_;
310
312 Teuchos::RCP<tsqr_type> tsqr_;
313
315 mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
316
318 bool ready_;
319
338 void
339 prepareTsqr (const MV& mv)
340 {
341 if (! ready_) {
342 prepareDistTsqr (mv);
343 ready_ = true;
344 }
345 }
346
353 void
354 prepareDistTsqr (const MV& mv)
355 {
356 using Teuchos::RCP;
357 using Teuchos::rcp;
358 using TSQR::Epetra::makeTsqrMessenger;
359 typedef TSQR::MessengerBase<scalar_type> base_mess_type;
360
361 // If mv falls out of scope, its Epetra_Comm may become invalid.
362 // Thus, we clone the input Epetra_Comm, so that the messenger
363 // owns the object.
364 RCP<const Epetra_Comm> comm = rcp (mv.Comm().Clone());
365 RCP<base_mess_type> messBase = makeTsqrMessenger<scalar_type> (comm);
366 distTsqr_->init (messBase);
367 }
368 };
369
370} // namespace Epetra
371
372#endif // defined(HAVE_TPETRA_EPETRA) && defined(HAVE_TPETRA_TSQR)
373
374#endif // EPETRA_TSQRADAPTOR_HPP
375
Function for wrapping Epetra_Comm in a communicator wrapper that TSQR can use.