Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_DefaultSerialDenseLinearOpWithSolve_def.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) 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 (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42
43#ifndef THYRA_DEFAULT_SERIAL_DENSE_LINEAR_OP_WITH_SOLVE_HPP
44#define THYRA_DEFAULT_SERIAL_DENSE_LINEAR_OP_WITH_SOLVE_HPP
45
46
47#include "Thyra_DefaultSerialDenseLinearOpWithSolve_decl.hpp"
48#include "Thyra_LinearOpWithSolveBase.hpp"
49#include "Thyra_DetachedMultiVectorView.hpp"
50#include "Thyra_MultiVectorStdOps.hpp"
51#include "Thyra_AssertOp.hpp"
52#include "Teuchos_Assert.hpp"
53
54
55namespace Thyra {
56
57
58// Constructors/initializers/accessors
59
60
61template<class Scalar>
63{}
64
65
66template<class Scalar>
68 const RCP<const MultiVectorBase<Scalar> > &M )
69{
70 using Teuchos::outArg;
71#ifdef TEUCHOS_DEBUG
72 TEUCHOS_ASSERT(!is_null(M));
73 TEUCHOS_ASSERT(isFullyInitialized(*M));
74 TEUCHOS_ASSERT(M->range()->hasInCoreView());
75 TEUCHOS_ASSERT(M->domain()->hasInCoreView());
76 THYRA_ASSERT_VEC_SPACES("", *M->range(), *M->domain());
77#endif
78 factorize(*M, outArg(LU_), outArg(ipiv_));
79 M_ = M;
80}
81
82template<class Scalar>
84{
85 return M_;
86}
87
88// Overridden from LinearOpBase
89
90
91template<class Scalar>
94{
95 if (!is_null(M_))
96 return M_->range();
97 return Teuchos::null;
98}
99
100
101template<class Scalar>
104{
105 if (!is_null(M_))
106 return M_->domain();
107 return Teuchos::null;
108}
109
110
111// protected
112
113
114// Overridden from LinearOpBase
115
116
117template<class Scalar>
119 EOpTransp M_trans) const
120{
121 return Thyra::opSupported(*M_, M_trans);
122}
123
124
125template<class Scalar>
127 const EOpTransp M_trans,
129 const Ptr<MultiVectorBase<Scalar> > &Y,
130 const Scalar alpha,
131 const Scalar beta
132 ) const
133{
134 Thyra::apply( *M_, M_trans, X, Y, alpha, beta );
135}
136
137
138// Overridden from LinearOpWithSolveBase
139
140
141template<class Scalar>
143 EOpTransp M_trans) const
144{
146 return ( ST::isComplex ? ( M_trans!=CONJ ) : true );
147}
148
149
150template<class Scalar>
152 EOpTransp M_trans, const SolveMeasureType& /* solveMeasureType */) const
153{
154 // We support all solve measures since we are a direct solver
155 return this->solveSupportsImpl(M_trans);
156}
157
158
159template<class Scalar>
162 const EOpTransp M_trans,
164 const Ptr<MultiVectorBase<Scalar> > &X,
165 const Ptr<const SolveCriteria<Scalar> > /* solveCriteria */
166 ) const
167{
168#ifdef TEUCHOS_DEBUG
170 "DefaultSerialDenseLinearOpWithSolve<Scalar>::solve(...)",
171 *this, M_trans, *X, &B );
172#endif
173 backsolve( LU_, ipiv_, M_trans, B, X );
174 SolveStatus<Scalar> solveStatus;
175 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;\
176 return solveStatus;
177}
178
179
180// private
181
182
183template<class Scalar>
187 const Ptr<Array<int> > &ipiv
188 )
189{
190 using Teuchos::outArg;
192 const int dim = dM.subDim();
193 ipiv->resize(dim);
195 RTOpPack::assign_entries<Scalar>( outArg(LU_tmp), dM.smv() );
196 int rank = -1;
197 RTOpPack::getrf<Scalar>( LU_tmp, (*ipiv)(), outArg(rank) );
198 TEUCHOS_ASSERT_EQUALITY( dim, rank );
199 *LU = LU_tmp; // Weak copy
200}
201
202
203template<class Scalar>
204void DefaultSerialDenseLinearOpWithSolve<Scalar>::backsolve(
206 const ArrayView<const int> ipiv,
207 const EOpTransp transp,
208 const MultiVectorBase<Scalar> &B,
209 const Ptr<MultiVectorBase<Scalar> > &X
210 )
211{
212 using Teuchos::outArg;
213 assign( X, B );
214 DetachedMultiVectorView<Scalar> dX(*X);
215 RTOpPack::getrs<Scalar>( LU, ipiv, convertToRTOpPackETransp(transp),
216 outArg(dX.smv()) );
217}
218
219
220} // end namespace Thyra
221
222
223#endif // THYRA_DEFAULT_SERIAL_DENSE_LINEAR_OP_WITH_SOLVE_HPP
Create an explicit non-mutable (const) view of a MultiVectorBase object.
Simple concreate subclass of LinearOpWithSolveBase for serial dense matrices implemented using LAPACK...
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
void initialize(const RCP< const MultiVectorBase< Scalar > > &M)
SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
Interface for a collection of column vectors called a multi-vector.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
@ SOLVE_STATUS_CONVERGED
The requested solution criteria has likely been achieved.
#define THYRA_ASSERT_VEC_SPACES(FUNC_NAME, VS1, VS2)
This is a very useful macro that should be used to validate that two vector spaces are compatible.
#define THYRA_ASSERT_LINEAR_OP_MULTIVEC_APPLY_SPACES(FUNC_NAME, M, M_T, X, Y)
This is a very useful macro that should be used to validate that the spaces for the multi-vector vers...
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
@ CONJ
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
Simple struct that defines the requested solution criteria for a solve.
Simple struct for the return status from a solve.
ESolveStatus solveStatus
The return status of the solve.