Teko Version of the Day
Loading...
Searching...
No Matches
Teko_StridedTpetraOperator.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_StridedTpetraOperator.hpp"
48#include "Teko_TpetraStridedMappingStrategy.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 "Tpetra_Vector.hpp"
61#include "MatrixMarket_Tpetra.hpp"
62
63#include "Teko_Utilities.hpp"
64
65namespace Teko {
66namespace TpetraHelpers {
67
68using Teuchos::RCP;
69using Teuchos::rcp;
70using Teuchos::rcp_dynamic_cast;
71
72StridedTpetraOperator::StridedTpetraOperator(int numVars,const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & content,
73 const std::string & label)
74 : Teko::TpetraHelpers::TpetraOperatorWrapper(), label_(label)
75{
76 std::vector<int> vars;
77
78 // build vector describing the sub maps
79 for(int i=0;i<numVars;i++) vars.push_back(1);
80
81 SetContent(vars,content);
82}
83
84StridedTpetraOperator::StridedTpetraOperator(const std::vector<int> & vars,const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & content,
85 const std::string & label)
86 : Teko::TpetraHelpers::TpetraOperatorWrapper(), label_(label)
87{
88 SetContent(vars,content);
89}
90
91void StridedTpetraOperator::SetContent(const std::vector<int> & vars,const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > & content)
92{
93 fullContent_ = content;
94 stridedMapping_ = rcp(new TpetraStridedMappingStrategy(vars,fullContent_->getDomainMap(),
95 *fullContent_->getDomainMap()->getComm()));
96 SetMapStrategy(stridedMapping_);
97
98 // build thyra operator
99 BuildBlockedOperator();
100}
101
102void StridedTpetraOperator::BuildBlockedOperator()
103{
104 TEUCHOS_ASSERT(stridedMapping_!=Teuchos::null);
105
106 // get a CRS matrix
107 const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > crsContent = rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(fullContent_);
108
109 // ask the strategy to build the Thyra operator for you
110 if(stridedOperator_==Teuchos::null) {
111 stridedOperator_ = stridedMapping_->buildBlockedThyraOp(crsContent,label_);
112 }
113 else {
114 const RCP<Thyra::BlockedLinearOpBase<ST> > blkOp = rcp_dynamic_cast<Thyra::BlockedLinearOpBase<ST> >(stridedOperator_,true);
115 stridedMapping_->rebuildBlockedThyraOp(crsContent,blkOp);
116 }
117
118 // set whatever is returned
119 SetOperator(stridedOperator_,false);
120
121 // reorder if neccessary
122 if(reorderManager_!=Teuchos::null)
123 Reorder(*reorderManager_);
124}
125
126const Teuchos::RCP<const Tpetra::Operator<ST,LO,GO,NT> > StridedTpetraOperator::GetBlock(int i,int j) const
127{
128 const RCP<const Thyra::BlockedLinearOpBase<ST> > blkOp
129 = Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(getThyraOp());
130
131 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blkOp->getBlock(i,j),true);
132 return tOp->getConstTpetraOperator();
133}
134
138void StridedTpetraOperator::Reorder(const BlockReorderManager & brm)
139{
140 reorderManager_ = rcp(new BlockReorderManager(brm));
141
142 // build reordered objects
143 RCP<const MappingStrategy> reorderMapping = rcp(new TpetraReorderedMappingStrategy(*reorderManager_,stridedMapping_));
144 RCP<const Thyra::BlockedLinearOpBase<ST> > blockOp
145 = rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<ST> >(stridedOperator_);
146
147 RCP<const Thyra::LinearOpBase<ST> > A = buildReorderedLinearOp(*reorderManager_,blockOp);
148
149 // set them as working values
150 SetMapStrategy(reorderMapping);
151 SetOperator(A,false);
152}
153
155void StridedTpetraOperator::RemoveReording()
156{
157 SetMapStrategy(stridedMapping_);
158 SetOperator(stridedOperator_,false);
159 reorderManager_ = Teuchos::null;
160}
161
164void StridedTpetraOperator::WriteBlocks(const std::string & prefix) const
165{
166 RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > blockOp
167 = rcp_dynamic_cast<Thyra::PhysicallyBlockedLinearOpBase<ST> >(stridedOperator_);
168
169 // get size of strided block operator
170 int rows = Teko::blockRowCount(blockOp);
171
172 for(int i=0;i<rows;i++) {
173 for(int j=0;j<rows;j++) {
174 // build the file name
175 std::stringstream ss;
176 ss << prefix << "_" << i << j << ".mm";
177
178 // get the row matrix object (Note: can't use "GetBlock" method b/c matrix might be reordered)
179 RCP<const Thyra::TpetraLinearOp<ST,LO,GO,NT> > tOp = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blockOp->getBlock(i,j));
180 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > mat
181 = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(tOp->getConstTpetraOperator());
182
183 // write to file
184 Tpetra::MatrixMarket::Writer<Tpetra::CrsMatrix<ST,LO,GO,NT> >::writeSparseFile(ss.str().c_str(),mat);
185 }
186 }
187}
188
196std::string StridedTpetraOperator::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 Thyra::TpetraLinearOp<ST,LO,GO,NT> > Aij = rcp_dynamic_cast<const Thyra::TpetraLinearOp<ST,LO,GO,NT> >(blockOp->getBlock(row,col),true);
211 RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > mat = Teuchos::rcp_dynamic_cast<const Tpetra::CrsMatrix<ST,LO,GO,NT> >(Aij->getConstTpetraOperator(),true);
212
213 // compute the norm
214 ST 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->getFrobeniusNorm();
224 break;
225 default:
226 TEUCHOS_TEST_FOR_EXCEPTION(true,std::runtime_error,
227 "StridedTpetraOperator::NormType incorrectly specified. Only Frobenius norm implemented for Tpetra matrices.");
228 }
229
230 ss << norm << " ";
231 }
232 ss << newline;
233 }
234
235 return ss.str();
236}
237
238bool StridedTpetraOperator::testAgainstFullOperator(int count,ST tol) const
239{
240 Tpetra::Vector<ST,LO,GO,NT> xf(getRangeMap());
241 Tpetra::Vector<ST,LO,GO,NT> xs(getRangeMap());
242 Tpetra::Vector<ST,LO,GO,NT> y(getDomainMap());
243
244 // test operator many times
245 bool result = true;
246 ST 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.randomize();
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 diffNorm = xs.norm2();
259 trueNorm = xf.norm2();
260
261 // check result
262 result &= (diffNorm/trueNorm < tol);
263 }
264
265 return result;
266}
267
268} // end namespace TpetraHelpers
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.