Teko Version of the Day
Loading...
Searching...
No Matches
Teko_EpetraHelpers.cpp
1/*
2// @HEADER
3//
4// ***********************************************************************
5//
6// Teko: A package for block and physics based preconditioning
7// Copyright 2010 Sandia Corporation
8//
9// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
40//
41// ***********************************************************************
42//
43// @HEADER
44
45*/
46
47#include "Teko_EpetraHelpers.hpp"
48
49// Thyra Includes
50#include "Thyra_EpetraLinearOp.hpp"
51#include "Thyra_BlockedLinearOpBase.hpp"
52#include "Thyra_DefaultMultipliedLinearOp.hpp"
53#include "Thyra_DefaultDiagonalLinearOp.hpp"
54#include "Thyra_DefaultZeroLinearOp.hpp"
55#include "Thyra_DefaultBlockedLinearOp.hpp"
56#include "Thyra_EpetraThyraWrappers.hpp"
57#include "Thyra_SpmdVectorBase.hpp"
58#include "Thyra_SpmdVectorSpaceBase.hpp"
59#include "Thyra_ScalarProdVectorSpaceBase.hpp"
60
61// Epetra includes
62#include "Epetra_Vector.h"
63
64// EpetraExt includes
65#include "EpetraExt_ProductOperator.h"
66#include "EpetraExt_MatrixMatrix.h"
67
68// Teko includes
69#include "Teko_EpetraOperatorWrapper.hpp"
70#include "Teko_Utilities.hpp"
71
72using Teuchos::RCP;
73using Teuchos::rcp;
74using Teuchos::rcpFromRef;
75using Teuchos::rcp_dynamic_cast;
76using Teuchos::null;
77
78namespace Teko {
79namespace Epetra {
80
91const Teuchos::RCP<const Thyra::LinearOpBase<double> > thyraDiagOp(const RCP<const Epetra_Vector> & ev,const Epetra_Map & map,
92 const std::string & lbl)
93{
94 const RCP<const Thyra::VectorBase<double> > thyraVec // need a Thyra::VectorBase object
95 = Thyra::create_Vector(ev,Thyra::create_VectorSpace(rcpFromRef(map)));
96 Teuchos::RCP<Thyra::LinearOpBase<double> > op
97 = Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<double>(thyraVec));
98 op->setObjectLabel(lbl);
99 return op;
100}
101
112const Teuchos::RCP<Thyra::LinearOpBase<double> > thyraDiagOp(const RCP<Epetra_Vector> & ev,const Epetra_Map & map,
113 const std::string & lbl)
114{
115 const RCP<Thyra::VectorBase<double> > thyraVec // need a Thyra::VectorBase object
116 = Thyra::create_Vector(ev,Thyra::create_VectorSpace(rcpFromRef(map)));
117 Teuchos::RCP<Thyra::LinearOpBase<double> > op
118 = Teuchos::rcp(new Thyra::DefaultDiagonalLinearOp<double>(thyraVec));
119 op->setObjectLabel(lbl);
120 return op;
121}
122
132void fillDefaultSpmdMultiVector(Teuchos::RCP<Thyra::DefaultSpmdMultiVector<double> > & spmdMV,
133 Teuchos::RCP<Epetra_MultiVector> & epetraMV)
134{
135 // first get desired range and domain
136 const RCP<const Thyra::SpmdVectorSpaceBase<double> > range = spmdMV->spmdSpace();
137 const RCP<const Thyra::ScalarProdVectorSpaceBase<double> > domain
138 = rcp_dynamic_cast<const Thyra::ScalarProdVectorSpaceBase<double> >(spmdMV->domain());
139
140 TEUCHOS_ASSERT(domain->dim()==epetraMV->NumVectors());
141
142 // New local view of raw data
143 double *localValues; int leadingDim;
144 if(epetraMV->ConstantStride() )
145 epetraMV->ExtractView( &localValues, &leadingDim );
146 else
147 TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
148
149 // Build the MultiVector
150 spmdMV->initialize(range, domain,
151 Teuchos::arcp(localValues,0,leadingDim*epetraMV->NumVectors(),false),
152 leadingDim);
153
154 // make sure the Epetra_MultiVector doesn't disappear prematurely
155 Teuchos::set_extra_data<RCP<Epetra_MultiVector> >(epetraMV,"Epetra_MultiVector",Teuchos::outArg(spmdMV));
156}
157
167void identityRowIndices(const Epetra_Map & rowMap, const Epetra_CrsMatrix & mat,std::vector<int> & outIndices)
168{
169 int maxSz = mat.GlobalMaxNumEntries();
170 std::vector<double> values(maxSz);
171 std::vector<int> indices(maxSz);
172
173 // loop over elements owned by this processor
174 for(int i=0;i<rowMap.NumMyElements();i++) {
175 bool rowIsIdentity = true;
176 int sz = 0;
177 int rowGID = rowMap.GID(i);
178 mat.ExtractGlobalRowCopy(rowGID,maxSz,sz,&values[0],&indices[0]);
179
180 // loop over the columns of this row
181 for(int j=0;j<sz;j++) {
182 int colGID = indices[j];
183
184 // look at row entries
185 if(colGID==rowGID) rowIsIdentity &= values[j]==1.0;
186 else rowIsIdentity &= values[j]==0.0;
187
188 // not a dirchlet row...quit
189 if(not rowIsIdentity)
190 break;
191 }
192
193 // save a row that is dirchlet
194 if(rowIsIdentity)
195 outIndices.push_back(rowGID);
196 }
197}
198
209void zeroMultiVectorRowIndices(Epetra_MultiVector & mv,const std::vector<int> & zeroIndices)
210{
211 int colCnt = mv.NumVectors();
212 std::vector<int>::const_iterator itr;
213
214 // loop over the indices to zero
215 for(itr=zeroIndices.begin();itr!=zeroIndices.end();++itr) {
216
217 // loop over columns
218 for(int j=0;j<colCnt;j++)
219 TEUCHOS_TEST_FOR_EXCEPT(mv.ReplaceGlobalValue(*itr,j,0.0));
220 }
221}
222
232ZeroedOperator::ZeroedOperator(const std::vector<int> & zeroIndices,
233 const Teuchos::RCP<const Epetra_Operator> & op)
234 : zeroIndices_(zeroIndices), epetraOp_(op)
235{
236 label_ = "zeroed( "+std::string(epetraOp_->Label())+" )";
237}
238
240int ZeroedOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
241{
242/*
243 Epetra_MultiVector temp(X);
244 zeroMultiVectorRowIndices(temp,zeroIndices_);
245 int result = epetraOp_->Apply(temp,Y);
246*/
247
248 int result = epetraOp_->Apply(X,Y);
249
250 // zero a few of the rows
251 zeroMultiVectorRowIndices(Y,zeroIndices_);
252
253 return result;
254}
255
256} // end namespace Epetra
257} // end namespace Teko
ZeroedOperator(const std::vector< int > &zeroIndices, const Teuchos::RCP< const Epetra_Operator > &op)
Constructor for a ZeroedOperator.
int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Perform a matrix-vector product with certain rows zeroed out.