EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_MMHelpers.h
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
42#ifndef EPETRAEXT_MMHELPERS_H
43#define EPETRAEXT_MMHELPERS_H
44
46#include "Epetra_ConfigDefs.h"
47#include "Epetra_DistObject.h"
48#include "Epetra_Map.h"
49#include "Teuchos_RCP.hpp"
50#include "Epetra_Comm.h"
51#include "Epetra_Import.h"
52#include "Epetra_CrsMatrix.h"
53#include <Teuchos_TimeMonitor.hpp>
54
55#include <vector>
56#include <set>
57#include <map>
58
59
61
62namespace EpetraExt {
63class LightweightCrsMatrix;
64
65//#define HAVE_EPETRAEXT_DEBUG // for extra sanity checks
66
67
68
69
70// ==============================================================
71//struct that holds views of the contents of a CrsMatrix. These
72//contents may be a mixture of local and remote rows of the
73//actual matrix.
75public:
77
78 virtual ~CrsMatrixStruct();
79
80 void deleteContents();
81
83
84 // The following class members get used in the transpose modes of the MMM
85 // but not in the A*B mode.
87 int** indices;
88 double** values;
89 bool* remote;
92
93 // Maps and matrices (all modes)
100
101 // The following class members are only used for A*B mode
102 std::vector<int> targetMapToOrigRow;
103 std::vector<int> targetMapToImportRow;
104};
105
107
108// ==============================================================
110 public:
111 virtual ~CrsWrapper(){}
112
113 virtual const Epetra_Map& RowMap() const = 0;
114
115 virtual bool Filled() = 0;
116
117#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
118 virtual int InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices) = 0;
119
120 virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices) = 0;
121#endif
122
123#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
124 virtual int InsertGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices) = 0;
125
126 virtual int SumIntoGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices) = 0;
127#endif
128};
129
130// ==============================================================
132 public:
135
136 const Epetra_Map& RowMap() const;
137
138 bool Filled();
139
140#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
141 int InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices);
142 int SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices);
143#endif
144
145#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
146 int InsertGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices);
147 int SumIntoGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices);
148#endif
149
150 private:
151 Epetra_CrsMatrix& ecrsmat_;
152};
153
154// ==============================================================
155template<typename int_type>
157 public:
159 virtual ~CrsWrapper_GraphBuilder();
160
161 const Epetra_Map& RowMap() const {return rowmap_; }
162
163 bool Filled();
164
165 int InsertGlobalValues(int_type GlobalRow, int NumEntries, double* Values, int_type* Indices);
166 int SumIntoGlobalValues(int_type GlobalRow, int NumEntries, double* Values, int_type* Indices);
167
168 std::map<int_type,std::set<int_type>*>& get_graph();
169
170 int get_max_row_length() { return max_row_length_; }
171
172 private:
173 std::map<int_type,std::set<int_type>*> graph_;
174 const Epetra_Map& rowmap_;
175 int max_row_length_;
176};
177
178// ==============================================================
179template<typename int_type>
180void insert_matrix_locations(CrsWrapper_GraphBuilder<int_type>& graphbuilder,
182
183template<typename int_type>
185 const std::vector<int_type>& proc_col_ranges,
186 std::vector<int_type>& send_rows,
187 std::vector<int>& rows_per_send_proc);
188
189#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
191 const std::vector<int>& proc_col_ranges,
192 std::vector<int>& send_rows,
193 std::vector<int>& rows_per_send_proc);
194#endif
195
196#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
198 const std::vector<long long>& proc_col_ranges,
199 std::vector<long long>& send_rows,
200 std::vector<int>& rows_per_send_proc);
201#endif
202
203template<typename int_type>
204std::pair<int_type,int_type> get_col_range(const Epetra_Map& emap);
205
206template<typename int_type>
207std::pair<int_type,int_type> get_col_range(const Epetra_CrsMatrix& mtx);
208
209// ==============================================================
211 friend class LightweightMap;
212 public:
215 long long IndexBase_;
216#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
217 std::vector<int> MyGlobalElements_int_;
218 Epetra_HashTable<int> * LIDHash_int_;
219#endif
220#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
221 std::vector<long long> MyGlobalElements_LL_;
222 Epetra_HashTable<long long> * LIDHash_LL_;
223#endif
224
225 // For "copy" constructor only...
227
228 };
229
231 public:
233#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
234 LightweightMap(int NumGlobalElements,int NumMyElements, const int * MyGlobalElements, int IndexBase, bool GenerateHash=true);
235#endif
236#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
237 LightweightMap(long long NumGlobalElements,int NumMyElements, const long long * MyGlobalElements, int IndexBase, bool GenerateHash=true);
238 LightweightMap(long long NumGlobalElements,int NumMyElements, const long long * MyGlobalElements, long long IndexBase, bool GenerateHash=true);
239#endif
240 LightweightMap(const Epetra_Map & Map);
241 LightweightMap(const LightweightMap & Map);
243
245#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
246 int LID(int GID) const;
247 int GID(int LID) const;
248#endif
249
250#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
251 int LID(long long GID) const;
252#endif
253 long long GID64(int LID) const;
254 int NumMyElements() const;
255
256#if defined(EPETRA_NO_32BIT_GLOBAL_INDICES) && defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
257 // default implementation so that no compiler/linker error in case neither 32 nor 64
258 // bit indices present.
259 int LID(long long GID) const { return -1; }
260#endif
261
262#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
263 int* MyGlobalElements() const;
264 int IndexBase() const {
265 if(IndexBase64() == (long long) static_cast<int>(IndexBase64()))
266 return (int) IndexBase64();
267 throw "EpetraExt::LightweightMap::IndexBase: IndexBase cannot fit an int.";
268 }
269 void MyGlobalElementsPtr(int *& MyGlobalElementList) const;
270#endif
271
272#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
273 long long* MyGlobalElements64() const;
274 void MyGlobalElementsPtr(long long *& MyGlobalElementList) const;
275#endif
276 long long IndexBase64() const {return Data_->IndexBase_;}
277
278 int MinLID() const;
279 int MaxLID() const;
280
281 bool GlobalIndicesInt() const { return IsInt; }
282 bool GlobalIndicesLongLong() const { return IsLongLong; }
283 private:
284 void CleanupData();
285 LightweightMapData *Data_;
286 bool IsLongLong;
287 bool IsInt;
288 //Epetra_BlockMapData* Data_;
289
290#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
291 void Construct_int(int NumGlobalElements,int NumMyElements, const int * MyGlobalElements, long long IndexBase, bool GenerateHash=true);
292#endif
293#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
294 void Construct_LL(long long NumGlobalElements,int NumMyElements, const long long * MyGlobalElements, long long IndexBase, bool GenerateHash=true);
295#endif
296};
297
298
299// ==============================================================
301 public:
302 RemoteOnlyImport(const Epetra_Import & Importer, LightweightMap & RemoteOnlyTargetMap);
304
305 int NumSameIDs() {return 0;}
306
307 int NumPermuteIDs() {return 0;}
308
309 int NumRemoteIDs() {return NumRemoteIDs_;}
310
311 int NumExportIDs() {return NumExportIDs_;}
312
313 int* ExportLIDs() {return ExportLIDs_;}
314
315 int* ExportPIDs() {return ExportPIDs_;}
316
317 int* RemoteLIDs() {return RemoteLIDs_;}
318
319 int* PermuteToLIDs() {return 0;}
320
321 int* PermuteFromLIDs() {return 0;}
322
323 int NumSend() {return NumSend_;}
324
325 Epetra_Distributor & Distributor() {return *Distor_;}
326
327 const Epetra_BlockMap & SourceMap() const {return *SourceMap_;}
328 const LightweightMap & TargetMap() const {return *TargetMap_;}
329
330 private:
331 int NumSend_;
332 int NumRemoteIDs_;
333 int NumExportIDs_;
334 int * ExportLIDs_;
335 int * ExportPIDs_;
336 int * RemoteLIDs_;
337 Epetra_Distributor* Distor_;
338 const Epetra_BlockMap* SourceMap_;
339 const LightweightMap *TargetMap_;
340};
341
342// ==============================================================
344 public:
345 LightweightCrsMatrix(const Epetra_CrsMatrix & A, RemoteOnlyImport & RowImporter, bool SortGhosts=false,const char * label=0);
346 LightweightCrsMatrix(const Epetra_CrsMatrix & A, Epetra_Import & RowImporter);
348
349 // Standard crs data structures
350 std::vector<int> rowptr_;
351 std::vector<int> colind_;
352 std::vector<double> vals_;
353
354#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
355 // Colind in LL-GID space (if needed)
356 std::vector<long long> colind_LL_;
357#endif
358
359 // Epetra Maps
360 bool use_lw;
365
366
367 // List of owning PIDs (from the DomainMap) as ordered by entries in the column map.
368 std::vector<int> ColMapOwningPIDs_;
369
370 // For building the final importer for C
371 std::vector<int> ExportLIDs_;
372 std::vector<int> ExportPIDs_;
373
374 private:
375
376 template <typename ImportType, typename int_type>
377 void Construct(const Epetra_CrsMatrix & A, ImportType & RowImporter,bool SortGhosts=false, const char * label=0);
378
379 // Templated versions of MakeColMapAndReindex (to prevent code duplication)
380 template <class GO>
381 int MakeColMapAndReindex(std::vector<int> owningPIDs,std::vector<GO> Gcolind,bool SortGhosts=false, const char * label=0);
382
383 template<typename int_type>
384 std::vector<int_type>& getcolind();
385
386 template<typename ImportType, typename int_type>
387 int PackAndPrepareReverseComm(const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter,
388 std::vector<int> &ReverseSendSizes, std::vector<int_type> &ReverseSendBuffer);
389
390 template<typename ImportType, typename int_type>
391 int MakeExportLists(const Epetra_CrsMatrix & SourceMatrix, ImportType & RowImporter,
392 std::vector<int> &ReverseRecvSizes, const int_type *ReverseRecvBuffer,
393 std::vector<int> & ExportPIDs, std::vector<int> & ExportLIDs);
394
395};
396
397#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
398template<> inline std::vector<int>& LightweightCrsMatrix::getcolind() { return colind_; }
399#endif
400#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
401template<> inline std::vector<long long>& LightweightCrsMatrix::getcolind() { return colind_LL_; }
402#endif
403
404// ==============================================================
405template<typename int_type>
407 const Epetra_Map& targetMap,
408 CrsMatrixStruct& Mview,
409 const Epetra_Import * prototypeImporter=0,bool SortGhosts=false,
410 const char * label=0)
411{
412 // The goal of this method to populare the Mview object with ONLY the rows of M
413 // that correspond to elements in 'targetMap.' There will be no population of the
414 // numEntriesPerRow, indices, values, remote or numRemote arrays.
415
416
417 // The prototype importer, if used, has to know who owns all of the PIDs in M's rowmap.
418 // The typical use of this is to provide A's column map when I&XV is called for B, for
419 // a C = A * B multiply.
420
421#ifdef ENABLE_MMM_TIMINGS
422 std::string tpref;
423 if(label) tpref = std::string(label);
424 using Teuchos::TimeMonitor;
425 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string("EpetraExt: MMM Ionly Setup"))));
426#endif
427
428 Mview.deleteContents();
429
430 Mview.origMatrix = &M;
431 const Epetra_Map& Mrowmap = M.RowMap();
432 int numProcs = Mrowmap.Comm().NumProc();
433 Mview.numRows = targetMap.NumMyElements();
434
435 Mview.origRowMap = &(M.RowMap());
436 Mview.rowMap = &targetMap;
437 Mview.colMap = &(M.ColMap());
438 Mview.domainMap = &(M.DomainMap());
439 Mview.importColMap = NULL;
440
441 int i;
442 int numRemote =0;
443 int N = Mview.rowMap->NumMyElements();
444
445 if(Mrowmap.SameAs(targetMap)) {
446 numRemote = 0;
447 Mview.targetMapToOrigRow.resize(N);
448 for(i=0;i<N; i++) Mview.targetMapToOrigRow[i]=i;
449 return 0;
450 }
451 else if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)){
452 numRemote = prototypeImporter->NumRemoteIDs();
453
454 Mview.targetMapToOrigRow.resize(N); Mview.targetMapToOrigRow.assign(N,-1);
455 Mview.targetMapToImportRow.resize(N); Mview.targetMapToImportRow.assign(N,-1);
456
457 const int * PermuteToLIDs = prototypeImporter->PermuteToLIDs();
458 const int * PermuteFromLIDs = prototypeImporter->PermuteFromLIDs();
459 const int * RemoteLIDs = prototypeImporter->RemoteLIDs();
460
461 for(i=0; i<prototypeImporter->NumSameIDs();i++)
462 Mview.targetMapToOrigRow[i] = i;
463
464 // NOTE: These are reversed on purpose since we're doing a reverse map.
465 for(i=0; i<prototypeImporter->NumPermuteIDs();i++)
466 Mview.targetMapToOrigRow[PermuteToLIDs[i]] = PermuteFromLIDs[i];
467
468 for(i=0; i<prototypeImporter->NumRemoteIDs();i++)
469 Mview.targetMapToImportRow[RemoteLIDs[i]] = i;
470
471 }
472 else
473 throw std::runtime_error("import_only: This routine only works if you either have the right map or no prototypeImporter");
474
475 if (numProcs < 2) {
476 if (Mview.numRemote > 0) {
477 std::cerr << "EpetraExt::MatrixMatrix::Multiply ERROR, numProcs < 2 but "
478 << "attempting to import remote matrix rows."<<std::endl;
479 return(-1);
480 }
481
482 //If only one processor we don't need to import any remote rows, so return.
483 return(0);
484 }
485
486#ifdef ENABLE_MMM_TIMINGS
487 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string("EpetraExt: MMM Ionly Import-1"))));
488#endif
489 const int * RemoteLIDs = prototypeImporter->RemoteLIDs();
490
491 //Create a map that describes the remote rows of M that we need.
492 int_type* MremoteRows = numRemote>0 ? new int_type[prototypeImporter->NumRemoteIDs()] : 0;
493 for(i=0; i<prototypeImporter->NumRemoteIDs(); i++)
494 MremoteRows[i] = (int_type) targetMap.GID64(RemoteLIDs[i]);
495
496 LightweightMap MremoteRowMap((int_type) -1, numRemote, MremoteRows, (int_type)Mrowmap.IndexBase64());
497
498#ifdef ENABLE_MMM_TIMINGS
499 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string("EpetraExt: MMM Ionly Import-2"))));
500#endif
501 //Create an importer with target-map MremoteRowMap and
502 //source-map Mrowmap.
503 RemoteOnlyImport * Rimporter=0;
504 Rimporter = new RemoteOnlyImport(*prototypeImporter,MremoteRowMap);
505
506#ifdef ENABLE_MMM_TIMINGS
507 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string("EpetraExt: MMM Ionly Import-3"))));
508#endif
509
510 Mview.importMatrix = new LightweightCrsMatrix(M,*Rimporter,SortGhosts,label);
511
512#ifdef ENABLE_MMM_TIMINGS
513 MM = Teuchos::rcp(new TimeMonitor(*TimeMonitor::getNewTimer(tpref + std::string("EpetraExt: MMM Ionly Import-4"))));
514#endif
515
516#ifdef ENABLE_MMM_STATISTICS
517 printMultiplicationStatistics(Rimporter, label + std::string(" I&X MMM"));
518#endif
519
520 // Cleanup
521 delete Rimporter;
522 delete [] MremoteRows;
523
524 return(0);
525}
526
527
528
529// Statistics printing routines for when ENABLE_MMM_STATISTICS is enabled
530void PrintMultiplicationStatistics(Epetra_Import * Transfer, const std::string &label);
531void PrintMultiplicationStatistics(Epetra_Export * Transfer, const std::string &label);
532
533
534}//namespace EpetraExt
535
536#endif
537
const Epetra_CrsMatrix * origMatrix
std::vector< int > targetMapToImportRow
LightweightCrsMatrix * importMatrix
const Epetra_BlockMap * importColMap
std::vector< int > targetMapToOrigRow
CrsWrapper_Epetra_CrsMatrix(Epetra_CrsMatrix &epetracrsmatrix)
int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)
int InsertGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)
std::map< int_type, std::set< int_type > * > & get_graph()
int InsertGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
int SumIntoGlobalValues(int_type GlobalRow, int NumEntries, double *Values, int_type *Indices)
CrsWrapper_GraphBuilder(const Epetra_Map &emap)
virtual bool Filled()=0
virtual const Epetra_Map & RowMap() const =0
virtual int InsertGlobalValues(long long GlobalRow, int NumEntries, double *Values, long long *Indices)=0
virtual int SumIntoGlobalValues(long long GlobalRow, int NumEntries, double *Values, long long *Indices)=0
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
std::vector< long long > colind_LL_
LightweightCrsMatrix(const Epetra_CrsMatrix &A, RemoteOnlyImport &RowImporter, bool SortGhosts=false, const char *label=0)
std::vector< long long > MyGlobalElements_LL_
Epetra_HashTable< int > * LIDHash_int_
Epetra_HashTable< long long > * LIDHash_LL_
LightweightMap & operator=(const LightweightMap &map)
long long GID64(int LID) const
long long * MyGlobalElements64() const
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
const LightweightMap & TargetMap() const
Epetra_Distributor & Distributor()
RemoteOnlyImport(const Epetra_Import &Importer, LightweightMap &RemoteOnlyTargetMap)
const Epetra_BlockMap & SourceMap() const
const Epetra_Map & RowMap() const
const Epetra_Map & ColMap() const
const Epetra_Map & DomainMap() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
void insert_matrix_locations(CrsWrapper_GraphBuilder< int_type > &graphbuilder, Epetra_CrsMatrix &C)
void pack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int > &proc_col_ranges, std::vector< int > &send_rows, std::vector< int > &rows_per_send_proc)
void PrintMultiplicationStatistics(Epetra_Import *Transfer, const std::string &label)
void printMultiplicationStatistics(Epetra_Import *Transfer, const std::string &label)
int import_only(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter)
std::pair< int_type, int_type > get_col_range(const Epetra_Map &emap)
int dumpCrsMatrixStruct(const CrsMatrixStruct &M)
void Tpack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int_type > &proc_col_ranges, std::vector< int_type > &send_rows, std::vector< int > &rows_per_send_proc)