Teko Version of the Day
Loading...
Searching...
No Matches
Teko_TpetraStridedMappingStrategy.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_TpetraStridedMappingStrategy.hpp"
48#include "Teko_InterlacedTpetra.hpp"
49#include "Teko_TpetraHelpers.hpp"
50
51#include "Thyra_TpetraThyraWrappers.hpp"
52#include "Thyra_TpetraLinearOp.hpp"
53#include "Thyra_DefaultProductMultiVector.hpp"
54#include "Thyra_DefaultProductVectorSpace.hpp"
55#include "Thyra_DefaultSpmdMultiVector.hpp"
56#include "Thyra_DefaultBlockedLinearOp.hpp"
57
58using Teuchos::RCP;
59using Teuchos::rcp;
60using Teuchos::rcp_dynamic_cast;
61
62namespace Teko {
63namespace TpetraHelpers {
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 Tpetra::Map<LO,GO,NT> to be broken up
73// comm - Teuchos::Comm<int> object related to the map
74//
75TpetraStridedMappingStrategy::TpetraStridedMappingStrategy(const std::vector<int> & vars,const RCP<const Tpetra::Map<LO,GO,NT> > & map,
76 const Teuchos::Comm<int> & 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 Tpetra::MultiVector<ST,LO,GO,NT> into a Thyra::MultiVectorBase with
85// blocking handled by the strides defined in the constructor.
86//
87// arguments:
88// X - source Tpetra::MultiVector<ST,LO,GO,NT>
89// thyra_X - destination Thyra::MultiVectorBase
90//
91void TpetraStridedMappingStrategy::copyTpetraIntoThyra(const Tpetra::MultiVector<ST,LO,GO,NT>& X,
92 const Teuchos::Ptr<Thyra::MultiVectorBase<ST> > & thyra_X) const
93{
94 int count = X.getNumVectors();
95
96 std::vector<RCP<Tpetra::MultiVector<ST,LO,GO,NT> > > subX;
97
98 // allocate vectors to copy into
99 Strided::buildSubVectors(blockMaps_,subX,count);
100
101 // copy source vector to X vector
102 Strided::one2many(subX,X,blockImport_);
103
104 // convert subX to an array of multi vectors
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 TpetraStridedMappingStrategy::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 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 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::Comm<int> object
153//
154void TpetraStridedMappingStrategy::buildBlockTransferData(const std::vector<int> & vars,const Teuchos::RCP<const Tpetra::Map<LO,GO,NT> > & baseMap, const Teuchos::Comm<int> & 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 - Tpetra::CrsMatrix<ST,LO,GO,NT> 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<ST> >
172TpetraStridedMappingStrategy::buildBlockedThyraOp(const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > & crsContent,const std::string & label) const
173{
174 int dim = blockMaps_.size();
175
176 RCP<Thyra::DefaultBlockedLinearOp<ST> > A = Thyra::defaultBlockedLinearOp<ST>();
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 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > Aij = Strided::buildSubBlock(i,j,crsContent,blockMaps_);
187 A->setNonconstBlock(i,j,Thyra::tpetraLinearOp<ST,LO,GO,NT>(Thyra::tpetraVectorSpace<ST,LO,GO,NT>(Aij->getRangeMap()),Thyra::tpetraVectorSpace<ST,LO,GO,NT>(Aij->getDomainMap()),Aij));
188 }
189 } // end for i
190 A->endBlockFill();
191
192 return A;
193}
194
195// Rebuilds a blocked Thyra operator that uses the strided
196// mapping strategy to define sub blocks.
197//
198// arguments:
199// crsContent - Tpetra::CrsMatrix<ST,LO,GO,NT> with FillComplete called, this
200// matrix is assumed to be square, with the same
201// range and domain maps
202// A - Destination block linear op composed of blocks of
203// Tpetra::CrsMatrix<ST,LO,GO,NT> at all relevant locations
204//
205void TpetraStridedMappingStrategy::rebuildBlockedThyraOp(const RCP<const Tpetra::CrsMatrix<ST,LO,GO,NT> > & crsContent,
206 const RCP<Thyra::BlockedLinearOpBase<ST> > & A) const
207{
208 int dim = blockMaps_.size();
209
210 for(int i=0;i<dim;i++) {
211 for(int j=0;j<dim;j++) {
212 // get Tpetra version of desired block
213 RCP<Thyra::LinearOpBase<ST> > Aij = A->getNonconstBlock(i,j);
214 RCP<Thyra::TpetraLinearOp<ST,LO,GO,NT> > tAij = rcp_dynamic_cast<Thyra::TpetraLinearOp<ST,LO,GO,NT> >(Aij,true);
215 RCP<Tpetra::CrsMatrix<ST,LO,GO,NT> > eAij = rcp_dynamic_cast<Tpetra::CrsMatrix<ST,LO,GO,NT> >(tAij->getTpetraOperator(),true);
216
217 // rebuild the blocks and place it the right location
218 Strided::rebuildSubBlock(i,j,crsContent,blockMaps_,*eAij);
219 }
220 } // end for i
221}
222
223} // end namespace TpetraHelpers
224} // end namespace Teko