Teko Version of the Day
Loading...
Searching...
No Matches
Teko_BlockedEpetraOperator.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_BlockedEpetraOperator.hpp"
48#include "Teko_BlockedMappingStrategy.hpp"
49#include "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 "EpetraExt_MultiVectorOut.h"
62#include "EpetraExt_RowMatrixOut.h"
63
64#include "Teko_Utilities.hpp"
65
66namespace Teko {
67namespace Epetra {
68
69using Teuchos::RCP;
70using Teuchos::rcp;
71using Teuchos::rcp_dynamic_cast;
72
73BlockedEpetraOperator::BlockedEpetraOperator(const std::vector<std::vector<int> > & vars,
74 const Teuchos::RCP<const Epetra_Operator> & content,
75 const std::string & label)
76 : Teko::Epetra::EpetraOperatorWrapper(), label_(label)
77{
78 SetContent(vars,content);
79}
80
81void BlockedEpetraOperator::SetContent(const std::vector<std::vector<int> > & vars,
82 const Teuchos::RCP<const Epetra_Operator> & content)
83{
84 fullContent_ = content;
85 blockedMapping_ = rcp(new BlockedMappingStrategy(vars,Teuchos::rcpFromRef(fullContent_->OperatorDomainMap()),
86 fullContent_->Comm()));
87 SetMapStrategy(blockedMapping_);
88
89 // build thyra operator
90 BuildBlockedOperator();
91}
92
93void BlockedEpetraOperator::BuildBlockedOperator()
94{
95 TEUCHOS_ASSERT(blockedMapping_!=Teuchos::null);
96
97 // get a CRS matrix
98 const RCP<const Epetra_CrsMatrix> crsContent
99 = rcp_dynamic_cast<const Epetra_CrsMatrix>(fullContent_);
100
101 // ask the strategy to build the Thyra operator for you
102 if(blockedOperator_==Teuchos::null) {
103 blockedOperator_ = blockedMapping_->buildBlockedThyraOp(crsContent,label_);
104 }
105 else {
106 const RCP<Thyra::BlockedLinearOpBase<double> > blkOp
107 = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<double> >(blockedOperator_,true);
108 blockedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
109 }
110
111 // set whatever is returned
112 SetOperator(blockedOperator_,false);
113
114 // reorder if neccessary
115 if(reorderManager_!=Teuchos::null)
116 Reorder(*reorderManager_);
117}
118
119const Teuchos::RCP<const Epetra_Operator> BlockedEpetraOperator::GetBlock(int i,int j) const
120{
121 const RCP<const Thyra::BlockedLinearOpBase<double> > blkOp
122 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(getThyraOp());
123
124 return Thyra::get_Epetra_Operator(*blkOp->getBlock(i,j));
125}
126
131{
132 reorderManager_ = rcp(new BlockReorderManager(brm));
133
134 // build reordered objects
135 RCP<const MappingStrategy> reorderMapping = rcp(new ReorderedMappingStrategy(*reorderManager_,blockedMapping_));
136 RCP<const Thyra::BlockedLinearOpBase<double> > blockOp
137 = rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<double> >(blockedOperator_);
138
139 RCP<const Thyra::LinearOpBase<double> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
140
141 // set them as working values
142 SetMapStrategy(reorderMapping);
143 SetOperator(A,false);
144}
145
148{
149 SetMapStrategy(blockedMapping_);
150 SetOperator(blockedOperator_,false);
151 reorderManager_ = Teuchos::null;
152}
153
156void BlockedEpetraOperator::WriteBlocks(const std::string & prefix) const
157{
158 RCP<Thyra::PhysicallyBlockedLinearOpBase<double> > blockOp
159 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<double> >(blockedOperator_);
160
161 // get size of blocked block operator
162 int rows = Teko::blockRowCount(blockOp);
163
164 for(int i=0;i<rows;i++) {
165 for(int j=0;j<rows;j++) {
166 // build the file name
167 std::stringstream ss;
168 ss << prefix << "_" << i << j << ".mm";
169
170 // get the row matrix object
171 RCP<const Epetra_RowMatrix> mat
172 = Teuchos::rcp_dynamic_cast<const Epetra_RowMatrix>(Thyra::get_Epetra_Operator(*blockOp->getBlock(i,j)));
173
174 // write to file
175 EpetraExt::RowMatrixToMatrixMarketFile(ss.str().c_str(),*mat);
176 }
177 }
178}
179
181{
182 Epetra_Vector xf(OperatorRangeMap());
183 Epetra_Vector xs(OperatorRangeMap());
184 Epetra_Vector y(OperatorDomainMap());
185
186 // test operator many times
187 bool result = true;
188 double diffNorm=0.0,trueNorm=0.0;
189 for(int i=0;i<count;i++) {
190 xf.PutScalar(0.0);
191 xs.PutScalar(0.0);
192 y.Random();
193
194 // apply operator
195 Apply(y,xs); // xs = A*y
196 fullContent_->Apply(y,xf); // xf = A*y
197
198 // compute norms
199 xs.Update(-1.0,xf,1.0);
200 xs.Norm2(&diffNorm);
201 xf.Norm2(&trueNorm);
202
203 // check result
204 result &= (diffNorm/trueNorm < tol);
205 }
206 return result;
207}
208
209} // end namespace Epetra
210} // end namespace Teko
int blockRowCount(const BlockedLinearOp &blo)
Get the row count in a block linear operator.
Class that describes how a flat blocked operator should be reordered.
virtual void SetContent(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< const Epetra_Operator > &content)
void Reorder(const BlockReorderManager &brm)
BlockedEpetraOperator(const std::vector< std::vector< int > > &vars, const Teuchos::RCP< const Epetra_Operator > &content, const std::string &label="<ANYM>")
void RemoveReording()
Remove any reordering on this object.
virtual void WriteBlocks(const std::string &prefix) const
bool testAgainstFullOperator(int count, double tol) const
Helps perform sanity checks.
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...
const RCP< const Thyra::LinearOpBase< double > > getThyraOp() const
Return the thyra operator associated with this wrapper.