EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_AMD_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 
42 #include "EpetraExt_ConfigDefs.h"
43 #ifdef HAVE_EXPERIMENTAL
44 #ifdef HAVE_GRAPH_REORDERINGS
45 
46 #include <EpetraExt_AMD_CrsGraph.h>
47 
48 #include <Epetra_Import.h>
49 #include <Epetra_CrsGraph.h>
50 #include <Epetra_Map.h>
51 
52 #include <vector>
53 
54 extern "C" {
55 #include <amd.h>
56 }
57 
58 namespace EpetraExt {
59 
62 {
63  if( NewMap_ ) delete NewMap_;
64 
65  if( NewGraph_ ) delete NewGraph_;
66 }
67 
70 operator()( OriginalTypeRef orig )
71 {
72  origObj_ = &orig;
73 
74  int n = orig.NumMyRows();
75  int nnz = orig.NumMyNonzeros();
76 
77  //create std CRS format
78  std::vector<int> ia(n+1,0);
79  std::vector<int> ja(nnz);
80  int cnt;
81  for( int i = 0; i < n; ++i )
82  {
83  int * tmpP = &ja[ia[i]];
84  orig.ExtractMyRowCopy( i, nnz-ia[i], cnt, tmpP );
85  ia[i+1] = ia[i] + cnt;
86  }
87 
88  //trim down to local only
89  std::vector<int> iat(n+1);
90  std::vector<int> jat(nnz);
91  int loc = 0;
92  for( int i = 0; i < n; ++i )
93  {
94  iat[i] = loc;
95  for( int j = ia[i]; j < ia[i+1]; ++j )
96  {
97  if( ja[j] < n )
98  jat[loc++] = ja[j];
99  else
100  break;
101  }
102  }
103  iat[n] = loc;
104 
105 
106  if( verbose_ )
107  {
108  std::cout << "Orig Graph\n";
109  std::cout << orig << std::endl;
110  std::cout << "-----------------------------------------\n";
111  std::cout << "CRS Format Graph\n";
112  std::cout << "-----------------------------------------\n";
113  for( int i = 0; i < n; ++i )
114  {
115  std::cout << i << ": " << iat[i+1] << ": ";
116  for( int j = iat[i]; j<iat[i+1]; ++j )
117  std::cout << " " << jat[j];
118  std::cout << std::endl;
119  }
120  std::cout << "-----------------------------------------\n";
121  }
122 
123  std::vector<int> perm(n);
124  std::vector<double> info(AMD_INFO);
125 
126  amd_order( n, &iat[0], &jat[0], &perm[0], NULL, &info[0] );
127 
128  if( info[AMD_STATUS] == AMD_INVALID )
129  std::cout << "AMD ORDERING: Invalid!!!!\n";
130 
131  if( verbose_ )
132  {
133  std::cout << "-----------------------------------------\n";
134  std::cout << "AMD Output\n";
135  std::cout << "-----------------------------------------\n";
136  std::cout << "STATUS: " << info[AMD_STATUS] << std::endl;
137  std::cout << "SYMM: " << info[AMD_SYMMETRY] << std::endl;
138  std::cout << "N: " << info[AMD_N] << std::endl;
139  std::cout << "NZ: " << info[AMD_NZ] << std::endl;
140  std::cout << "SYMM: " << info[AMD_SYMMETRY] << std::endl;
141  std::cout << "NZDIAG: " << info[AMD_NZDIAG] << std::endl;
142  std::cout << "NZ A+At: " << info[AMD_NZ_A_PLUS_AT] << std::endl;
143  std::cout << "NDENSE: " << info[AMD_SYMMETRY] << std::endl;
144  std::cout << "Perm\n";
145  for( int i = 0; i<n; ++i )
146  std::cout << perm[i] << std::endl;
147  std::cout << "-----------------------------------------\n";
148  }
149 
150  //Generate New Domain and Range Maps
151  //for now, assume they start out as identical
152  const Epetra_BlockMap & OldMap = orig.RowMap();
153  int nG = orig.NumGlobalRows();
154 
155  std::vector<int> newElements( n );
156  for( int i = 0; i < n; ++i )
157  newElements[i] = OldMap.GID( perm[i] );
158 
159  NewMap_ = new Epetra_Map( nG, n, &newElements[0], OldMap.IndexBase(), OldMap.Comm() );
160 
161  if( verbose_ )
162  {
163  std::cout << "Old Map\n";
164  std::cout << OldMap << std::endl;
165  std::cout << "New Map\n";
166  std::cout << *NewMap_ << std::endl;
167  }
168 
169  //Generate New Graph
170  NewGraph_ = new Epetra_CrsGraph( Copy, *NewMap_, 0 );
171  Epetra_Import Importer( *NewMap_, OldMap );
172  NewGraph_->Import( orig, Importer, Insert );
173  NewGraph_->FillComplete();
174 
175  if( verbose_ )
176  {
177  std::cout << "New CrsGraph\n";
178  std::cout << *NewGraph_ << std::endl;
179  }
180 
181  newObj_ = NewGraph_;
182 
183  return *NewGraph_;
184 }
185 
186 } //namespace EpetraExt
187 
188 #endif //HAVE_GRAPH_REORDERINGS
189 #endif //HAVE_EXPERIMENTAL
~CrsGraph_AMD()
EpetraExt::CrsGraph_AMD Destructor.
int IndexBase() const
int GID(int LID) const
const Epetra_Comm & Comm() const
NewTypeRef operator()(OriginalTypeRef orig)
EpetraExt::CrsGraph_AMD Transform Operator.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)