EpetraExt_PutBlockMap.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #include "EpetraExt_PutBlockMap.h"
00030 #include "Epetra_Comm.h"
00031 #include "Epetra_BlockMap.h"
00032 #include "Epetra_Map.h"
00033 #include "Epetra_IntVector.h"
00034 #include "Epetra_IntSerialDenseVector.h"
00035 #include "Epetra_Import.h"
00036 
00037 using namespace Matlab;
00038 namespace Matlab {
00039 
00040 int CopyBlockMap(mxArray* matlabA, const Epetra_BlockMap& map) {
00041 
00042   int valueCount = 0;
00043   const Epetra_Comm & comm = map.Comm();
00044   int numProc = comm.NumProc();
00045   bool doSizes = !map.ConstantElementSize();
00046 
00047   if (numProc==1) {
00048     int * myElements = map.MyGlobalElements();
00049     int * elementSizeList = 0;
00050     if (doSizes) elementSizeList = map.ElementSizeList();
00051     return(DoCopyBlockMap(matlabA, valueCount, map.NumGlobalElements(), myElements, elementSizeList, doSizes));
00052   }
00053 
00054   int numRows = map.NumMyElements();
00055   
00056   Epetra_Map allGidsMap(-1, numRows, 0,comm);
00057   
00058   Epetra_IntVector allGids(allGidsMap);
00059   for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
00060   
00061   Epetra_IntVector allSizes(allGidsMap);
00062   for (int i=0; i<numRows; i++) allSizes[i] = map.ElementSize(i);
00063   
00064   // Now construct a Map on PE 0 by strip-mining the rows of the input matrix map.
00065   int numChunks = numProc;
00066   int stripSize = allGids.GlobalLength()/numChunks;
00067   int remainder = allGids.GlobalLength()%numChunks;
00068   int curStart = 0;
00069   int curStripSize = 0;
00070   Epetra_IntSerialDenseVector importGidList;
00071   Epetra_IntSerialDenseVector importSizeList;
00072   int numImportGids = 0;
00073   if (comm.MyPID()==0) {
00074     importGidList.Size(stripSize+1); // Set size of vector to max needed
00075     if (doSizes) importSizeList.Size(stripSize+1); // Set size of vector to max needed
00076   }
00077   for (int i=0; i<numChunks; i++) {
00078     if (comm.MyPID()==0) { // Only PE 0 does this part
00079       curStripSize = stripSize;
00080       if (i<remainder) curStripSize++; // handle leftovers
00081       for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
00082       curStart += curStripSize;
00083     }
00084     // The following import map will be non-trivial only on PE 0.
00085     Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
00086     Epetra_Import gidImporter(importGidMap, allGidsMap);
00087     
00088     Epetra_IntVector importGids(importGidMap);
00089     if (importGids.Import(allGids, gidImporter, Insert)) return(-1); 
00090     Epetra_IntVector importSizes(importGidMap);
00091     if (doSizes) if (importSizes.Import(allSizes, gidImporter, Insert)) return(-1); 
00092     
00093     // importGids (and importSizes, if non-trivial block map)
00094     // now have a list of GIDs (and sizes, respectively) for the current strip of map.
00095     
00096     int * myElements = importGids.Values();
00097     int * elementSizeList = 0;
00098     if (doSizes) elementSizeList = importSizes.Values();
00099     // Finally we are ready to write this strip of the map to file
00100     if (comm.MyPID()==0) {
00101       DoCopyBlockMap(matlabA, valueCount, importGids.MyLength(), myElements, elementSizeList, doSizes);
00102     }
00103   }
00104   return(0);
00105 }
00106 
00107 int DoCopyBlockMap(mxArray* matlabA, int& valueCount, int length, const int * v1, const int * v2, bool doSizes) {
00108   // declare and get initial values of all matlabA pointers
00109   double* matlabAvaluesPtr = mxGetPr(matlabA);
00110   int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
00111   int* matlabArowIndicesPtr = mxGetIr(matlabA);
00112 
00113   // set all matlabA pointers to the proper offset
00114   matlabAvaluesPtr += valueCount;
00115   matlabArowIndicesPtr += valueCount;
00116   int numGidsDone = valueCount;
00117   if (doSizes) {
00118     numGidsDone /= 2;
00119   }
00120   
00121   matlabAcolumnIndicesPtr += numGidsDone;
00122   
00123   for (int i=0; i<length; i++) {
00124     *matlabAcolumnIndicesPtr++ = valueCount;
00125     *matlabArowIndicesPtr++ = 0;
00126     *matlabAvaluesPtr++ = v1[i]; // row GID of block map
00127 
00128     if (doSizes) {
00129       *matlabAvaluesPtr++ = v2[i]; // size of block
00130       valueCount += 2;
00131       *matlabArowIndicesPtr++ = 1;
00132     }
00133     else {
00134       valueCount++;
00135     }
00136   }
00137   
00138   return(0);
00139 }
00140 } // namespace Matlab

Generated on Tue Jul 13 09:23:06 2010 for EpetraExt by  doxygen 1.4.7