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