Teko Version of the Day
Loading...
Searching...
No Matches
Teko_StridedEpetraOperator.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 "Epetra/Teko_StridedEpetraOperator.hpp"
48#include "Epetra/Teko_StridedMappingStrategy.hpp"
49#include "Epetra/Teko_ReorderedMappingStrategy.hpp"
50
51#include "Teuchos_VerboseObject.hpp"
52
53#include "Thyra_LinearOpBase.hpp"
54#include "Thyra_EpetraLinearOp.hpp"
55#include "Thyra_EpetraThyraWrappers.hpp"
56#include "Thyra_DefaultProductMultiVector.hpp"
57#include "Thyra_DefaultProductVectorSpace.hpp"
58#include "Thyra_DefaultBlockedLinearOp.hpp"
59#include "Thyra_get_Epetra_Operator.hpp"
60
61#include "Epetra_Vector.h"
62#include "EpetraExt_MultiVectorOut.h"
63#include "EpetraExt_RowMatrixOut.h"
64
65#include "Teko_Utilities.hpp"
66
67namespace Teko {
68namespace Epetra {
69
70using Teuchos::RCP;
71using Teuchos::rcp;
72using Teuchos::rcp_dynamic_cast;
73
74StridedEpetraOperator::StridedEpetraOperator(int numVars,const Teuchos::RCP<const Epetra_Operator> & content,
75 const std::string & label)
76 : Teko::Epetra::EpetraOperatorWrapper(), label_(label)
77{
78 std::vector<int> vars;
79
80 // build vector describing the sub maps
81 for(int i=0;i<numVars;i++) vars.push_back(1);
82
83 SetContent(vars,content);
84}
85
86StridedEpetraOperator::StridedEpetraOperator(const std::vector<int> & vars,const Teuchos::RCP<const Epetra_Operator> & content,
87 const std::string & label)
88 : Teko::Epetra::EpetraOperatorWrapper(), label_(label)
89{
90 SetContent(vars,content);
91}
92
93void StridedEpetraOperator::SetContent(const std::vector<int> & vars,const Teuchos::RCP<const Epetra_Operator> & content)
94{
95 fullContent_ = content;
96 stridedMapping_ = rcp(new StridedMappingStrategy(vars,Teuchos::rcpFromRef(fullContent_->OperatorDomainMap()),
97 fullContent_->Comm()));
98 SetMapStrategy(stridedMapping_);
99
100 // build thyra operator
101 BuildBlockedOperator();
102}
103
104void StridedEpetraOperator::BuildBlockedOperator()
105{
106 TEUCHOS_ASSERT(stridedMapping_!=Teuchos::null);
107
108 // get a CRS matrix
109 const RCP<const Epetra_CrsMatrix> crsContent = rcp_dynamic_cast<const Epetra_CrsMatrix>(fullContent_);
110
111 // ask the strategy to build the Thyra operator for you
112 if(stridedOperator_==Teuchos::null) {
113 stridedOperator_ = stridedMapping_->buildBlockedThyraOp(crsContent,label_);
114 }
115 else {
116 const RCP<Thyra::BlockedLinearOpBase<double> > blkOp = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(stridedOperator_,true);
117 stridedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
118 }
119
120 // set whatever is returned
121 SetOperator(stridedOperator_,false);
122
123 // reorder if neccessary
124 if(reorderManager_!=Teuchos::null)
125 Reorder(*reorderManager_);
126}
127
128const Teuchos::RCP<const Epetra_Operator> StridedEpetraOperator::GetBlock(int i,int j) const
129{
130 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
131 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
132
133 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i,j));
134}
135
139void StridedEpetraOperator::Reorder(const BlockReorderManager & brm)
140{
141 reorderManager_ = rcp(new BlockReorderManager(brm));
142
143 // build reordered objects
144 RCP<const MappingStrategy> reorderMapping = rcp(new ReorderedMappingStrategy(*reorderManager_,stridedMapping_));
145 RCP<const Thyra::BlockedLinearOpBase<double> > blockOp
146 = rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(stridedOperator_);
147
148 RCP<const Thyra::LinearOpBase<double> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
149
150 // set them as working values
151 SetMapStrategy(reorderMapping);
152 SetOperator(A,false);
153}
154
156void StridedEpetraOperator::RemoveReording()
157{
158 SetMapStrategy(stridedMapping_);
159 SetOperator(stridedOperator_,false);
160 reorderManager_ = Teuchos::null;
161}
162
165void StridedEpetraOperator::WriteBlocks(const std::string & prefix) const
166{
167 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockOp
168 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> >(stridedOperator_);
169
170 // get size of strided block operator
171 int rows = Teko::blockRowCount(blockOp);
172
173 for(int i=0;i<rows;i++) {
174 for(int j=0;j<rows;j++) {
175 // build the file name
176 std::stringstream ss;
177 ss << prefix << "_" << i << j << ".mm";
178
179 // get the row matrix object (Note: can't use "GetBlock" method b/c matrix might be reordered)
180 RCP<const Epetra_RowMatrix> mat
181 = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(Thyra::get_Epetra_Operator(*blockOp->getBlock(i,j)));
182
183 // write to file
184 EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*mat);
185 }
186 }
187}
188
196std::string StridedEpetraOperator::PrintNorm(const eNormType & nrmType,const char newline)
197{
198 BlockedLinearOp blockOp = toBlockedLinearOp(stridedOperator_);
199
200 // get size of strided block operator
201 int rowCount = Teko::blockRowCount(blockOp);
202 int colCount = Teko::blockRowCount(blockOp);
203
204 std::stringstream ss;
205 ss.precision(4);
206 ss << std::scientific;
207 for(int row=0;row<rowCount;row++) {
208 for(int col=0;col<colCount;col++) {
209 // get the row matrix object (Note: can't use "GetBlock" method b/c matrix might be reordered)
210 RCP<const Epetra_CrsMatrix> mat = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(
211 Thyra::get_Epetra_Operator(*blockOp->getBlock(row,col)));
212
213 // compute the norm
214 double norm = 0.0;
215 switch(nrmType) {
216 case Inf:
217 norm = mat->NormInf();
218 break;
219 case One:
220 norm = mat->NormOne();
221 break;
222 case Frobenius:
223 norm = mat->NormFrobenius();
224 break;
225 default:
226 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
227 "StridedEpetraOperator::eNormType incorrectly specified");
228 }
229
230 ss << norm << " ";
231 }
232 ss << newline;
233 }
234
235 return ss.str();
236}
237
238bool StridedEpetraOperator::testAgainstFullOperator(int count,double tol) const
239{
240 Epetra_Vector xf(OperatorRangeMap());
241 Epetra_Vector xs(OperatorRangeMap());
242 Epetra_Vector y(OperatorDomainMap());
243
244 // test operator many times
245 bool result = true;
246 double diffNorm=0.0,trueNorm=0.0;
247 for(int i=0;i<count;i++) {
248 xf.PutScalar(0.0);
249 xs.PutScalar(0.0);
250 y.Random();
251
252 // apply operator
253 Apply(y,xs); // xs = A*y
254 fullContent_->Apply(y,xf); // xf = A*y
255
256 // compute norms
257 xs.Update(-1.0,xf,1.0);
258 xs.Norm2(&diffNorm);
259 xf.Norm2(&trueNorm);
260
261 // check result
262 result &= (diffNorm/trueNorm < tol);
263 }
264
265 return result;
266}
267
268} // end namespace Epetra
269} // end namespace Teko
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
BlockedLinearOp toBlockedLinearOp(LinearOp &clo)
Converts a LinearOp to a BlockedLinearOp.