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