EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_PutMultiVector.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 
43 #include "Epetra_Comm.h"
44 #include "Epetra_BlockMap.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_Vector.h"
47 #include "Epetra_IntVector.h"
49 #include "Epetra_Import.h"
50 
51 using namespace Matlab;
52 namespace Matlab {
53 int CopyMultiVector(double** matlabApr, const Epetra_MultiVector& A) {
54 
55  Epetra_BlockMap bmap = A.Map();
56  const Epetra_Comm & comm = bmap.Comm();
57  int numProc = comm.NumProc();
58 
59  if (numProc==1)
60  DoCopyMultiVector(matlabApr, A);
61  else {
62 
63  // In the case of more than one column in the multivector, and writing to MatrixMarket
64  // format, we call this function recursively, passing each vector of the multivector
65  // individually so that we can get all of it written to file before going on to the next
66  // multivector
67  if (A.NumVectors() > 1) {
68  for (int i=0; i < A.NumVectors(); i++)
69  if (CopyMultiVector(matlabApr, *(A(i)))) return(-1);
70  return(0);
71  }
72 
73  Epetra_Map map(-1, bmap.NumMyPoints(), 0, comm);
74  // Create a veiw of this multivector using a map (instead of block map)
75  Epetra_MultiVector A1(View, map, A.Pointers(), A.NumVectors());
76  int numRows = map.NumMyElements();
77 
78  Epetra_Map allGidsMap(-1, numRows, 0,comm);
79 
80  Epetra_IntVector allGids(allGidsMap);
81  for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
82 
83  // Now construct a MultiVector on PE 0 by strip-mining the rows of the input matrix A.
84  int numChunks = numProc;
85  int stripSize = allGids.GlobalLength()/numChunks;
86  int remainder = allGids.GlobalLength()%numChunks;
87  int curStart = 0;
88  int curStripSize = 0;
89  Epetra_IntSerialDenseVector importGidList;
90  int numImportGids = 0;
91  if (comm.MyPID()==0)
92  importGidList.Size(stripSize+1); // Set size of vector to max needed
93  for (int i=0; i<numChunks; i++) {
94  if (comm.MyPID()==0) { // Only PE 0 does this part
95  curStripSize = stripSize;
96  if (i<remainder) curStripSize++; // handle leftovers
97  for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
98  curStart += curStripSize;
99  }
100  // The following import map will be non-trivial only on PE 0.
101  Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
102  Epetra_Import gidImporter(importGidMap, allGidsMap);
103  Epetra_IntVector importGids(importGidMap);
104  if (importGids.Import(allGids, gidImporter, Insert)) return(-1);
105 
106  // importGids now has a list of GIDs for the current strip of matrix rows.
107  // Use these values to build another importer that will get rows of the matrix.
108 
109  // The following import map will be non-trivial only on PE 0.
110  Epetra_Map importMap(-1, importGids.MyLength(), importGids.Values(), 0, comm);
111  Epetra_Import importer(importMap, map);
112  Epetra_MultiVector importA(importMap, A1.NumVectors());
113  if (importA.Import(A1, importer, Insert)) return(-1);
114 
115  // Finally we are ready to write this strip of the matrix to ostream
116  if (DoCopyMultiVector(matlabApr, importA)) return(-1);
117  }
118  }
119  return(0);
120 }
121 
122 int DoCopyMultiVector(double** matlabApr, const Epetra_MultiVector& A) {
123 
124  int ierr = 0;
125  int length = A.GlobalLength();
126  int numVectors = A.NumVectors();
127  const Epetra_Comm & comm = A.Map().Comm();
128  if (comm.MyPID()!=0) {
129  if (A.MyLength()!=0) ierr = -1;
130  }
131  else {
132  if (length!=A.MyLength()) ierr = -1;
133  double* matlabAvalues = *matlabApr;
134  double* Aptr = A.Values();
135  memcpy((void *)matlabAvalues, (void *)Aptr, sizeof(*Aptr) * length * numVectors);
136  *matlabApr += length;
137  }
138  int ierrGlobal;
139  comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
140  return(ierrGlobal);
141 }
142 } // namespace Matlab
int CopyMultiVector(double **matlabApr, const Epetra_MultiVector &A)
int DoCopyMultiVector(double **matlabApr, const Epetra_MultiVector &A)
int Size(int Length_in)
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
int * Values() const
int GlobalLength() const
virtual int MyPID() const =0
virtual const Epetra_BlockMap & Map() const =0
const Epetra_Comm & Comm() const
virtual int NumProc() const =0
int NumMyPoints() const
int MyLength() const