Teko Version of the Day
Loading...
Searching...
No Matches
Teko_BlockedTpetraOperator.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_BlockedTpetraOperator.hpp"
48#include "Teko_TpetraBlockedMappingStrategy.hpp"
49#include "Teko_TpetraReorderedMappingStrategy.hpp"
50
51#include "Teuchos_VerboseObject.hpp"
52
53#include "Thyra_LinearOpBase.hpp"
54#include "Thyra_TpetraLinearOp.hpp"
55#include "Thyra_TpetraThyraWrappers.hpp"
56#include "Thyra_DefaultProductMultiVector.hpp"
57#include "Thyra_DefaultProductVectorSpace.hpp"
58#include "Thyra_DefaultBlockedLinearOp.hpp"
59
60#include "MatrixMarket_Tpetra.hpp"
61
62#include "Teko_Utilities.hpp"
63
64namespace Teko {
65namespace TpetraHelpers {
66
67using Teuchos::RCP;
68using Teuchos::rcp;
69using Teuchos::rcp_dynamic_cast;
70
71BlockedTpetraOperator::BlockedTpetraOperator(const std::vector<std::vector<GO> > & vars,
72 const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & content,
73 const std::string & label)
74 : Teko::TpetraHelpers::TpetraOperatorWrapper(), label_(label)
75{
76 SetContent(vars,content);
77}
78
79void BlockedTpetraOperator::SetContent(const std::vector<std::vector<GO> > & vars,
80 const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & content)
81{
82 fullContent_ = content;
83 blockedMapping_ = rcp(new TpetraBlockedMappingStrategy(vars,fullContent_->getDomainMap(),
84 *fullContent_->getDomainMap()->getComm()));
85 SetMapStrategy(blockedMapping_);
86
87 // build thyra operator
88 BuildBlockedOperator();
89}
90
91void BlockedTpetraOperator::BuildBlockedOperator()
92{
93 TEUCHOS_ASSERT(blockedMapping_!=Teuchos::null);
94
95 // get a CRS matrix
96 const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsContent
97 = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(fullContent_);
98
99 // ask the strategy to build the Thyra operator for you
100 if(blockedOperator_==Teuchos::null) {
101 blockedOperator_ = blockedMapping_->buildBlockedThyraOp(crsContent,label_);
102 }
103 else {
104 const RCP<Thyra::BlockedLinearOpBase<ST> > blkOp
105 = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<ST> >(blockedOperator_,true);
106 blockedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
107 }
108
109 // set whatever is returned
110 SetOperator(blockedOperator_,false);
111
112 // reorder if neccessary
113 if(reorderManager_!=Teuchos::null)
114 Reorder(*reorderManager_);
115}
116
117const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > BlockedTpetraOperator::GetBlock(int i,int j) const
118{
119 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
120 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());
121
122 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blkOp->getBlock(i,j),true);
123 return tOp->getConstTpetraOperator();
124}
125
130{
131 reorderManager_ = rcp(new BlockReorderManager(brm));
132
133 // build reordered objects
134 RCP<const MappingStrategy> reorderMapping = rcp(new TpetraReorderedMappingStrategy(*reorderManager_,blockedMapping_));
135 RCP<const Thyra::BlockedLinearOpBase<ST> > blockOp
136 = rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(blockedOperator_);
137
138 RCP<const Thyra::LinearOpBase<ST> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
139
140 // set them as working values
141 SetMapStrategy(reorderMapping);
142 SetOperator(A,false);
143}
144
147{
148 SetMapStrategy(blockedMapping_);
149 SetOperator(blockedOperator_,false);
150 reorderManager_ = Teuchos::null;
151}
152
155void BlockedTpetraOperator::WriteBlocks(const std::string & prefix) const
156{
157 RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > blockOp
158 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<ST> >(blockedOperator_);
159
160 // get size of blocked block operator
161 int rows = Teko::blockRowCount(blockOp);
162
163 for(int i=0;i<rows;i++) {
164 for(int j=0;j<rows;j++) {
165 // build the file name
166 std::stringstream ss;
167 ss << prefix << "_" << i << j << ".mm";
168
169 // get the row matrix object
170 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blockOp->getBlock(i,j));
171 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > mat
172 = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator());
173
174 // write to file
175 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<ST,LO,GO,NT> >::writeSparseFile(ss.str().c_str(),mat);
176 }
177 }
178}
179
181{
182 Tpetra::Vector<ST,LO,GO,NT> xf(getRangeMap());
183 Tpetra::Vector<ST,LO,GO,NT> xs(getRangeMap());
184 Tpetra::Vector<ST,LO,GO,NT> y(getDomainMap());
185
186 // test operator many times
187 bool result = true;
188 ST 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.randomize();
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 diffNorm = xs.norm2();
201 trueNorm = xf.norm2();
202
203 // check result
204 result &= (diffNorm/trueNorm < tol);
205 }
206
207 return result;
208}
209
210} // end namespace TpetraHelpers
211} // 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< GO > > &vars, const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &content)
virtual void WriteBlocks(const std::string &prefix) const
BlockedTpetraOperator(const std::vector< std::vector< GO > > &vars, const Teuchos::RCP< const Tpetra::Operator< ST, LO, GO, NT > > &content, const std::string &label="<ANYM>")
bool testAgainstFullOperator(int count, ST tol) const
Helps perform sanity checks.
void RemoveReording()
Remove any reordering on this object.
Implements the Epetra_Operator interface with a Thyra LinearOperator. This enables the use of absrtac...
const RCP< const Thyra::LinearOpBase< ST > > getThyraOp() const
Return the thyra operator associated with this wrapper.