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