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