EpetraExt Development
Loading...
Searching...
No Matches
EpetraExt_BTF_CrsGraph.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
43#include <Epetra_Import.h>
44#include <Epetra_CrsGraph.h>
45#include <Epetra_Map.h>
46
47#include <vector>
48
49using std::vector;
50using std::cout;
51using std::endl;
52
53#define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
54#define GENBTF_F77 F77_FUNC(genbtf,GENBTF)
55extern "C" {
56extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
57extern void GENBTF_F77( int*, int*, int*, int*, int*, int*, int*, int*, int*,
58 int*, int*, int*, int*, int*, int*, int*, int*, int*,
59 int*, int*, int* );
60}
61
62namespace EpetraExt {
63
66{
67 if( NewRowMap_ ) delete NewRowMap_;
68 if( NewDomainMap_ ) delete NewDomainMap_;
69
70 if( NewGraph_ ) delete NewGraph_;
71}
72
75operator()( OriginalTypeRef orig )
76{
77 origObj_ = &orig;
78
79 if( orig.RowMap().DistributedGlobal() )
80 { cout << "FAIL for Global!\n"; abort(); }
81 if( orig.IndicesAreGlobal() )
82 { cout << "FAIL for Global Indices!\n"; abort(); }
83
84 int n = orig.NumMyRows();
85 int nnz = orig.NumMyNonzeros();
86
87 //create std CRS format
88 vector<int> ia(n+1,0);
89 vector<int> ja(nnz);
90 int cnt;
91 for( int i = 0; i < n; ++i )
92 {
93 int * tmpP = &ja[ia[i]];
94 orig.ExtractMyRowCopy( i, nnz-ia[i], cnt, tmpP );
95 ia[i+1] = ia[i] + cnt;
96 }
97
98 //convert to Fortran indexing
99 for( int i = 0; i < n+1; ++i ) ++ia[i];
100 for( int i = 0; i < nnz; ++i ) ++ja[i];
101
102#ifdef BTF_VERBOSE
103 cout << "-----------------------------------------\n";
104 cout << "CRS Format Graph\n";
105 cout << "-----------------------------------------\n";
106 for( int i = 0; i < n; ++i )
107 {
108 cout << i << ": " << ia[i+1] << ": ";
109 for( int j = ia[i]-1; j<ia[i+1]-1; ++j )
110 cout << " " << ja[j];
111 cout << endl;
112 }
113 cout << "-----------------------------------------\n";
114#endif
115
116 vector<int> iat(n+1);
117 vector<int> jat(nnz);
118 int * jaf = &ja[0];
119 int * iaf = &ia[0];
120 int * jatf = &jat[0];
121 int * iatf = &iat[0];
122 MATTRANS_F77( &n, &n, jaf, iaf, jatf, iatf );
123
124#ifdef BTF_VERBOSE
125 cout << "-----------------------------------------\n";
126 cout << "CCS Format Graph\n";
127 cout << "-----------------------------------------\n";
128 for( int i = 0; i < n; ++i )
129 {
130 cout << i << ": " << iat[i+1] << ": ";
131 for( int j = iat[i]-1; j<iat[i+1]-1; ++j )
132 cout << " " << jat[j];
133 cout << endl;
134 }
135 cout << "-----------------------------------------\n";
136#endif
137
138 vector<int> w(10*n);
139
140 vector<int> rowperm(n);
141 vector<int> colperm(n);
142
143 //horizontal block
144 int nhrows, nhcols, hrzcmp;
145 //square block
146 int nsrows, sqcmpn;
147 //vertial block
148 int nvrows, nvcols, vrtcmp;
149
150 vector<int> rcmstr(n+1);
151 vector<int> ccmstr(n+1);
152
153 int msglvl = 0;
154 int output = 6;
155
156 GENBTF_F77( &n, &n, &iat[0], &jat[0], &ia[0], &ja[0], &w[0],
157 &rowperm[0], &colperm[0], &nhrows, &nhcols,
158 &hrzcmp, &nsrows, &sqcmpn, &nvrows, &nvcols, &vrtcmp,
159 &rcmstr[0], &ccmstr[0], &msglvl, &output );
160
161 //convert back to C indexing
162 for( int i = 0; i < n; ++i )
163 {
164 --rowperm[i];
165 --colperm[i];
166 }
167 for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
168 {
169 --rcmstr[i];
170 --ccmstr[i];
171 }
172
173#ifdef BTF_VERBOSE
174 cout << "-----------------------------------------\n";
175 cout << "BTF Output\n";
176 cout << "-----------------------------------------\n";
177 cout << "RowPerm and ColPerm\n";
178 for( int i = 0; i<n; ++i )
179 cout << rowperm[i] << "\t" << colperm[i] << endl;
180 if( hrzcmp )
181 {
182 cout << "Num Horizontal: Rows, Cols, Comps\n";
183 cout << nhrows << "\t" << nhcols << "\t" << hrzcmp << endl;
184 }
185 cout << "Num Square: Rows, Comps\n";
186 cout << nsrows << "\t" << sqcmpn << endl;
187 if( vrtcmp )
188 {
189 cout << "Num Vertical: Rows, Cols, Comps\n";
190 cout << nvrows << "\t" << nvcols << "\t" << vrtcmp << endl;
191 }
192 cout << "Row, Col of upper left pt in blocks\n";
193 for( int i = 0; (i<n+1) && (rcmstr[i]!=n+1); ++i )
194 cout << i << " " << rcmstr[i] << " " << ccmstr[i] << endl;
195 cout << "-----------------------------------------\n";
196#endif
197
198 if( hrzcmp || vrtcmp )
199 { cout << "FAILED! hrz cmp's:" << hrzcmp << " vrtcmp: " << vrtcmp << endl;
200 exit(0); }
201
202 //convert rowperm to OLD->NEW
203 //reverse ordering of permutation to get upper triangular
204 vector<int> rowperm_t( n );
205 vector<int> colperm_t( n );
206 for( int i = 0; i < n; ++i )
207 {
208// rowperm_t[ rowperm[i] ] = n-i;
209// colperm[i] = n-colperm[i];
210 rowperm_t[i] = rowperm[(n-1)-i];
211 colperm_t[i] = colperm[(n-1)-i];
212 }
213
214 //Generate New Domain and Range Maps
215 //for now, assume they start out as identical
216 const Epetra_BlockMap & OldMap = orig.RowMap();
217 vector<int> myElements( n );
218 OldMap.MyGlobalElements( &myElements[0] );
219
220 vector<int> newDomainElements( n );
221 vector<int> newRangeElements( n );
222 for( int i = 0; i < n; ++i )
223 {
224 newDomainElements[ i ] = myElements[ rowperm_t[i] ];
225 newRangeElements[ i ] = myElements[ colperm_t[i] ];
226cout << i << "\t" << rowperm_t[i] << "\t" << colperm[i] << "\t" << myElements[i] << endl;
227 }
228
229 NewRowMap_ = new Epetra_Map( n, n, &newDomainElements[0], OldMap.IndexBase(), OldMap.Comm() );
230 NewDomainMap_ = new Epetra_Map( n, n, &newRangeElements[0], OldMap.IndexBase(), OldMap.Comm() );
231
232#ifdef BTF_VERBOSE
233 cout << "New Row Map\n";
234 cout << *RowMap << endl;
235 cout << "New Domain Map\n";
236 cout << *DomainMap << endl;
237#endif
238
239 //Generate New Graph
240 NewGraph_ = new Epetra_CrsGraph( Copy, *NewRowMap_, 0 );
241 Epetra_Import Importer( *NewRowMap_, OldMap );
242 NewGraph_->Import( orig, Importer, Insert );
243 NewGraph_->FillComplete( *NewDomainMap_, *NewRowMap_ );
244
245#ifdef BTF_VERBOSE
246 cout << "New CrsGraph\n";
247 cout << *NewGraph_ << endl;
248#endif
249
250 newObj_ = NewGraph_;
251
252 return *NewGraph_;
253}
254
255} //namespace EpetraExt
#define GENBTF_F77
#define MATTRANS_F77
Insert
Copy
NewTypeRef operator()(OriginalTypeRef orig)
Construction of BTF ordered Epetra_CrsGraph from orig object.
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.