Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
55 namespace Thyra {
56 
57 
58 // Constructors/initializers/accessors
59 
60 
61 template<class Scalar>
63 {}
64 
65 
66 template<class Scalar>
68  const RCP<const MultiVectorBase<Scalar> > &M )
69 {
70  using Teuchos::outArg;
71 #ifdef TEUCHOS_DEBUG
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 
82 template<class Scalar>
84 {
85  return M_;
86 }
87 
88 // Overridden from LinearOpBase
89 
90 
91 template<class Scalar>
94 {
95  if (!is_null(M_))
96  return M_->range();
97  return Teuchos::null;
98 }
99 
100 
101 template<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 
117 template<class Scalar>
119  EOpTransp M_trans) const
120 {
121  return Thyra::opSupported(*M_, M_trans);
122 }
123 
124 
125 template<class Scalar>
127  const EOpTransp M_trans,
128  const MultiVectorBase<Scalar> &X,
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 
141 template<class Scalar>
143  EOpTransp M_trans) const
144 {
146  return ( ST::isComplex ? ( M_trans!=CONJ ) : true );
147 }
148 
149 
150 template<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 
159 template<class Scalar>
162  const EOpTransp M_trans,
163  const MultiVectorBase<Scalar> &B,
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 
183 template<class Scalar>
185  const MultiVectorBase<Scalar> &M,
187  const Ptr<Array<int> > &ipiv
188  )
189 {
190  using Teuchos::outArg;
192  const int dim = dM.subDim();
193  ipiv->resize(dim);
194  RTOpPack::SubMultiVectorView<Scalar> LU_tmp(dim, 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 
203 template<class Scalar>
204 void 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
#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...
Create an explicit non-mutable (const) view of a MultiVectorBase object.
SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
bool is_null(const boost::shared_ptr< T > &p)
EOpTransp
Enumeration for determining how a linear operator is applied. `*.
#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...
Simple concreate subclass of LinearOpWithSolveBase for serial dense matrices implemented using LAPACK...
Use the non-transposed operator with complex-conjugate elements (same as NOTRANS for real scalar type...
bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
void initialize(const RCP< const MultiVectorBase< Scalar > > &M)
Interface for a collection of column vectors called a multi-vector.
Simple struct for the return status from a solve.
ESolveStatus solveStatus
The return status of the solve.
void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
The requested solution criteria has likely been achieved.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
Simple struct that defines the requested solution criteria for a solve.