EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
41 #include <EpetraExt_BTF_CrsGraph.h>
42 
43 #include <Epetra_Import.h>
44 #include <Epetra_CrsGraph.h>
45 #include <Epetra_Map.h>
46 
47 #include <vector>
48 
49 using std::vector;
50 using std::cout;
51 using std::endl;
52 
53 #define MATTRANS_F77 F77_FUNC(mattrans,MATTRANS)
54 #define GENBTF_F77 F77_FUNC(genbtf,GENBTF)
55 extern "C" {
56 extern void MATTRANS_F77( int*, int*, int*, int*, int*, int* );
57 extern 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 
62 namespace EpetraExt {
63 
66 {
67  if( NewRowMap_ ) delete NewRowMap_;
68  if( NewDomainMap_ ) delete NewDomainMap_;
69 
70  if( NewGraph_ ) delete NewGraph_;
71 }
72 
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] ];
226 cout << 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
int MyGlobalElements(int *MyGlobalElementList) const
NewTypeRef operator()(OriginalTypeRef orig)
Construction of BTF ordered Epetra_CrsGraph from orig object.
int IndexBase() const
#define GENBTF_F77
const Epetra_Comm & Comm() const
#define MATTRANS_F77
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)