EpetraExt Development
Loading...
Searching...
No Matches
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"
48#include "Epetra_IntSerialDenseVector.h"
49#include "Epetra_Import.h"
50
51using namespace Matlab;
52namespace Matlab {
53int 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
122int 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
Insert
View
virtual int NumProc() const=0
const Epetra_BlockMap & Map() const
int Size(int Length_in)
int GlobalLength() const
int NumVectors() const
int MyLength() const
double * Values() const
double ** Pointers() const
int DoCopyMultiVector(double **matlabApr, const Epetra_MultiVector &A)
int CopyMultiVector(double **matlabApr, const Epetra_MultiVector &A)