EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_BTF_CrsMatrix.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
42
44
45#include <Epetra_Import.h>
46#include <Epetra_CrsMatrix.h>
47#include <Epetra_CrsGraph.h>
48#include <Epetra_Map.h>
49
50#include <vector>
51
52using std::vector;
53using std::cout;
54using std::endl;
55
56#define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
57#define GENBTF_F77 F77_FUNC(genbtf,GENBTF)
58
59extern "C" {
60extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
61extern void GENBTF_F77( int*, int*, int*, int*, int*, int*, int*, int*, int*,
62 int*, int*, int*, int*, int*, int*, int*, int*, int*,
63 int*, int*, int* );
64}
65
66namespace EpetraExt {
67
70{
71 if( NewRowMap_ ) delete NewRowMap_;
72 if( NewColMap_ ) delete NewColMap_;
73
74 if( Importer_ ) delete Importer_;
75
76 if( NewMatrix_ ) delete NewMatrix_;
77 if( NewGraph_ ) delete NewGraph_;
78}
79
83{
84 origObj_ = &orig;
85
86 if( orig.RowMap().DistributedGlobal() )
87 { cout << "FAIL for Global!\n"; abort(); }
88 if( orig.IndicesAreGlobal() )
89 { cout << "FAIL for Global Indices!\n"; abort(); }
90
91 int n = orig.NumMyRows();
92 int nnz = orig.NumMyNonzeros();
93
94 if( verbose_ )
95 {
96 cout << "Orig Matrix:\n";
97 cout << orig << endl;
98 }
99
100 //create std CRS format
101 //also create graph without zero elements
102 vector<int> ia(n+1,0);
103 int maxEntries = orig.MaxNumEntries();
104 vector<int> ja_tmp(maxEntries);
105 vector<double> jva_tmp(maxEntries);
106 vector<int> ja(nnz);
107 int cnt;
108
109 const Epetra_BlockMap & OldRowMap = orig.RowMap();
110 const Epetra_BlockMap & OldColMap = orig.ColMap();
111 Epetra_CrsGraph strippedGraph( Copy, OldRowMap, OldColMap, 0 );
112
113 for( int i = 0; i < n; ++i )
114 {
115 orig.ExtractMyRowCopy( i, maxEntries, cnt, &jva_tmp[0], &ja_tmp[0] );
116 ia[i+1] = ia[i];
117 for( int j = 0; j < cnt; ++j )
118 if( fabs(jva_tmp[j]) > threshold_ )
119 ja[ ia[i+1]++ ] = ja_tmp[j];
120
121 int new_cnt = ia[i+1] - ia[i];
122 strippedGraph.InsertMyIndices( i, new_cnt, &ja[ ia[i] ] );
123 }
124 nnz = ia[n];
125 strippedGraph.FillComplete();
126
127 if( verbose_ )
128 {
129 cout << "Stripped Graph\n";
130 cout << strippedGraph;
131 }
132
133 vector<int> iat(n+1,0);
134 vector<int> jat(nnz);
135 for( int i = 0; i < n; ++i )
136 for( int j = ia[i]; j < ia[i+1]; ++j )
137 ++iat[ ja[j]+1 ];
138 for( int i = 0; i < n; ++i )
139 iat[i+1] += iat[i];
140 for( int i = 0; i < n; ++i )
141 for( int j = ia[i]; j < ia[i+1]; ++j )
142 jat[ iat[ ja[j] ]++ ] = i;
143 for( int i = 0; i < n; ++i )
144 iat[n-i] = iat[n-i-1];
145 iat[0] = 0;
146
147 //convert to Fortran indexing
148 for( int i = 0; i < n+1; ++i ) ++ia[i];
149 for( int i = 0; i < nnz; ++i ) ++ja[i];
150 for( int i = 0; i < n+1; ++i ) ++iat[i];
151 for( int i = 0; i < nnz; ++i ) ++jat[i];
152
153 if( verbose_ )
154 {
155 cout << "-----------------------------------------\n";
156 cout << "CRS Format Graph\n";
157 cout << "-----------------------------------------\n";
158 for( int i = 0; i < n; ++i )
159 {
160 cout << i+1 << ": " << ia[i+1] << ": ";
161 for( int j = ia[i]-1; j<ia[i+1]-1; ++j )
162 cout << " " << ja[j];
163 cout << endl;
164 }
165 cout << "-----------------------------------------\n";
166 }
167
168/*
169 vector<int> iat(n+1);
170 vector<int> jat(nnz);
171 int * jaf = &ja[0];
172 int * iaf = &ia[0];
173 int * jatf = &jat[0];
174 int * iatf = &iat[0];
175 MATTRANS_F77( &n, &n, &ja[0], &ia[0], &jat[0], &iat[0] );
176*/
177
178 if( verbose_ )
179 {
180 cout << "-----------------------------------------\n";
181 cout << "CCS Format Graph\n";
182 cout << "-----------------------------------------\n";
183 for( int i = 0; i < n; ++i )
184 {
185 cout << i+1 << ": " << iat[i+1] << ": ";
186 for( int j = iat[i]-1; j<iat[i+1]-1; ++j )
187 cout << " " << jat[j];
188 cout << endl;
189 }
190 cout << "-----------------------------------------\n";
191 }
192
193 vector<int> w(10*n);
194
195 vector<int> rowperm(n);
196 vector<int> colperm(n);
197
198 //horizontal block
199 int nhrows, nhcols, hrzcmp;
200 //square block
201 int nsrows, sqcmpn;
202 //vertial block
203 int nvrows, nvcols, vrtcmp;
204
205 vector<int> rcmstr(n+1);
206 vector<int> ccmstr(n+1);
207
208 int msglvl = 0;
209 int output = 6;
210
211 GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
212 &rowperm[0], &colperm[0], &nhrows, &nhcols,
213 &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
214 &rcmstr[0], &ccmstr[0], &msglvl, &output );
215
216 //convert back to C indexing
217 for( int i = 0; i < n; ++i )
218 {
219 --rowperm[i];
220 --colperm[i];
221 }
222 for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
223 {
224 --rcmstr[i];
225 --ccmstr[i];
226 }
227
228 if( verbose_ )
229 {
230 cout << "-----------------------------------------\n";
231 cout << "BTF Output\n";
232 cout << "-----------------------------------------\n";
233 cout << "RowPerm and ColPerm\n";
234 for( int i = 0; i<n; ++i )
235 cout << rowperm[i] << "\t" << colperm[i] << endl;
236 if( hrzcmp )
237 {
238 cout << "Num Horizontal: Rows, Cols, Comps\n";
239 cout << nhrows << "\t" << nhcols << "\t" << hrzcmp << endl;
240 }
241 cout << "Num Square: Rows, Comps\n";
242 cout << nsrows << "\t" << sqcmpn << endl;
243 if( vrtcmp )
244 {
245 cout << "Num Vertical: Rows, Cols, Comps\n";
246 cout << nvrows << "\t" << nvcols << "\t" << vrtcmp << endl;
247 }
248 cout << "Row, Col of upper left pt in blocks\n";
249 for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
250 cout << i << " " << rcmstr[i] << " " << ccmstr[i] << endl;
251 cout << "-----------------------------------------\n";
252 }
253
254 if( hrzcmp || vrtcmp )
255 {
256 cout << "FAILED! hrz cmp's:" << hrzcmp << " vrtcmp: " << vrtcmp << endl;
257 exit(0);
258 }
259
260 //convert rowperm to OLD->NEW
261 //reverse ordering of permutation to get upper triangular
262 vector<int> rowperm_t( n );
263 vector<int> colperm_t( n );
264 for( int i = 0; i < n; ++i )
265 {
266// rowperm_t[ rowperm[i] ] = i;
267 rowperm_t[i] = rowperm[i];
268 colperm_t[i] = colperm[i];
269 }
270
271 //Generate New Domain and Range Maps
272 //for now, assume they start out as identical
273 vector<int> myElements( n );
274 OldRowMap.MyGlobalElements( &myElements[0] );
275
276 vector<int> newDomainElements( n );
277 vector<int> newRangeElements( n );
278 for( int i = 0; i < n; ++i )
279 {
280 newDomainElements[ i ] = myElements[ rowperm_t[i] ];
281 newRangeElements[ i ] = myElements[ colperm_t[i] ];
282 }
283
284 NewRowMap_ = new Epetra_Map( n, n, &newDomainElements[0], OldRowMap.IndexBase(), OldRowMap.Comm() );
285 NewColMap_ = new Epetra_Map( n, n, &newRangeElements[0], OldColMap.IndexBase(), OldColMap.Comm() );
286
287 if( verbose_ )
288 {
289 cout << "New Row Map\n";
290 cout << *NewRowMap_ << endl;
291 cout << "New ColMap\n";
292 cout << *NewColMap_ << endl;
293 }
294
295 //Generate New Graph
296 NewGraph_ = new Epetra_CrsGraph( Copy, *NewRowMap_, *NewColMap_, 0 );
297 Importer_ = new Epetra_Import( *NewRowMap_, OldRowMap );
298 NewGraph_->Import( strippedGraph, *Importer_, Insert );
299 NewGraph_->FillComplete();
300 if( verbose_ )
301 {
302 cout << "NewGraph\n";
303 cout << *NewGraph_;
304 }
305
306 NewMatrix_ = new Epetra_CrsMatrix( Copy, *NewGraph_ );
307 NewMatrix_->Import( orig, *Importer_, Insert );
308 NewMatrix_->FillComplete();
309
310 if( verbose_ )
311 {
312 cout << "New CrsMatrix\n";
313 cout << *NewMatrix_ << endl;
314 }
315
316 newObj_ = NewMatrix_;
317
318 return *NewMatrix_;
319}
320
321bool
323fwd()
324{
325 int ret = NewMatrix_->Import( *origObj_, *Importer_, Insert );
326 if (ret<0) return false;
327 return true;
328}
329
330bool
332rvs()
333{
334 int ret = origObj_->Export( *NewMatrix_, *Importer_, Insert );
335 if (ret<0) return false;
336 return true;
337}
338
339} //namespace EpetraExt
#define GENBTF_F77
#define GENBTF_F77
#define MATTRANS_F77
Insert
Copy
NewTypeRef operator()(OriginalTypeRef orig)
Analysis of transform operation on original object and construction of new object.
bool rvs()
Reverse transfer of data from new object created in the operator() method call to the orig object inp...
bool fwd()
Forward transfer of data from orig object input in the operator() method call to the new object creat...
int MyGlobalElements(int *MyGlobalElementList) const
int IndexBase() const
const Epetra_Comm & Comm() const
int InsertMyIndices(int LocalRow, int NumIndices, int *Indices)
int FillComplete(bool OptimizeDataStorage=true)
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.