EpetraExt Development
EpetraExt_PutBlockMap.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "EpetraExt_PutBlockMap.h"
00043 #include "Epetra_Comm.h"
00044 #include "Epetra_BlockMap.h"
00045 #include "Epetra_Map.h"
00046 #include "Epetra_IntVector.h"
00047 #include "Epetra_IntSerialDenseVector.h"
00048 #include "Epetra_Import.h"
00049 
00050 using namespace Matlab;
00051 namespace Matlab {
00052 
00053 int CopyBlockMap(mxArray* matlabA, const Epetra_BlockMap& map) {
00054 
00055   int valueCount = 0;
00056   const Epetra_Comm & comm = map.Comm();
00057   int numProc = comm.NumProc();
00058   bool doSizes = !map.ConstantElementSize();
00059 
00060   if (numProc==1) {
00061     int * myElements = map.MyGlobalElements();
00062     int * elementSizeList = 0;
00063     if (doSizes) elementSizeList = map.ElementSizeList();
00064     return(DoCopyBlockMap(matlabA, valueCount, map.NumGlobalElements(), myElements, elementSizeList, doSizes));
00065   }
00066 
00067   int numRows = map.NumMyElements();
00068   
00069   Epetra_Map allGidsMap(-1, numRows, 0,comm);
00070   
00071   Epetra_IntVector allGids(allGidsMap);
00072   for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
00073   
00074   Epetra_IntVector allSizes(allGidsMap);
00075   for (int i=0; i<numRows; i++) allSizes[i] = map.ElementSize(i);
00076   
00077   // Now construct a Map on PE 0 by strip-mining the rows of the input matrix map.
00078   int numChunks = numProc;
00079   int stripSize = allGids.GlobalLength()/numChunks;
00080   int remainder = allGids.GlobalLength()%numChunks;
00081   int curStart = 0;
00082   int curStripSize = 0;
00083   Epetra_IntSerialDenseVector importGidList;
00084   Epetra_IntSerialDenseVector importSizeList;
00085   int numImportGids = 0;
00086   if (comm.MyPID()==0) {
00087     importGidList.Size(stripSize+1); // Set size of vector to max needed
00088     if (doSizes) importSizeList.Size(stripSize+1); // Set size of vector to max needed
00089   }
00090   for (int i=0; i<numChunks; i++) {
00091     if (comm.MyPID()==0) { // Only PE 0 does this part
00092       curStripSize = stripSize;
00093       if (i<remainder) curStripSize++; // handle leftovers
00094       for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
00095       curStart += curStripSize;
00096     }
00097     // The following import map will be non-trivial only on PE 0.
00098     Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
00099     Epetra_Import gidImporter(importGidMap, allGidsMap);
00100     
00101     Epetra_IntVector importGids(importGidMap);
00102     if (importGids.Import(allGids, gidImporter, Insert)) return(-1); 
00103     Epetra_IntVector importSizes(importGidMap);
00104     if (doSizes) if (importSizes.Import(allSizes, gidImporter, Insert)) return(-1); 
00105     
00106     // importGids (and importSizes, if non-trivial block map)
00107     // now have a list of GIDs (and sizes, respectively) for the current strip of map.
00108     
00109     int * myElements = importGids.Values();
00110     int * elementSizeList = 0;
00111     if (doSizes) elementSizeList = importSizes.Values();
00112     // Finally we are ready to write this strip of the map to file
00113     if (comm.MyPID()==0) {
00114       DoCopyBlockMap(matlabA, valueCount, importGids.MyLength(), myElements, elementSizeList, doSizes);
00115     }
00116   }
00117   return(0);
00118 }
00119 
00120 int DoCopyBlockMap(mxArray* matlabA, int& valueCount, int length, const int * v1, const int * v2, bool doSizes) {
00121   // declare and get initial values of all matlabA pointers
00122   double* matlabAvaluesPtr = mxGetPr(matlabA);
00123   int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
00124   int* matlabArowIndicesPtr = mxGetIr(matlabA);
00125 
00126   // set all matlabA pointers to the proper offset
00127   matlabAvaluesPtr += valueCount;
00128   matlabArowIndicesPtr += valueCount;
00129   int numGidsDone = valueCount;
00130   if (doSizes) {
00131     numGidsDone /= 2;
00132   }
00133   
00134   matlabAcolumnIndicesPtr += numGidsDone;
00135   
00136   for (int i=0; i<length; i++) {
00137     *matlabAcolumnIndicesPtr++ = valueCount;
00138     *matlabArowIndicesPtr++ = 0;
00139     *matlabAvaluesPtr++ = v1[i]; // row GID of block map
00140 
00141     if (doSizes) {
00142       *matlabAvaluesPtr++ = v2[i]; // size of block
00143       valueCount += 2;
00144       *matlabArowIndicesPtr++ = 1;
00145     }
00146     else {
00147       valueCount++;
00148     }
00149   }
00150   
00151   return(0);
00152 }
00153 } // namespace Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines