EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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
43#ifdef HAVE_EXPERIMENTAL
44#ifdef HAVE_GRAPH_REORDERINGS
45
47
48#include <Epetra_Import.h>
49#include <Epetra_CrsGraph.h>
50#include <Epetra_Map.h>
51
52#include <vector>
53
54extern "C" {
55#include <amd.h>
56}
57
58namespace EpetraExt {
59
62{
63 if( NewMap_ ) delete NewMap_;
64
65 if( NewGraph_ ) delete NewGraph_;
66}
67
70operator()( 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 );
174
175 if( verbose_ )
176 {
177 std::cout << "New CrsGraph\n";
178 std::cout << *NewGraph_ << std::endl;
179 }
180
182
183 return *NewGraph_;
184}
185
186} //namespace EpetraExt
187
188#endif //HAVE_GRAPH_REORDERINGS
189#endif //HAVE_EXPERIMENTAL
NewTypeRef operator()(OriginalTypeRef orig)
EpetraExt::CrsGraph_AMD Transform Operator.
~CrsGraph_AMD()
EpetraExt::CrsGraph_AMD Destructor.
int GID(int LID) const
int IndexBase() const
const Epetra_Comm & Comm() const
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.