Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosMatOrthoManager.hpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Belos: Block Linear Solvers Package
5// Copyright 2004 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// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
47#ifndef BELOS_MATORTHOMANAGER_HPP
48#define BELOS_MATORTHOMANAGER_HPP
49
69#include "BelosConfigDefs.hpp"
70#include "BelosTypes.hpp"
71#include "BelosOrthoManager.hpp"
74
75namespace Belos {
76
77 template <class ScalarType, class MV, class OP>
78 class MatOrthoManager : public OrthoManager<ScalarType,MV> {
79 protected:
80 Teuchos::RCP<const OP> _Op;
81 bool _hasOp;
82
83 public:
85
86
87 MatOrthoManager(Teuchos::RCP<const OP> Op = Teuchos::null) : _Op(Op), _hasOp(Op!=Teuchos::null) {};
88
90 virtual ~MatOrthoManager() {};
92
94
95
97 void setOp( Teuchos::RCP<const OP> Op ) {
98 _Op = Op;
99 _hasOp = (_Op != Teuchos::null);
100 };
101
103 Teuchos::RCP<const OP> getOp() const { return _Op; }
104
106
107
109
110
115 void innerProd( const MV& X, const MV& Y,
116 Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const {
117 typedef Teuchos::ScalarTraits<ScalarType> SCT;
120
121 Teuchos::RCP<const MV> P,Q;
122 Teuchos::RCP<MV> R;
123
124 if (_hasOp) {
125 // attempt to minimize the amount of work in applying
126 if ( MVT::GetNumberVecs(X) < MVT::GetNumberVecs(Y) ) {
127 R = MVT::Clone(X,MVT::GetNumberVecs(X));
128 OPT::Apply(*_Op,X,*R);
129 P = R;
130 Q = Teuchos::rcp( &Y, false );
131 }
132 else {
133 P = Teuchos::rcp( &X, false );
134 R = MVT::Clone(Y,MVT::GetNumberVecs(Y));
135 OPT::Apply(*_Op,Y,*R);
136 Q = R;
137 }
138 }
139 else {
140 P = Teuchos::rcp( &X, false );
141 Q = Teuchos::rcp( &Y, false );
142 }
143
144 MVT::MvTransMv(SCT::one(),*P,*Q,Z);
145 }
146
153 void innerProd( const MV& X, const MV& Y, Teuchos::RCP<const MV> MY,
154 Teuchos::SerialDenseMatrix<int,ScalarType>& Z ) const {
155 typedef Teuchos::ScalarTraits<ScalarType> SCT;
157
158 Teuchos::RCP<MV> P,Q;
159
160 if ( MY == Teuchos::null ) {
161 innerProd(X,Y,Z);
162 }
163 else if ( _hasOp ) {
164 // the user has done the matrix vector for us
165 MVT::MvTransMv(SCT::one(),X,*MY,Z);
166 }
167 else {
168 // there is no matrix vector
169 MVT::MvTransMv(SCT::one(),X,Y,Z);
170 }
171 }
172
175 void norm( const MV& X, std::vector< typename Teuchos::ScalarTraits<ScalarType>::magnitudeType >& normvec ) const {
176 norm(X,Teuchos::null,normvec);
177 }
178
196 void
197 norm (const MV& X,
198 Teuchos::RCP<const MV> MX,
199 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>& normvec) const
200 {
201 typedef Teuchos::ScalarTraits<ScalarType> SCT;
202 typedef Teuchos::ScalarTraits<typename SCT::magnitudeType> MT;
205
206 int nvecs = MVT::GetNumberVecs(X);
207
208 // Make sure that normvec has enough entries to hold the norms
209 // of all the columns of X. std::vector<T>::size_type is
210 // unsigned, so do the appropriate cast to avoid signed/unsigned
211 // comparisons that trigger compiler warnings.
212 if (normvec.size() < static_cast<size_t>(nvecs))
213 normvec.resize (nvecs);
214
215 if (!_hasOp) {
216 // X == MX, since the operator M is the identity.
217 MX = Teuchos::rcp(&X, false);
218 MVT::MvNorm(X, normvec);
219 }
220 else {
221 // The caller didn't give us a previously computed MX, so
222 // apply the operator. We assign to MX only after applying
223 // the operator, so that if the application fails, MX won't be
224 // modified.
225 if(MX == Teuchos::null) {
226 Teuchos::RCP<MV> tempVec = MVT::Clone(X,nvecs);
227 OPT::Apply(*_Op,X,*tempVec);
228 MX = tempVec;
229 }
230 else {
231 // The caller gave us a previously computed MX. Make sure
232 // that it has at least as many columns as X.
233 const int numColsMX = MVT::GetNumberVecs(*MX);
234 TEUCHOS_TEST_FOR_EXCEPTION(numColsMX < nvecs, std::invalid_argument,
235 "MatOrthoManager::norm(X, MX, normvec): "
236 "MX has fewer columns than X: "
237 "MX has " << numColsMX << " columns, "
238 "and X has " << nvecs << " columns.");
239 }
240
241 std::vector<ScalarType> dotvec(nvecs);
242 MVT::MvDot(X,*MX,dotvec);
243 for (int i=0; i<nvecs; i++) {
244 normvec[i] = MT::squareroot( SCT::magnitude(dotvec[i]) );
245 }
246 }
247 }
248
249
271 virtual void project ( MV &X, Teuchos::RCP<MV> MX,
272 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
273 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
274
275
276
279 virtual void project ( MV &X,
280 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
281 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
282 project(X,Teuchos::null,C,Q);
283 }
284
306 virtual int normalize ( MV &X, Teuchos::RCP<MV> MX,
307 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const = 0;
308
309
312 virtual int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
313 return normalize(X,Teuchos::null,B);
314 }
315
316
317 protected:
318 virtual int
320 Teuchos::RCP<MV> MX,
321 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
322 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
323 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const = 0;
324
325 virtual int
327 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
328 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
329 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
330 {
331 return this->projectAndNormalizeWithMxImpl (X, Teuchos::null, C, B, Q);
332 }
333
334 public:
335
370 int
372 Teuchos::RCP<MV> MX,
373 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
374 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
375 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
376 {
377 return this->projectAndNormalizeWithMxImpl (X, MX, C, B, Q);
378 }
379
381
383
386 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
387 orthonormError(const MV &X) const {
388 return orthonormError(X,Teuchos::null);
389 }
390
394 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
395 orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const = 0;
396
399 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
400 orthogError(const MV &X1, const MV &X2) const {
401 return orthogError(X1,Teuchos::null,X2);
402 }
403
408 virtual typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
409 orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const = 0;
410
412
413 };
414
415} // end of Belos namespace
416
417
418#endif
419
420// end of file BelosMatOrthoManager.hpp
Belos header file which uses auto-configuration information to include necessary C++ headers.
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Collection of types and exceptions used within the Belos solvers.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
void innerProd(const MV &X, const MV &Y, Teuchos::RCP< const MV > MY, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator....
virtual int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method calls normalize(X,Teuchos::null,B); see documentation for that function.
virtual ~MatOrthoManager()
Destructor.
MatOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null)
Default constructor.
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator.
int projectAndNormalize(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for .
virtual int projectAndNormalizeWithMxImpl(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const =0
virtual Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogError(const MV &X1, Teuchos::RCP< const MV > MX1, const MV &X2) const =0
This method computes the error in orthogonality of two multivectors. The method has the option of exp...
virtual void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const =0
Given a list of (mutually and internally) orthonormal bases Q, this method takes a multivector X and ...
virtual Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
void norm(const MV &X, Teuchos::RCP< const MV > MX, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec) const
Compute norm of each column of X.
void norm(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec) const
Provides the norm induced by innerProd().
virtual Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X, Teuchos::RCP< const MV > MX) const =0
This method computes the error in orthonormality of a multivector. The method has the option of explo...
virtual int normalize(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const =0
This method takes a multivector X and attempts to compute an orthonormal basis for ,...
virtual int projectAndNormalizeImpl(MV &X, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
void setOp(Teuchos::RCP< const OP > Op)
Set operator.
virtual void project(MV &X, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
This method calls project(X,Teuchos::null,C,Q); see documentation for that function.
Teuchos::RCP< const OP > _Op
virtual Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors. This method.
Teuchos::RCP< const OP > getOp() const
Get operator.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...