EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_PutMultiVector.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //                     New_Package Example Package
00005 //                 Copyright (2004) 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_PutMultiVector.h"
00030 #include "Epetra_Comm.h"
00031 #include "Epetra_BlockMap.h"
00032 #include "Epetra_Map.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_IntVector.h"
00035 #include "Epetra_IntSerialDenseVector.h"
00036 #include "Epetra_Import.h"
00037 
00038 using namespace Matlab;
00039 namespace Matlab {
00040 int CopyMultiVector(double** matlabApr, const Epetra_MultiVector& A) {
00041 
00042   Epetra_BlockMap bmap = A.Map();
00043   const Epetra_Comm & comm = bmap.Comm();
00044   int numProc = comm.NumProc();
00045 
00046   if (numProc==1)
00047     DoCopyMultiVector(matlabApr, A);
00048   else {
00049 
00050     // In the case of more than one column in the multivector, and writing to MatrixMarket
00051     // format, we call this function recursively, passing each vector of the multivector
00052     // individually so that we can get all of it written to file before going on to the next 
00053     // multivector
00054     if (A.NumVectors() > 1) {
00055       for (int i=0; i < A.NumVectors(); i++)
00056   if (CopyMultiVector(matlabApr, *(A(i)))) return(-1);
00057       return(0);
00058     }
00059 
00060     Epetra_Map map(-1, bmap.NumMyPoints(), 0, comm);
00061     // Create a veiw of this multivector using a map (instead of block map)
00062     Epetra_MultiVector A1(View, map, A.Pointers(), A.NumVectors());
00063     int numRows = map.NumMyElements();
00064     
00065     Epetra_Map allGidsMap(-1, numRows, 0,comm);
00066     
00067     Epetra_IntVector allGids(allGidsMap);
00068     for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
00069     
00070     // Now construct a MultiVector on PE 0 by strip-mining the rows of the input matrix A.
00071     int numChunks = numProc;
00072     int stripSize = allGids.GlobalLength()/numChunks;
00073     int remainder = allGids.GlobalLength()%numChunks;
00074     int curStart = 0;
00075     int curStripSize = 0;
00076     Epetra_IntSerialDenseVector importGidList;
00077     int numImportGids = 0;
00078     if (comm.MyPID()==0) 
00079       importGidList.Size(stripSize+1); // Set size of vector to max needed
00080     for (int i=0; i<numChunks; i++) {
00081       if (comm.MyPID()==0) { // Only PE 0 does this part
00082   curStripSize = stripSize;
00083   if (i<remainder) curStripSize++; // handle leftovers
00084   for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
00085   curStart += curStripSize;
00086       }
00087       // The following import map will be non-trivial only on PE 0.
00088       Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
00089       Epetra_Import gidImporter(importGidMap, allGidsMap);
00090       Epetra_IntVector importGids(importGidMap);
00091       if (importGids.Import(allGids, gidImporter, Insert)) return(-1); 
00092 
00093       // importGids now has a list of GIDs for the current strip of matrix rows.
00094       // Use these values to build another importer that will get rows of the matrix.
00095 
00096       // The following import map will be non-trivial only on PE 0.
00097       Epetra_Map importMap(-1, importGids.MyLength(), importGids.Values(), 0, comm);
00098       Epetra_Import importer(importMap, map);
00099       Epetra_MultiVector importA(importMap, A1.NumVectors());
00100       if (importA.Import(A1, importer, Insert)) return(-1); 
00101 
00102       // Finally we are ready to write this strip of the matrix to ostream
00103       if (DoCopyMultiVector(matlabApr, importA)) return(-1);
00104     }
00105   }
00106   return(0);
00107 }
00108 
00109 int DoCopyMultiVector(double** matlabApr, const Epetra_MultiVector& A) {
00110 
00111   int ierr = 0;
00112   int length = A.GlobalLength();
00113   int numVectors = A.NumVectors();
00114   const Epetra_Comm & comm = A.Map().Comm();
00115   if (comm.MyPID()!=0) {
00116     if (A.MyLength()!=0) ierr = -1;
00117   }
00118   else {
00119     if (length!=A.MyLength()) ierr = -1;
00120     double* matlabAvalues = *matlabApr;
00121     double* Aptr = A.Values();
00122     memcpy((void *)matlabAvalues, (void *)Aptr, sizeof(*Aptr) * length * numVectors);
00123     *matlabApr += length;   
00124   }
00125   int ierrGlobal;
00126   comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
00127   return(ierrGlobal);
00128 }
00129 } // namespace Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines