47 #include "Teko_EpetraThyraConverter.hpp"
50 #include "Teuchos_Array.hpp"
51 #include "Teuchos_ArrayRCP.hpp"
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"
65 using Teuchos::ptr_dynamic_cast;
68 using Teuchos::rcp_dynamic_cast;
69 using Teuchos::rcpFromRef;
78 void blockEpetraToThyra(
int numVectors,
const double* epetraData,
int leadingDim,
79 const Teuchos::Ptr<Thyra::MultiVectorBase<double> >& mv,
int& localDim) {
83 const Ptr<Thyra::ProductMultiVectorBase<double> > prodMV =
84 ptr_dynamic_cast<Thyra::ProductMultiVectorBase<double> >(mv);
85 if (prodMV == Teuchos::null) {
87 const Ptr<Thyra::SpmdMultiVectorBase<double> > spmdX =
88 ptr_dynamic_cast<Thyra::SpmdMultiVectorBase<double> >(mv,
true);
89 const RCP<const Thyra::SpmdVectorSpaceBase<double> > spmdVS = spmdX->spmdSpace();
91 int localSubDim = spmdVS->localSubDim();
93 Thyra::Ordinal thyraLeadingDim = 0;
97 Teuchos::ArrayRCP<double> thyraData_arcp;
98 Teuchos::ArrayView<double> thyraData;
99 spmdX->getNonconstLocalData(Teuchos::outArg(thyraData_arcp), Teuchos::outArg(thyraLeadingDim));
100 thyraData = thyraData_arcp();
102 for (
int i = 0; i < localSubDim; i++) {
104 for (
int v = 0; v < numVectors; v++)
105 thyraData[i + thyraLeadingDim * v] = epetraData[i + leadingDim * v];
111 localDim = localSubDim;
117 const double* localData = epetraData;
120 for (
int blkIndex = 0; blkIndex < prodMV->productSpace()->numBlocks(); blkIndex++) {
122 const RCP<Thyra::MultiVectorBase<double> > blockVec =
123 prodMV->getNonconstMultiVectorBlock(blkIndex);
126 blockEpetraToThyra(numVectors, localData, leadingDim, blockVec.ptr(), subDim);
140 void blockEpetraToThyra(
const Epetra_MultiVector& epetraX,
141 const Teuchos::Ptr<Thyra::MultiVectorBase<double> >& thyraX) {
142 TEUCHOS_ASSERT(thyraX->range()->dim() == epetraX.GlobalLength());
145 int leadingDim = 0, numVectors = 0, localDim = 0;
146 double* epetraData = 0;
147 epetraX.ExtractView(&epetraData, &leadingDim);
149 numVectors = epetraX.NumVectors();
151 blockEpetraToThyra(numVectors, epetraData, leadingDim, thyraX.ptr(), localDim);
153 TEUCHOS_ASSERT(localDim == epetraX.MyLength());
156 void blockThyraToEpetra(
int numVectors,
double* epetraData,
int leadingDim,
157 const Teuchos::RCP<
const Thyra::MultiVectorBase<double> >& tX,
162 const RCP<const Thyra::ProductMultiVectorBase<double> > prodX =
163 rcp_dynamic_cast<
const Thyra::ProductMultiVectorBase<double> >(tX);
164 if (prodX == Teuchos::null) {
168 RCP<const Thyra::SpmdMultiVectorBase<double> > spmdX =
169 rcp_dynamic_cast<
const Thyra::SpmdMultiVectorBase<double> >(tX,
true);
170 RCP<const Thyra::SpmdVectorSpaceBase<double> > spmdVS = spmdX->spmdSpace();
172 int localSubDim = spmdVS->localSubDim();
174 Thyra::Ordinal thyraLeadingDim = 0;
178 Teuchos::ArrayView<const double> thyraData;
179 Teuchos::ArrayRCP<const double> thyraData_arcp;
180 spmdX->getLocalData(Teuchos::outArg(thyraData_arcp), Teuchos::outArg(thyraLeadingDim));
181 thyraData = thyraData_arcp();
183 for (
int i = 0; i < localSubDim; i++) {
185 for (
int v = 0; v < numVectors; v++)
186 epetraData[i + leadingDim * v] = thyraData[i + thyraLeadingDim * v];
190 localDim = localSubDim;
195 const RCP<const Thyra::ProductVectorSpaceBase<double> > prodVS = prodX->productSpace();
198 double* localData = epetraData;
201 for (
int blkIndex = 0; blkIndex < prodVS->numBlocks(); blkIndex++) {
205 blockThyraToEpetra(numVectors, localData, leadingDim, prodX->getMultiVectorBlock(blkIndex),
223 void blockThyraToEpetra(
const Teuchos::RCP<
const Thyra::MultiVectorBase<double> >& thyraX,
224 Epetra_MultiVector& epetraX) {
226 int numVectors = thyraX->domain()->dim();
229 TEUCHOS_ASSERT(numVectors == epetraX.NumVectors());
230 TEUCHOS_ASSERT(thyraX->range()->dim() == epetraX.GlobalLength());
233 int leadingDim = 0, localDim = 0;
234 double* epetraData = 0;
235 epetraX.ExtractView(&epetraData, &leadingDim);
238 blockThyraToEpetra(numVectors, epetraData, leadingDim, thyraX, localDim);
241 TEUCHOS_ASSERT(localDim == epetraX.Map().NumMyElements());
244 void thyraVSToEpetraMap(std::vector<int>& myIndicies,
int blockOffset,
245 const Thyra::VectorSpaceBase<double>& vs,
int& localDim) {
249 const RCP<const Thyra::ProductVectorSpaceBase<double> > prodVS =
250 rcp_dynamic_cast<
const Thyra::ProductVectorSpaceBase<double> >(rcpFromRef(vs));
253 if (prodVS == Teuchos::null) {
257 const RCP<const Thyra::SpmdVectorSpaceBase<double> > spmdVS =
258 rcp_dynamic_cast<
const Thyra::SpmdVectorSpaceBase<double> >(rcpFromRef(vs));
259 TEUCHOS_TEST_FOR_EXCEPTION(spmdVS == Teuchos::null, std::runtime_error,
260 "thyraVSToEpetraMap requires all subblocks to be SPMD");
263 int localOffset = spmdVS->localOffset();
264 int localSubDim = spmdVS->localSubDim();
267 for (
int i = 0; i < localSubDim; i++) myIndicies.push_back(blockOffset + localOffset + i);
269 localDim += localSubDim;
275 for (
int blkIndex = 0; blkIndex < prodVS->numBlocks(); blkIndex++) {
279 thyraVSToEpetraMap(myIndicies, blockOffset, *prodVS->getBlock(blkIndex), subDim);
281 blockOffset += prodVS->getBlock(blkIndex)->dim();
289 const RCP<Epetra_Map> thyraVSToEpetraMap(
const Thyra::VectorSpaceBase<double>& vs,
290 const RCP<const Epetra_Comm>& comm) {
292 std::vector<int> myGIDs;
295 thyraVSToEpetraMap(myGIDs, 0, vs, localDim);
297 TEUCHOS_ASSERT(myGIDs.size() == (
unsigned int)localDim);
300 return rcp(
new Epetra_Map(vs.dim(), myGIDs.size(), &(myGIDs[0]), 0, *comm));