Stratimikos Version of the Day
Loading...
Searching...
No Matches
Thyra_TsqrAdaptor.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Stratimikos: Thyra-based strategies for linear solvers
5// Copyright (2006) 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 Roscoe A. Bartlett (rabartl@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef __Thyra_TsqrAdaptor_hpp
43#define __Thyra_TsqrAdaptor_hpp
44
45#include "BelosConfigDefs.hpp"
46
47// BelosThyraAdapter.hpp only includes this file if HAVE_BELOS_TSQR is
48// defined. Thus, it's OK to include TSQR header files here.
49
50#include "Thyra_MultiVectorBase.hpp"
51#include "Thyra_SpmdVectorSpaceBase.hpp"
52
53#ifdef HAVE_MPI
54# include "Teuchos_DefaultMpiComm.hpp"
55#endif // HAVE_MPI
56#include "Teuchos_DefaultSerialComm.hpp"
57#include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
58
59#include <stdexcept>
60
61
62namespace Thyra {
63
85 template<class Scalar>
86 class TsqrAdaptor : public Teuchos::ParameterListAcceptorDefaultBase {
87 public:
88 typedef Thyra::MultiVectorBase<Scalar> MV;
89 typedef Scalar scalar_type;
90 typedef int ordinal_type; // MultiVectorBase really does use int for this
91 typedef Teuchos::SerialDenseMatrix<ordinal_type, scalar_type> dense_matrix_type;
92 typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
93
100 TsqrAdaptor (const Teuchos::RCP<Teuchos::ParameterList>& /* plist */)
101 {
102 throw std::logic_error ("Thyra adaptor for TSQR not implemented");
103 }
104
107 {
108 throw std::logic_error ("Thyra adaptor for TSQR not implemented");
109 }
110
111 Teuchos::RCP<const Teuchos::ParameterList>
112 getValidParameters () const
113 {
114 throw std::logic_error ("Thyra adaptor for TSQR not implemented");
115 }
116
117 void
118 setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& /* plist */)
119 {
120 throw std::logic_error ("Thyra adaptor for TSQR not implemented");
121 }
122
144 void
145 factorExplicit (MV& /* A */,
146 MV& /* Q */,
147 dense_matrix_type& /* R */,
148 const bool /* forceNonnegativeDiagonal */ = false)
149 {
150 throw std::logic_error ("Thyra adaptor for TSQR not implemented");
151 }
152
183 int
184 revealRank (MV& /* Q */,
185 dense_matrix_type& /* R */,
186 const magnitude_type& /* tol */)
187 {
188 throw std::logic_error ("Thyra adaptor for TSQR not implemented");
189 }
190
191 private:
202 static Teuchos::RCP<const Teuchos::Comm<int> >
203 getComm (const MV& X)
204 {
205 using Teuchos::RCP;
206 using Teuchos::rcp;
207 using Teuchos::rcp_dynamic_cast;
208 using Teuchos::rcp_implicit_cast;
209 typedef Thyra::VectorSpaceBase<Scalar> space_base_type;
210 typedef Thyra::SpmdVectorSpaceBase<Scalar> space_type;
211
212 // Thyra stores the communicator in the "vector space," but only
213 // if that vector space is an SpmdVectorSpaceBase.
214 RCP<const space_base_type> rangeBase = X.range ();
215 TEUCHOS_TEST_FOR_EXCEPTION(rangeBase.is_null (), std::runtime_error, "X.range() is null.");
216 RCP<const space_type> range = rcp_dynamic_cast<const space_type> (rangeBase);
217 TEUCHOS_TEST_FOR_EXCEPTION(range.is_null (), std::runtime_error, "X.range() is not an SpmdVectorSpaceBase.");
218
219 // Thyra annoyingly uses a (possibly) different template
220 // parameter for its Teuchos::Comm than everybody else. The
221 // least hackish way to work around this is to convert the Comm
222 // to one of two subclasses (MpiComm or SerialComm). If it's an
223 // MpiComm, we can extract the RCP<const OpaqueWrapper<MPI_Comm>
224 // > and make a new MpiComm<int> from it. If it's a SerialComm,
225 // just create a new SerialComm<int>. If it's neither of those,
226 // then I have no idea what to do. Note that MpiComm is only
227 // defined if HAVE_MPI is defined.
228 RCP<const Teuchos::Comm<Thyra::Ordinal> > thyraComm = range->getComm ();
229#ifdef HAVE_MPI
230 RCP<const Teuchos::MpiComm<Thyra::Ordinal> > thyraMpiComm =
231 rcp_dynamic_cast<const Teuchos::MpiComm<Thyra::Ordinal> > (thyraComm);
232 if (thyraMpiComm.is_null ()) {
233 RCP<const Teuchos::SerialComm<Thyra::Ordinal> > thyraSerialComm =
234 rcp_dynamic_cast<const Teuchos::SerialComm<Thyra::Ordinal> > (thyraComm);
235 TEUCHOS_TEST_FOR_EXCEPTION(
236 thyraSerialComm.is_null (), std::runtime_error,
237 "Thyra's communicator is neither an MpiComm nor a SerialComm. "
238 "Sorry, I have no idea what to do with it in that case.");
239 // It's a SerialComm. Make a SerialComm of our own.
240 // SerialComm instances are all the same, so there's no need
241 // to keep the original one.
242 return rcp_implicit_cast<const Teuchos::Comm<int> > (rcp (new Teuchos::SerialComm<int>));
243 }
244 else { // Yippie, we have an MpiComm.
245 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = thyraMpiComm->getRawMpiComm ();
246 // NOTE (mfh 18 Jun 2013) Since the error handler is attached
247 // to the MPI_Comm, not to the Teuchos widget, we don't have
248 // to set the error handler again on the new MpiComm object.
249 return rcp_implicit_cast<const Teuchos::Comm<int> > (rcp (new Teuchos::MpiComm<int> (rawMpiComm)));
250 }
251#else // NOT HAVE_MPI
252 // Either it's a SerialComm or I don't know what to do with it.
253 RCP<const Teuchos::SerialComm<Thyra::Ordinal> > thyraSerialComm =
254 rcp_dynamic_cast<const Teuchos::SerialComm<Thyra::Ordinal> > (thyraComm);
255 TEUCHOS_TEST_FOR_EXCEPTION(
256 thyraSerialComm.is_null (), std::runtime_error,
257 "Thyra's communicator is not a SerialComm, and MPI is not enabled, so "
258 "it can't be an MpiComm either. That means it must be some other "
259 "subclass of Comm, about which I don't know. "
260 "Sorry, I have no idea what to do with it in that case.");
261 // It's a SerialComm. Make a SerialComm of our own.
262 // SerialComm instances are all the same, so there's no need
263 // to keep the original one.
264 return rcp_implicit_cast<const Teuchos::Comm<int> > (rcp (new Teuchos::SerialComm<int>));
265#endif // HAVE_MPI
266 }
267
280 void
281 prepareDistTsqr (const MV& /* X */) {}
282
302 void
303 prepareTsqr (const MV& /* X */) {}
304 };
305
306} // namespace Tpetra
307
308#endif // __Thyra_TsqrAdaptor_hpp
309
Stub adaptor from Thyra::MultiVectorBase to TSQR.
int revealRank(MV &, dense_matrix_type &, const magnitude_type &)
Rank-revealing decomposition.
TsqrAdaptor()
Constructor (that uses default parameters).
TsqrAdaptor(const Teuchos::RCP< Teuchos::ParameterList > &)
Constructor (that accepts a parameter list).
void factorExplicit(MV &, MV &, dense_matrix_type &, const bool=false)
Compute QR factorization [Q,R] = qr(A,0).

Generated for Stratimikos by doxygen 1.9.6