Teko Version of the Day
Teko_EpetraThyraConverter.cpp
00001 /*
00002 // @HEADER
00003 // 
00004 // ***********************************************************************
00005 // 
00006 //      Teko: A package for block and physics based preconditioning
00007 //                  Copyright 2010 Sandia Corporation 
00008 //  
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //  
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //  
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //  
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //  
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission. 
00026 //  
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //  
00039 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
00040 // 
00041 // ***********************************************************************
00042 // 
00043 // @HEADER
00044 
00045 */
00046 
00047 #include "Teko_EpetraThyraConverter.hpp"
00048 
00049 // Teuchos includes
00050 #include "Teuchos_Array.hpp"
00051 #include "Teuchos_ArrayRCP.hpp"
00052 
00053 // Thyra includes
00054 #include "Thyra_DefaultProductVectorSpace.hpp"
00055 #include "Thyra_DefaultProductMultiVector.hpp"
00056 #include "Thyra_SpmdMultiVectorBase.hpp"
00057 #include "Thyra_SpmdVectorSpaceBase.hpp"
00058 #include "Thyra_MultiVectorStdOps.hpp"
00059 
00060 #include <iostream>
00061 #include <vector>
00062 
00063 using Teuchos::RCP;
00064 using Teuchos::Ptr;
00065 using Teuchos::rcp;
00066 using Teuchos::rcpFromRef;
00067 using Teuchos::rcp_dynamic_cast;
00068 using Teuchos::ptr_dynamic_cast;
00069 using Teuchos::null;
00070 
00071 namespace Teko {
00072 namespace Epetra {
00073 
00074 // const Teuchos::RCP<const Thyra::MultiVectorBase<double> > 
00075 // blockEpetraToThyra(int numVectors,const double * epetraData,int leadingDim,const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > & vs,int & localDim)
00076 
00077 void blockEpetraToThyra(int numVectors,const double * epetraData,int leadingDim,const Teuchos::Ptr<Thyra::MultiVectorBase<double> > & mv,int & localDim)
00078 {
00079    localDim = 0;
00080 
00081    // check the base case
00082    const Ptr<Thyra::ProductMultiVectorBase<double> > prodMV 
00083          = ptr_dynamic_cast<Thyra::ProductMultiVectorBase<double> > (mv);
00084    if(prodMV==Teuchos::null) {
00085       // VS object must be a SpmdMultiVector object
00086       const Ptr<Thyra::SpmdMultiVectorBase<double> > spmdX  = ptr_dynamic_cast<Thyra::SpmdMultiVectorBase<double> >(mv,true);
00087       const RCP<const Thyra::SpmdVectorSpaceBase<double> > spmdVS = spmdX->spmdSpace();
00088 
00089       int localSubDim = spmdVS->localSubDim();
00090 
00091       Thyra::Ordinal thyraLeadingDim=0;
00092       // double * thyraData=0;
00093       // spmdX->getLocalData(&thyraData,&thyraLeadingDim);
00094 
00095       Teuchos::ArrayRCP<double> thyraData_arcp;
00096       Teuchos::ArrayView<double> thyraData;
00097       spmdX->getNonconstLocalData(Teuchos::outArg(thyraData_arcp),Teuchos::outArg(thyraLeadingDim));
00098       thyraData = thyraData_arcp(); // build array view
00099 
00100       for(int i=0;i<localSubDim;i++) {
00101          // copy each vector
00102          for(int v=0;v<numVectors;v++)
00103             thyraData[i+thyraLeadingDim*v] = epetraData[i+leadingDim*v];
00104       }
00105 
00106       // spmdX->commitLocalData(&thyraData[0]);
00107 
00108       // set the local dimension
00109       localDim = localSubDim;
00110 
00111       return;
00112    }
00113 
00114    // this keeps track of current location in the epetraData vector
00115    const double * localData = epetraData;
00116 
00117    // loop over all the blocks in the vector space
00118    for(int blkIndex=0;blkIndex<prodMV->productSpace()->numBlocks();blkIndex++) {
00119       int subDim = 0;
00120       const RCP<Thyra::MultiVectorBase<double> > blockVec = prodMV->getNonconstMultiVectorBlock(blkIndex);
00121 
00122       // perorm the recusive copy
00123       blockEpetraToThyra(numVectors, localData,leadingDim,blockVec.ptr(),subDim);
00124 
00125       // shift to the next block
00126       localData += subDim;
00127 
00128       // account for the size of this subblock
00129       localDim += subDim;
00130    }
00131 }
00132 
00133 // Convert a Epetra_MultiVector with assumed block structure dictated by the
00134 // vector space into a Thyra::MultiVectorBase object.
00135 // const Teuchos::RCP<const Thyra::MultiVectorBase<double> > blockEpetraToThyra(const Epetra_MultiVector & e,const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > & vs)
00136 void blockEpetraToThyra(const Epetra_MultiVector & epetraX,const Teuchos::Ptr<Thyra::MultiVectorBase<double> > & thyraX) 
00137 {
00138    TEUCHOS_ASSERT(thyraX->range()->dim()==epetraX.GlobalLength());
00139 
00140    // extract local information from the Epetra_MultiVector
00141    int leadingDim=0,numVectors=0,localDim=0;
00142    double * epetraData=0;
00143    epetraX.ExtractView(&epetraData,&leadingDim);
00144 
00145    numVectors = epetraX.NumVectors();
00146 
00147    blockEpetraToThyra(numVectors,epetraData,leadingDim,thyraX.ptr(),localDim);   
00148 
00149    TEUCHOS_ASSERT(localDim==epetraX.MyLength());
00150 }
00151 
00152 void blockThyraToEpetra(int numVectors,double * epetraData,int leadingDim,const Teuchos::RCP<const Thyra::MultiVectorBase<double> > & tX,int & localDim)
00153 {
00154    localDim = 0;
00155 
00156    // check the base case
00157    const RCP<const Thyra::ProductMultiVectorBase<double> > prodX
00158          = rcp_dynamic_cast<const Thyra::ProductMultiVectorBase<double> > (tX);
00159    if(prodX==Teuchos::null) {
00160       // the base case
00161 
00162       // VS object must be a SpmdMultiVector object
00163       RCP<const Thyra::SpmdMultiVectorBase<double> > spmdX  = rcp_dynamic_cast<const Thyra::SpmdMultiVectorBase<double> >(tX,true);
00164       RCP<const Thyra::SpmdVectorSpaceBase<double> > spmdVS = spmdX->spmdSpace();
00165 
00166       int localSubDim = spmdVS->localSubDim();
00167 
00168       Thyra::Ordinal thyraLeadingDim=0;
00169       // const double * thyraData=0;
00170       // spmdX->getLocalData(&thyraData,&thyraLeadingDim);
00171     
00172       Teuchos::ArrayView<const double> thyraData;
00173       Teuchos::ArrayRCP<const double> thyraData_arcp;
00174       spmdX->getLocalData(Teuchos::outArg(thyraData_arcp),Teuchos::outArg(thyraLeadingDim));
00175       thyraData = thyraData_arcp(); // grab the array view
00176 
00177       for(int i=0;i<localSubDim;i++) {
00178          // copy each vector
00179          for(int v=0;v<numVectors;v++)
00180             epetraData[i+leadingDim*v] = thyraData[i+thyraLeadingDim*v];
00181       }
00182 
00183       // set the local dimension
00184       localDim = localSubDim;
00185 
00186       return;
00187    }
00188 
00189    const RCP<const Thyra::ProductVectorSpaceBase<double> > prodVS = prodX->productSpace();
00190 
00191    // this keeps track of current location in the epetraData vector
00192    double * localData = epetraData;
00193 
00194    // loop over all the blocks in the vector space
00195    for(int blkIndex=0;blkIndex<prodVS->numBlocks();blkIndex++) {
00196       int subDim = 0;
00197 
00198       // construct the block vector
00199       blockThyraToEpetra(numVectors, localData,leadingDim,prodX->getMultiVectorBlock(blkIndex),subDim);
00200 
00201       // shift to the next block
00202       localData += subDim;
00203 
00204       // account for the size of this subblock
00205       localDim += subDim;
00206    }
00207 
00208    return;
00209 }
00210 
00211 // Convert a Thyra::MultiVectorBase object to a Epetra_MultiVector object with
00212 // the map defined by the Epetra_Map.
00213 // const Teuchos::RCP<const Epetra_MultiVector> 
00214 // blockThyraToEpetra(const Teuchos::RCP<const Thyra::MultiVectorBase<double> > & tX,const RCP<const Epetra_Map> & map)
00215 void blockThyraToEpetra(const Teuchos::RCP<const Thyra::MultiVectorBase<double> > & thyraX,Epetra_MultiVector & epetraX)
00216 {
00217    // build an Epetra_MultiVector object
00218    int numVectors = thyraX->domain()->dim();
00219 
00220    // make sure the number of vectors are the same
00221    TEUCHOS_ASSERT(numVectors==epetraX.NumVectors());
00222    TEUCHOS_ASSERT(thyraX->range()->dim()==epetraX.GlobalLength());
00223 
00224    // extract local information from the Epetra_MultiVector
00225    int leadingDim=0,localDim=0;
00226    double * epetraData=0;
00227    epetraX.ExtractView(&epetraData,&leadingDim);
00228 
00229    // perform recursive copy
00230    blockThyraToEpetra(numVectors,epetraData,leadingDim,thyraX,localDim);
00231 
00232    // sanity check
00233    TEUCHOS_ASSERT(localDim==epetraX.Map().NumMyElements());
00234 }
00235 
00236 void thyraVSToEpetraMap(std::vector<int> & myIndicies, int blockOffset, const Thyra::VectorSpaceBase<double> & vs, int & localDim)
00237 {
00238    // zero out set local dimension
00239    localDim = 0;
00240 
00241    const RCP<const Thyra::ProductVectorSpaceBase<double> > prodVS
00242          = rcp_dynamic_cast<const Thyra::ProductVectorSpaceBase<double> >(rcpFromRef(vs));
00243 
00244    // is more recursion needed?
00245    if(prodVS==Teuchos::null) {
00246       // base case
00247 
00248       // try to cast to an SPMD capable vector space
00249       const RCP<const Thyra::SpmdVectorSpaceBase<double> > spmdVS 
00250             = rcp_dynamic_cast<const Thyra::SpmdVectorSpaceBase<double> >(rcpFromRef(vs));
00251       TEST_FOR_EXCEPTION(spmdVS==Teuchos::null,std::runtime_error,
00252                          "thyraVSToEpetraMap requires all subblocks to be SPMD");
00253 
00254       // get local data storage information
00255       int localOffset = spmdVS->localOffset();
00256       int localSubDim = spmdVS->localSubDim();
00257 
00258       // add indicies to matrix
00259       for(int i=0;i<localSubDim;i++)
00260          myIndicies.push_back(blockOffset+localOffset+i);
00261 
00262       localDim += localSubDim;
00263 
00264       return;
00265    }
00266 
00267    // loop over all the blocks in the vector space
00268    for(int blkIndex=0;blkIndex<prodVS->numBlocks();blkIndex++) {
00269       int subDim = 0;
00270 
00271       // construct the block vector
00272       thyraVSToEpetraMap(myIndicies, blockOffset,*prodVS->getBlock(blkIndex),subDim);
00273 
00274       blockOffset += prodVS->getBlock(blkIndex)->dim();
00275 
00276       // account for the size of this subblock
00277       localDim += subDim;
00278    }
00279 }
00280 
00281 // From a Thyra vector space create a compatable Epetra_Map
00282 const RCP<Epetra_Map> thyraVSToEpetraMap(const Thyra::VectorSpaceBase<double> & vs,const RCP<const Epetra_Comm> & comm)
00283 {
00284    int localDim = 0;
00285    std::vector<int> myGIDs;
00286    
00287    // call recursive routine that constructs the mapping
00288    thyraVSToEpetraMap(myGIDs,0,vs,localDim);
00289 
00290    TEUCHOS_ASSERT(myGIDs.size()==(unsigned int) localDim);
00291 
00292    // create the map
00293    return rcp(new Epetra_Map(vs.dim(), myGIDs.size(), &(myGIDs[0]), 0, *comm));
00294 }
00295 
00296 } // end namespace Epetra
00297 } // end namespace Teko
 All Classes Files Functions Variables