EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_Transpose_RowMatrix.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
44#include <Epetra_Export.h>
45#include <Epetra_CrsGraph.h>
46#include <Epetra_CrsMatrix.h>
47#include <Epetra_Map.h>
48#include <Epetra_Import.h>
49#include <Epetra_Export.h>
50#include <Epetra_Comm.h>
51
52#include <Teuchos_TimeMonitor.hpp>
53#include <vector>
54
55//#define ENABLE_TRANSPOSE_TIMINGS
56
57namespace EpetraExt {
58
59// Provide a "resize" operation for double*'s.
60inline void resize_doubles(int nold,int nnew,double*& d){
61 if(nnew > nold){
62 double *tmp = new double[nnew];
63 for(int i=0; i<nold; i++)
64 tmp[i]=d[i];
65 delete [] d;
66 d=tmp;
67 }
68}
69
70
73{
74 if( TransposeMatrix_ ) delete TransposeMatrix_;
75
76 if( !OrigMatrixIsCrsMatrix_ )
77 {
78 delete [] Indices_;
79 delete [] Values_;
80 }
81}
82
83//=========================================================================
85{
86#ifdef ENABLE_TRANSPOSE_TIMINGS
87 Teuchos::Time myTime("global");
88 Teuchos::TimeMonitor MM(myTime);
89 Teuchos::RCP<Teuchos::Time> mtime;
90 mtime=MM.getNewTimer("Transpose: CreateTransposeLocal 1");
91 mtime->start();
92#endif
93
94 int i,j,err;
95 const Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<const Epetra_CrsMatrix*>(&orig);
96 if(OrigCrsMatrix) OrigMatrixIsCrsMatrix_ = true;
97 else OrigMatrixIsCrsMatrix_ = false;
98
99 const Epetra_Map & TransMap = orig.RowMatrixColMap();
100 int TransNnz = orig.NumMyNonzeros();
101 int NumIndices;
102
103 Epetra_CrsMatrix *TempTransA1 = new Epetra_CrsMatrix(Copy, TransMap,orig.RowMatrixRowMap(),0);
104 Epetra_IntSerialDenseVector & TransRowptr = TempTransA1->ExpertExtractIndexOffset();
105 Epetra_IntSerialDenseVector & TransColind = TempTransA1->ExpertExtractIndices();
106 double *& TransVals = TempTransA1->ExpertExtractValues();
107 NumMyRows_ = orig.NumMyRows();
108 NumMyCols_ = orig.NumMyCols();
109
110 TransRowptr.Resize(NumMyCols_+1);
111 TransColind.Resize(TransNnz);
112 resize_doubles(0,TransNnz,TransVals);
113 std::vector<int> CurrentStart(NumMyCols_,0);
114
115 // Count up nnz using the Rowptr to count the number of non-nonzeros.
116 if (OrigMatrixIsCrsMatrix_)
117 {
118 const Epetra_CrsGraph & OrigGraph = OrigCrsMatrix->Graph(); // Get matrix graph
119
120 for (i=0; i<NumMyRows_; i++)
121 {
122 err = OrigGraph.ExtractMyRowView(i, NumIndices, Indices_); // Get view of ith row
123 if (err != 0) throw OrigGraph.ReportError("ExtractMyRowView failed",err);
124 for (j=0; j<NumIndices; j++) ++CurrentStart[Indices_[j]];
125 }
126 }
127 else // Original is not a CrsMatrix
128 {
129 MaxNumEntries_ = 0;
130 MaxNumEntries_ = orig.MaxNumEntries();
131 delete [] Indices_; delete [] Values_;
132 Indices_ = new int[MaxNumEntries_];
133 Values_ = new double[MaxNumEntries_];
134
135 for (i=0; i<NumMyRows_; i++)
136 {
137 err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
138 if (err != 0) {
139 std::cerr << "ExtractMyRowCopy failed."<<std::endl;
140 throw err;
141 }
142 for (j=0; j<NumIndices; j++) ++CurrentStart[Indices_[j]];
143 }
144 }
145
146 // Scansum the rowptr; reset currentstart
147 TransRowptr[0] = 0;
148 for (i=1;i<NumMyCols_+1; i++) TransRowptr[i] = CurrentStart[i-1] + TransRowptr[i-1];
149 for (i=0;i<NumMyCols_; i++) CurrentStart[i] = TransRowptr[i];
150
151 // Now copy values and global indices into newly create transpose storage
152 for (i=0; i<NumMyRows_; i++)
153 {
154 if (OrigMatrixIsCrsMatrix_)
155 err = OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_);
156 else
157 err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
158 if (err != 0) {
159 std::cerr << "ExtractMyRowCopy failed."<<std::endl;
160 throw err;
161 }
162
163 for (j=0; j<NumIndices; j++)
164 {
165 int idx = CurrentStart[Indices_[j]];
166 TransColind[idx] = i;
167 TransVals[idx] = Values_[j];
168 ++CurrentStart[Indices_[j]];// increment the counter into the local row
169 }
170 }
171
172#ifdef ENABLE_TRANSPOSE_TIMINGS
173 mtime->stop();
174 mtime=MM.getNewTimer("Transpose: CreateTransposeLocal 2");
175 mtime->start();
176#endif
177
178 // Prebuild the importers and exporters the no-communication way, flipping the importers
179 // and exporters around.
180 Epetra_Import * myimport = 0;
181 Epetra_Export * myexport = 0;
182 if(OrigMatrixIsCrsMatrix_ && OrigCrsMatrix->Importer())
183 myexport = new Epetra_Export(*OrigCrsMatrix->Importer());
184 if(OrigMatrixIsCrsMatrix_ && OrigCrsMatrix->Exporter())
185 myimport = new Epetra_Import(*OrigCrsMatrix->Exporter());
186
187#ifdef ENABLE_TRANSPOSE_TIMINGS
188 mtime->stop();
189 mtime=MM.getNewTimer("Transpose: CreateTransposeLocal 3");
190 mtime->start();
191#endif
192
193 // Call ExpertStaticFillComplete
194 err = TempTransA1->ExpertStaticFillComplete(orig.OperatorRangeMap(),orig.OperatorDomainMap(),myimport,myexport);
195 if (err != 0) {
196 throw TempTransA1->ReportError("ExpertStaticFillComplete failed.",err);
197 }
198
199#ifdef ENABLE_TRANSPOSE_TIMINGS
200 mtime->stop();
201#endif
202
203 return TempTransA1;
204}
205
206
207//=========================================================================
211{
212 origObj_ = &orig;
213
214 if( !TransposeRowMap_ )
215 {
216 if( IgnoreNonLocalCols_ )
217 TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorRangeMap()); // Should be replaced with refcount =
218 else
219 TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorDomainMap()); // Should be replaced with refcount =
220 }
221
222 NumMyRows_ = orig.NumMyRows();
223 NumMyCols_ = orig.NumMyCols();
224
225 // This routine will work for any RowMatrix object, but will attempt cast the matrix to a CrsMatrix if
226 // possible (because we can then use a View of the matrix and graph, which is much cheaper).
227
228 // First get the local indices to count how many nonzeros will be in the
229 // transpose graph on each processor
230 Epetra_CrsMatrix * TempTransA1 = CreateTransposeLocal(orig);
231
232 if(!TempTransA1->Exporter()) {
233 // Short Circuit: There is no need to make another matrix since TransposeRowMap_== TransMap
234 newObj_ = TransposeMatrix_ = TempTransA1;
235 return *newObj_;
236 }
237
238#ifdef ENABLE_TRANSPOSE_TIMINGS
239 Teuchos::Time myTime("global");
240 Teuchos::TimeMonitor MM(myTime);
241 Teuchos::RCP<Teuchos::Time> mtime;
242 mtime=MM.getNewTimer("Transpose: Final FusedExport");
243 mtime->start();
244#endif
245
246 // Now that transpose matrix with shared rows is entered, create a new matrix that will
247 // get the transpose with uniquely owned rows (using the same row distribution as A).
248 TransposeMatrix_ = new Epetra_CrsMatrix(*TempTransA1,*TempTransA1->Exporter(),NULL,0,TransposeRowMap_, false);
249
250#ifdef ENABLE_TRANSPOSE_TIMINGS
251 mtime->stop();
252#endif
253
254 newObj_ = TransposeMatrix_;
255 delete TempTransA1;
256 return *newObj_;
257}
258
259//=========================================================================
261{
262 Epetra_CrsMatrix * TempTransA1 = CreateTransposeLocal(*origObj_);
263 const Epetra_Export * TransposeExporter=0;
264 bool DeleteExporter = false;
265
266 if(TempTransA1->Exporter()) TransposeExporter = TempTransA1->Exporter();
267 else {
268 DeleteExporter=true;
269 TransposeExporter = new Epetra_Export(TransposeMatrix_->DomainMap(),TransposeMatrix_->RowMap());
270 }
271
272 TransposeMatrix_->PutScalar(0.0); // Zero out all values of the matrix
273
274 EPETRA_CHK_ERR(TransposeMatrix_->Export(*TempTransA1, *TransposeExporter, Add));
275
276 if(DeleteExporter) delete TransposeExporter;
277 delete TempTransA1;
278 return 0;
279}
280
281//=========================================================================
283{
284 EPETRA_CHK_ERR(-1); //Not Implemented Yet
285 return false;
286}
287
288} // namespace EpetraExt
Add
Copy
Epetra_CrsMatrix * CreateTransposeLocal(OriginalTypeRef orig)
Local-only transpose operator. Don't use this unless you're sure you know what you're doing.
NewTypeRef operator()(OriginalTypeRef orig)
Transpose Transform Operator.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
int Resize(int Length_in)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
void resize_doubles(int nold, int nnew, double *&d)