Teko Version of the Day
Loading...
Searching...
No Matches
Teko_TpetraThyraConverter.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_TpetraThyraConverter.hpp"
48#include "Tpetra_Core.hpp"
49
50// Teuchos includes
51#include "Teuchos_Array.hpp"
52#include "Teuchos_ArrayRCP.hpp"
53
54// Thyra includes
55#include "Thyra_DefaultProductVectorSpace.hpp"
56#include "Thyra_DefaultProductMultiVector.hpp"
57#include "Thyra_SpmdMultiVectorBase.hpp"
58#include "Thyra_SpmdVectorSpaceBase.hpp"
59#include "Thyra_MultiVectorStdOps.hpp"
60
61#include <iostream>
62#include <vector>
63
64using Teuchos::RCP;
65using Teuchos::Ptr;
66using Teuchos::rcp;
67using Teuchos::rcpFromRef;
68using Teuchos::rcp_dynamic_cast;
69using Teuchos::ptr_dynamic_cast;
70using Teuchos::null;
71
72namespace Teko {
73namespace TpetraHelpers {
74
75// const Teuchos::RCP<const Thyra::MultiVectorBase<double> >
76// blockTpetraToThyra(int numVectors,const double * tpetraData,int leadingDim,const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > & vs,int & localDim)
77
78void blockTpetraToThyra(int numVectors,Teuchos::ArrayRCP<const ST> tpetraData,int leadingDim,const Teuchos::Ptr<Thyra::MultiVectorBase<ST> > & mv,int & localDim)
79{
80 localDim = 0;
81
82 // check the base case
83 const Ptr<Thyra::ProductMultiVectorBase<ST> > prodMV
84 = ptr_dynamic_cast<Thyra::ProductMultiVectorBase<ST> > (mv);
85 if(prodMV==Teuchos::null) {
86 // VS object must be a SpmdMultiVector object
87 const Ptr<Thyra::SpmdMultiVectorBase<ST> > spmdX = ptr_dynamic_cast<Thyra::SpmdMultiVectorBase<ST> >(mv,true);
88 const RCP<const Thyra::SpmdVectorSpaceBase<ST> > spmdVS = spmdX->spmdSpace();
89
90 int localSubDim = spmdVS->localSubDim();
91
92 Thyra::Ordinal thyraLeadingDim=0;
93
94 Teuchos::ArrayRCP<ST> thyraData_arcp;
95 Teuchos::ArrayView<ST> thyraData;
96 spmdX->getNonconstLocalData(Teuchos::outArg(thyraData_arcp),Teuchos::outArg(thyraLeadingDim));
97 thyraData = thyraData_arcp(); // build array view
98
99 for(int i=0;i<localSubDim;i++) {
100 // copy each vector
101 for(int v=0;v<numVectors;v++)
102 thyraData[i+thyraLeadingDim*v] = tpetraData[i+leadingDim*v];
103 }
104
105 // set the local dimension
106 localDim = localSubDim;
107
108 return;
109 }
110
111 // this keeps track of current location in the tpetraData vector
112 Teuchos::ArrayRCP<const ST> localData = tpetraData;
113
114 // loop over all the blocks in the vector space
115 for(int blkIndex=0;blkIndex<prodMV->productSpace()->numBlocks();blkIndex++) {
116 int subDim = 0;
117 const RCP<Thyra::MultiVectorBase<ST> > blockVec = prodMV->getNonconstMultiVectorBlock(blkIndex);
118
119 // perform the recusive copy
120 blockTpetraToThyra(numVectors, localData,leadingDim,blockVec.ptr(),subDim);
121
122 // shift to the next block
123 localData += subDim;
124
125 // account for the size of this subblock
126 localDim += subDim;
127 }
128}
129
130// Convert a Tpetra_MultiVector with assumed block structure dictated by the
131// vector space into a Thyra::MultiVectorBase object.
132// const Teuchos::RCP<const Thyra::MultiVectorBase<double> > blockTpetraToThyra(const Tpetra_MultiVector & e,const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > & vs)
133void blockTpetraToThyra(const Tpetra::MultiVector<ST,LO,GO,NT> & tpetraX,const Teuchos::Ptr<Thyra::MultiVectorBase<ST> > & thyraX)
134{
135 TEUCHOS_ASSERT((Tpetra::global_size_t) thyraX->range()->dim()==tpetraX.getGlobalLength());
136
137 // extract local information from the Tpetra_MultiVector
138 LO leadingDim=0,localDim=0;
139 leadingDim = tpetraX.getStride();
140 Teuchos::ArrayRCP<const ST> tpetraData = tpetraX.get1dView();
141
142 int numVectors = tpetraX.getNumVectors();
143
144 blockTpetraToThyra(numVectors,tpetraData,leadingDim,thyraX.ptr(),localDim);
145
146 TEUCHOS_ASSERT((size_t) localDim==tpetraX.getLocalLength());
147}
148
149void blockThyraToTpetra(LO numVectors,Teuchos::ArrayRCP<ST> tpetraData,LO leadingDim,const Teuchos::RCP<const Thyra::MultiVectorBase<ST> > & tX,LO & localDim)
150{
151 localDim = 0;
152
153 // check the base case
154 const RCP<const Thyra::ProductMultiVectorBase<ST> > prodX
155 = rcp_dynamic_cast<const Thyra::ProductMultiVectorBase<ST> > (tX);
156 if(prodX==Teuchos::null) {
157 // the base case
158
159 // VS object must be a SpmdMultiVector object
160 RCP<const Thyra::SpmdMultiVectorBase<ST> > spmdX = rcp_dynamic_cast<const Thyra::SpmdMultiVectorBase<ST> >(tX,true);
161 RCP<const Thyra::SpmdVectorSpaceBase<ST> > spmdVS = spmdX->spmdSpace();
162
163 Thyra::Ordinal thyraLeadingDim=0;
164 Teuchos::ArrayView<const ST> thyraData;
165 Teuchos::ArrayRCP<const ST> thyraData_arcp;
166 spmdX->getLocalData(Teuchos::outArg(thyraData_arcp),Teuchos::outArg(thyraLeadingDim));
167 thyraData = thyraData_arcp(); // grab the array view
168
169 LO localSubDim = spmdVS->localSubDim();
170 for(LO i=0;i<localSubDim;i++) {
171 // copy each vector
172 for(LO v=0;v<numVectors;v++){
173 tpetraData[i+leadingDim*v] = thyraData[i+thyraLeadingDim*v];
174 }
175 }
176
177 // set the local dimension
178 localDim = localSubDim;
179
180 return;
181 }
182
183 const RCP<const Thyra::ProductVectorSpaceBase<ST> > prodVS = prodX->productSpace();
184
185 // this keeps track of current location in the tpetraData vector
186 Teuchos::ArrayRCP<ST> localData = tpetraData;
187
188 // loop over all the blocks in the vector space
189 for(int blkIndex=0;blkIndex<prodVS->numBlocks();blkIndex++) {
190 int subDim = 0;
191
192 // construct the block vector
193 blockThyraToTpetra(numVectors, localData,leadingDim,prodX->getMultiVectorBlock(blkIndex),subDim);
194
195 // shift to the next block
196 localData += subDim;
197
198 // account for the size of this subblock
199 localDim += subDim;
200 }
201
202 return;
203}
204
205// Convert a Thyra::MultiVectorBase object to a Tpetra_MultiVector object with
206// the map defined by the Tpetra_Map.
207// const Teuchos::RCP<const Tpetra_MultiVector>
208// blockThyraToTpetra(const Teuchos::RCP<const Thyra::MultiVectorBase<double> > & tX,const RCP<const Tpetra_Map> & map)
209void blockThyraToTpetra(const Teuchos::RCP<const Thyra::MultiVectorBase<ST> > & thyraX,Tpetra::MultiVector<ST,LO,GO,NT> & tpetraX)
210{
211 // build an Tpetra_MultiVector object
212 LO numVectors = thyraX->domain()->dim();
213
214 // make sure the number of vectors are the same
215 TEUCHOS_ASSERT((size_t) numVectors==tpetraX.getNumVectors());
216 TEUCHOS_ASSERT((Tpetra::global_size_t) thyraX->range()->dim()==tpetraX.getGlobalLength());
217
218 // extract local information from the Tpetra_MultiVector
219 LO leadingDim=0,localDim=0;
220 leadingDim = tpetraX.getStride();
221 Teuchos::ArrayRCP<ST> tpetraData = tpetraX.get1dViewNonConst();
222
223 // perform recursive copy
224 blockThyraToTpetra(numVectors,tpetraData,leadingDim,thyraX,localDim);
225
226 // sanity check
227 TEUCHOS_ASSERT((size_t) localDim==tpetraX.getLocalLength());
228}
229
230void thyraVSToTpetraMap(std::vector<GO> & myIndicies, int blockOffset, const Thyra::VectorSpaceBase<ST> & vs, int & localDim)
231{
232 // zero out set local dimension
233 localDim = 0;
234
235 const RCP<const Thyra::ProductVectorSpaceBase<ST> > prodVS
236 = rcp_dynamic_cast<const Thyra::ProductVectorSpaceBase<ST> >(rcpFromRef(vs));
237
238 // is more recursion needed?
239 if(prodVS==Teuchos::null) {
240 // base case
241
242 // try to cast to an SPMD capable vector space
243 const RCP<const Thyra::SpmdVectorSpaceBase<ST> > spmdVS
244 = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<ST> >(rcpFromRef(vs));
245 TEUCHOS_TEST_FOR_EXCEPTION(spmdVS==Teuchos::null,std::runtime_error,
246 "thyraVSToTpetraMap requires all subblocks to be SPMD");
247
248 // get local data storage information
249 int localOffset = spmdVS->localOffset();
250 int localSubDim = spmdVS->localSubDim();
251
252 // add indicies to matrix
253 for(int i=0;i<localSubDim;i++)
254 myIndicies.push_back(blockOffset+localOffset+i);
255
256 localDim += localSubDim;
257
258 return;
259 }
260
261 // loop over all the blocks in the vector space
262 for(int blkIndex=0;blkIndex<prodVS->numBlocks();blkIndex++) {
263 int subDim = 0;
264
265 // construct the block vector
266 thyraVSToTpetraMap(myIndicies, blockOffset,*prodVS->getBlock(blkIndex),subDim);
267
268 blockOffset += prodVS->getBlock(blkIndex)->dim();
269
270 // account for the size of this subblock
271 localDim += subDim;
272 }
273}
274
275// From a Thyra vector space create a compatable Tpetra_Map
276const RCP<Tpetra::Map<LO,GO,NT> > thyraVSToTpetraMap(const Thyra::VectorSpaceBase<ST> & vs,const RCP<const Teuchos::Comm<Thyra::Ordinal> > & /* comm */)
277{
278 int localDim = 0;
279 std::vector<GO> myGIDs;
280
281 // call recursive routine that constructs the mapping
282 thyraVSToTpetraMap(myGIDs,0,vs,localDim);
283
284 TEUCHOS_ASSERT(myGIDs.size() == (size_t) localDim);
285
286 // FIXME (mfh 12 Jul 2018) This ignores the input comm, so it can't
287 // be right.
288
289 // create the map
290 return rcp(new Tpetra::Map<LO,GO,NT>(vs.dim(), Teuchos::ArrayView<const GO>(myGIDs), 0, Tpetra::getDefaultComm()));
291}
292
293} // end namespace Tpetra
294} // end namespace Teko