Teko  Version of the Day
 All Classes Files Functions Variables Pages
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 
65 namespace Teko {
66 namespace TpetraHelpers {
67 
68 using Teuchos::RCP;
69 using Teuchos::rcp;
70 using Teuchos::rcp_dynamic_cast;
71 
72 StridedTpetraOperator::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 
84 StridedTpetraOperator::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 
91 void 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 
102 void 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 
126 const 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 
138 void 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 
155 void StridedTpetraOperator::RemoveReording()
156 {
157  SetMapStrategy(stridedMapping_);
158  SetOperator(stridedOperator_,false);
159  reorderManager_ = Teuchos::null;
160 }
161 
164 void 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 
196 std::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 
238 bool 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