Teko Version of the Day
Loading...
Searching...
No Matches
Teko_TpetraBlockedMappingStrategy.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_TpetraBlockedMappingStrategy.hpp"
48#include "Teko_TpetraHelpers.hpp"
49
50#include "Thyra_TpetraThyraWrappers.hpp"
51#include "Thyra_TpetraLinearOp.hpp"
52#include "Thyra_DefaultProductMultiVector.hpp"
53#include "Thyra_DefaultProductVectorSpace.hpp"
54#include "Thyra_DefaultSpmdMultiVector.hpp"
55#include "Thyra_DefaultBlockedLinearOp.hpp"
56
57using Teuchos::RCP;
58using Teuchos::rcp;
59using Teuchos::rcp_dynamic_cast;
60
61namespace Teko {
62namespace TpetraHelpers {
63
64// Creates a strided mapping strategy. This class is useful
65// for breaking up nodally ordered matrices (i.e. the unknowns
66// in a FEM problem are ordered [u0,v0,p0,u1,v1,p1,...]). Current
67// implimentation only supports a fixed number of variables
68//
69// arguments:
70// vars - Number of different variables
71// map - original Tpetra::Map<LO,GO,NT> to be broken up
72// comm - Teuchos::RCP<Teuchos::Comm<int> > object related to the map
73//
74TpetraBlockedMappingStrategy::TpetraBlockedMappingStrategy(const std::vector<std::vector<GO> > & vars,
75 const Teuchos::RCP<const Tpetra::Map<LO,GO,NT> > & map, const Teuchos::Comm<int> & comm)
76{
77 rangeMap_ = map;
78 domainMap_ = map;
79 buildBlockTransferData(vars, rangeMap_,comm);
80}
81
82// Virtual function defined in MappingStrategy. This copies
83// an Tpetra::MultiVector<ST,LO,GO,NT> into a Thyra::MultiVectorBase with
84// blocking handled by the strides defined in the constructor.
85//
86// arguments:
87// X - source Tpetra::MultiVector<ST,LO,GO,NT>
88// thyra_X - destination Thyra::MultiVectorBase
89//
90void TpetraBlockedMappingStrategy::copyTpetraIntoThyra(const Tpetra::MultiVector<ST,LO,GO,NT>& X,
91 const Teuchos::Ptr<Thyra::MultiVectorBase<ST> > & thyra_X) const
92{
93 int count = X.getNumVectors();
94
95 std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > > subX;
96
97 // allocate vectors to copy into
98 Blocking::buildSubVectors(blockMaps_,subX,count);
99
100 // copy source vector to X vector
101 Blocking::one2many(subX,X,blockImport_);
102
103 // convert subX to an array of multi vectors
104 Teuchos::Array<RCP<Thyra::MultiVectorBase<ST> > > thyra_subX;
105 Teuchos::Ptr<Thyra::ProductMultiVectorBase<ST> > prod_X
106 = Teuchos::ptr_dynamic_cast<Thyra::ProductMultiVectorBase<ST> >(thyra_X);
107 for(unsigned int i=0;i<blockMaps_.size();i++) {
108 RCP<Thyra::TpetraMultiVector<ST,LO,GO,NT> > vec = rcp_dynamic_cast<Thyra::TpetraMultiVector<ST,LO,GO,NT> >(prod_X->getNonconstMultiVectorBlock(i),true);
109
110 fillDefaultSpmdMultiVector(vec,subX[i]);
111 }
112}
113
114// Virtual function defined in MappingStrategy. This copies
115// an Tpetra::MultiVector<ST,LO,GO,NT> into a Thyra::MultiVectorBase with
116// blocking handled by the strides defined in the constructor.
117//
118// arguments:
119// thyra_Y - source Thyra::MultiVectorBase
120// Y - destination Tpetra::MultiVector<ST,LO,GO,NT>
121//
122void TpetraBlockedMappingStrategy::copyThyraIntoTpetra(const RCP<const Thyra::MultiVectorBase<ST> > & thyra_Y,
123 Tpetra::MultiVector<ST,LO,GO,NT>& Y) const
124{
125 std::vector<RCP<const Tpetra::MultiVector<ST,LO,GO,NT> > > subY;
126 RCP<const Thyra::DefaultProductMultiVector<ST> > prod_Y
127 = rcp_dynamic_cast<const Thyra::DefaultProductMultiVector<ST> >(thyra_Y);
128
129 // convert thyra product vector to subY
130 for(unsigned int i=0;i<blockMaps_.size();i++){
131 RCP<const Thyra::TpetraMultiVector<ST,LO,GO,NT> > tmv = rcp_dynamic_cast<const Thyra::TpetraMultiVector<ST,LO,GO,NT> >(prod_Y->getMultiVectorBlock(i),true);
132 subY.push_back(tmv->getConstTpetraMultiVector());
133 }
134
135 // endow the subVectors with required information about the maps
136 // Blocking::associateSubVectors(blockMaps_,subY);
137
138 // copy solution vectors to Y vector
139 Blocking::many2one(Y,subY,blockExport_);
140}
141
142// this is the core routine that builds the maps
143// and importers/exporters neccessary for all the
144// transfers. Currently it simply calls out to the
145// interlaced tpetra functions. (Comment: this
146// routine should probably be private or protected
147// ... it is basically the meat of the constructor)
148//
149// arguments:
150// vars - Vector describing the blocking of variables
151// baseMap - basic map to use in the transfers
152// comm - Teuchos::RCP<Teuchos::Comm<int> > object
153//
154void TpetraBlockedMappingStrategy::buildBlockTransferData(const std::vector<std::vector<GO> > & vars,
155 const Teuchos::RCP<const Tpetra::Map<LO,GO,NT> > & baseMap, const Teuchos::Comm<int> & comm)
156{
157 // build block for each vector
158 for(std::size_t i=0;i<vars.size();i++) {
159 // build maps and exporters/importers
160 Blocking::MapPair mapPair = Blocking::buildSubMap(vars[i],comm);
161 Blocking::ImExPair iePair = Blocking::buildExportImport(*baseMap, mapPair);
162
163 blockMaps_.push_back(mapPair);
164 blockImport_.push_back(iePair.first);
165 blockExport_.push_back(iePair.second);
166 }
167}
168
169// Builds a blocked Thyra operator that uses the strided
170// mapping strategy to define sub blocks.
171//
172// arguments:
173// mat - Tpetra::CrsMatrix<ST,LO,GO,NT> with FillComplete called, this
174// matrix is assumed to be square, with the same
175// range and domain maps
176// returns: Blocked Thyra linear operator with sub blocks
177// defined by this mapping strategy
178//
179const Teuchos::RCP<Thyra::BlockedLinearOpBase<ST> >
180TpetraBlockedMappingStrategy::buildBlockedThyraOp(const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > & crsContent,const std::string & label) const
181{
182 int dim = blockMaps_.size();
183
184 RCP<Thyra::DefaultBlockedLinearOp<ST> > A = Thyra::defaultBlockedLinearOp<ST>();
185
186 A->beginBlockFill(dim,dim);
187 for(int i=0;i<dim;i++) {
188 for(int j=0;j<dim;j++) {
189 // label block correctly
190 std::stringstream ss;
191 ss << label << "_" << i << "," << j;
192
193 // build the blocks and place it the right location
194 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > blk = Blocking::buildSubBlock(i,j,crsContent,blockMaps_);
195 A->setNonconstBlock(i,j,Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(blk->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(blk->getDomainMap()),blk));
196 }
197 } // end for i
198 A->endBlockFill();
199
200 return A;
201}
202
203// Rebuilds a blocked Thyra operator that uses the strided
204// mapping strategy to define sub blocks.
205//
206// arguments:
207// crsContent - Tpetra::CrsMatrix<ST,LO,GO,NT> with FillComplete called, this
208// matrix is assumed to be square, with the same
209// range and domain maps
210// A - Destination block linear op composed of blocks of
211// Tpetra::CrsMatrix<ST,LO,GO,NT> at all relevant locations
212//
213void TpetraBlockedMappingStrategy::rebuildBlockedThyraOp(const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > & crsContent,
214 const RCP<Thyra::BlockedLinearOpBase<ST> > & A) const
215{
216 int dim = blockMaps_.size();
217
218 for(int i=0;i<dim;i++) {
219 for(int j=0;j<dim;j++) {
220 // get Tpetra version of desired block
221 RCP<Thyra::LinearOpBase<ST> > Aij = A->getNonconstBlock(i,j);
222 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tAij = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(Aij,true);
223 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > eAij = rcp_dynamic_cast<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tAij->getTpetraOperator(),true);
224
225 // rebuild the blocks and place it the right location
226 Blocking::rebuildSubBlock(i,j,crsContent,blockMaps_,*eAij);
227 }
228 } // end for i
229}
230
231} // end namespace Tpetra
232} // end namespace Teko