Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_SLU.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Amesos: Direct Sparse Solver Package
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// This library is free software; you can redistribute it and/or modify
11// it under the terms of the GNU Lesser General Public License as
12// published by the Free Software Foundation; either version 2.1 of the
13// License, or (at your option) any later version.
14//
15// This library is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18// Lesser General Public License for more details.
19//
20// You should have received a copy of the GNU Lesser General Public
21// License along with this library; if not, write to the Free Software
22// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23// USA
24// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25//
26// ***********************************************************************
27// @HEADER
28
29#include "Epetra_SLU.h"
30
31namespace SLU
32{
33 //FIXME
34 typedef int int_t;
35
36extern "C" {
37#include "dsp_defs.h"
38}
39}
40
41#include "Epetra_CrsGraph.h"
42#include "Epetra_LinearProblem.h"
43#include "Epetra_MultiVector.h"
44#include "Epetra_CrsMatrix.h"
45
46struct SLUData
47{
48 SLU::SuperMatrix A, B, X, L, U;
49 SLU::factor_param_t iparam;
50 SLU::mem_usage_t mem_usage;
51};
52
53Epetra_SLU::Epetra_SLU( Epetra_LinearProblem * Problem,
54 int fill_fac,
55 int panel_size,
56 int relax )
57: count_(0)
58{
59 using namespace SLU;
60
61 data_ = new SLUData();
62
63 data_->iparam.panel_size = panel_size;
64 data_->iparam.relax = relax;
65 data_->iparam.diag_pivot_thresh = -1;
66 data_->iparam.drop_tol = -1;
67
68 A_ = dynamic_cast<Epetra_CrsMatrix *> (Problem->GetOperator());
69 X_ = Problem->GetLHS();
70 B_ = Problem->GetRHS();
71
72 AG_ = new Epetra_CrsGraph( A_->Graph() );
73
74 int NumIndices;
75 int NumMyCols = AG_->NumMyCols();
76 int NumMyEqs = A_->NumMyRows();
77 int GlobalMaxNumIndices = AG_->GlobalMaxNumIndices();
78 int * xIndices;
79
80 TransNumNZ_ = new int[NumMyCols];
81 for( int i = 0; i < NumMyCols; i++ )
82 TransNumNZ_[i] = 0;
83 for( int i = 0; i < NumMyEqs; i++ )
84 {
85 assert( AG_->ExtractMyRowView( i, NumIndices, xIndices ) == 0 );
86 for( int j = 0; j < NumIndices; j++ )
87 ++TransNumNZ_[ xIndices[j] ];
88 }
89 TotNumNZ_ = 0;
90 for( int i = 0; i < NumMyCols; i++ )
92
93 TransIndices_ = new int*[ NumMyCols ];
94 TransValues_ = new double*[ NumMyCols ];
95
96 RowIndices_ = new int[TotNumNZ_];
97 RowValues_ = new double[TotNumNZ_];
98 ColPointers_ = new int[NumMyCols+1];
99
100 // Set pointers into the RowIndices and Values arrays, define ColPointers
101 NumIndices = 0;
102 for(int i=0; i<NumMyCols; i++) {
103 ColPointers_[i] = NumIndices;
104 TransIndices_[i] = RowIndices_ + NumIndices;
105 TransValues_[i] = RowValues_ + NumIndices;
106 NumIndices += TransNumNZ_[i];
107 }
108 ColPointers_[NumMyCols] = NumIndices;
109
110 Copy();
111
112 int m = A_->NumGlobalRows();
113 int n = A_->NumGlobalCols();
114
115 /* Create matrix A in the format expected by SuperLU. */
116 dCreate_CompCol_Matrix( &(data_->A), m, n, TotNumNZ_, RowValues_,
117 RowIndices_, ColPointers_, SLU_NC, SLU_D, SLU_GE );
118
119 /* Create right-hand side matrix B. */
120 double * rhs_x;
121 double * rhs_b;
122 int LDA_x, LDA_b;
123 X_->ExtractView( &rhs_x, &LDA_x );
124 dCreate_Dense_Matrix( &(data_->X), m, 1, rhs_x, m, SLU_DN, SLU_D, SLU_GE);
125 B_->ExtractView( &rhs_b, &LDA_b );
126 dCreate_Dense_Matrix( &(data_->B), m, 1, rhs_b, m, SLU_DN, SLU_D, SLU_GE);
127
128 perm_r_ = new int[m];
129 perm_c_ = new int[n];
130
131 etree_ = new int[n];
132
133 ferr_ = new double[1];
134 berr_ = new double[1];
135
136 R_ = new double[m];
137 C_ = new double[n];
138}
139
141{
142 SLU::Destroy_SuperMatrix_Store( &(data_->A) );
143 SLU::Destroy_SuperMatrix_Store( &(data_->B) );
144 SLU::Destroy_SuperMatrix_Store( &(data_->X) );
145
146 delete [] TransIndices_;
147 delete [] TransValues_;
148
149 delete [] TransNumNZ_;
150
151 delete [] RowIndices_;
152 delete [] RowValues_;
153 delete [] ColPointers_;
154
155 delete [] perm_r_;
156 delete [] perm_c_;
157
158 delete [] etree_;
159
160 delete [] ferr_;
161 delete [] berr_;
162
163 delete [] R_;
164 delete [] C_;
165
166// Destroy_CompCol_Matrix(data_->A);
167// SLU::Destroy_SuperNode_Matrix( &(data_->L) );
168// SLU::Destroy_CompCol_Matrix( &(data_->U) );
169
170 delete data_;
171
172 delete AG_;
173}
174
176{
177 int NumMyCols = AG_->NumMyCols();
178 int NumMyEqs = A_->NumMyRows();
179 int GlobalMaxNumIndices = AG_->GlobalMaxNumIndices();
180 int NumIndices;
181 int * xIndices;
182 double * xValues;
183
184 for ( int i = 0; i < NumMyCols; i++ )
185 TransNumNZ_[i] = 0;
186
187 for ( int i = 0; i < NumMyEqs; i++ )
188 {
189 assert(A_->ExtractMyRowView( i, NumIndices, xValues, xIndices) == 0 );
190 int ii = A_->GRID(i);
191 for ( int j = 0; j < NumIndices; j++ )
192 {
193 int TransRow = xIndices[j];
194 int loc = TransNumNZ_[TransRow];
195 TransIndices_[TransRow][loc] = ii;
196 TransValues_[TransRow][loc] = xValues[j];
197 ++TransNumNZ_[TransRow]; // increment counter into current transpose row
198 }
199 }
200}
201
202int Epetra_SLU::Solve( bool Verbose,
203 bool Equil,
204 bool Factor,
205 int perm_type,
206 double pivot_thresh,
207 bool Refact,
208 bool Trans )
209{
210 Copy();
211
212 using namespace SLU;
213
214 int m = A_->NumGlobalRows();
215
216 int permt = perm_type;
217 if( m < 3 ) permt = 0;
218
219 /*
220 * Get column permutation vector perm_c[], according to permc_spec:
221 * permc_spec = 0: use the natural ordering
222 * permc_spec = 1: use minimum degree ordering on structure of A'*A
223 * permc_spec = 2: use minimum degree ordering on structure of A'+A
224 */
225 if( !count_ ) get_perm_c( permt, &(data_->A), perm_c_ );
226
227 if( Verbose ) cout << "MATRIX COPIED!" << endl;
228
229 int info = 0;
230
231 char fact, trans, refact, equed;
232
233 if( Trans ) trans = 'T';
234 else trans = 'N';
235
236 if( Equil ) fact = 'E';
237 else fact = 'N';
238
239 if( !count_ || !Refact )
240 {
241 refact = 'N';
242 count_ = 0;
243 }
244 else
245 {
246 refact = 'Y';
247 if( !Factor ) fact = 'F';
248 }
249
250 if( Equil ) equed = 'B';
251 else equed = 'N';
252
253 data_->iparam.diag_pivot_thresh = pivot_thresh;
254
255 double rpg, rcond;
256
257 if( Verbose ) cout << "TRANS: " << trans << endl;
258 if( Verbose ) cout << "REFACT: " << refact << endl;
259
260 dgssvx( &fact, &trans, &refact, &(data_->A), &(data_->iparam), perm_c_,
261 perm_r_, etree_, &equed, R_, C_, &(data_->L), &(data_->U),
262 NULL, 0, &(data_->B), &(data_->X), &rpg, &rcond,
263 ferr_, berr_, &(data_->mem_usage), &info );
264
265 if( info ) cout << "WARNING: SuperLU returned with error code = " << info << endl;
266
267 if( Verbose )
268 {
269 cout << "SYSTEM DIRECT SOLVED!" << endl;
270
271 cout << "SuperLU INFO: " << info << "\n\n";
272 cout << " ncols = " << A_->NumGlobalCols() << endl;
273
274 cout << "SuperLU Memory Usage\n";
275 cout << "--------------------\n";
276 cout << "LU: " << data_->mem_usage.for_lu << endl;
277 cout << "Total: " << data_->mem_usage.total_needed << endl;
278 cout << "Expands: " << data_->mem_usage.expansions << endl;
279 cout << "--------------------\n\n";
280
281 if (m<200) dPrint_CompCol_Matrix("A", &(data_->A));
282// if (m<200) dPrint_CompCol_Matrix("U", &(data_->U));
283// if (m<200) dPrint_SuperNode_Matrix("L", &(data_->L));
284// if (m<200) PrintInt10("\nperm_r", m, perm_r);
285 if (m<200) dPrint_Dense_Matrix("B", &(data_->B));
286 if (m<200) dPrint_Dense_Matrix("X", &(data_->X));
287 }
288
289 count_++;
290
291 return 0;
292}
293
int * ColPointers_
Definition Epetra_SLU.h:114
int * TransNumNZ_
Definition Epetra_SLU.h:106
double * RowValues_
Definition Epetra_SLU.h:113
Epetra_MultiVector * B_
Definition Epetra_SLU.h:94
double * berr_
Definition Epetra_SLU.h:124
Epetra_CrsGraph * AG_
Definition Epetra_SLU.h:99
int * perm_r_
Definition Epetra_SLU.h:116
int * etree_
Definition Epetra_SLU.h:121
double * R_
Definition Epetra_SLU.h:96
int * RowIndices_
Definition Epetra_SLU.h:112
Epetra_CrsMatrix * A_
Definition Epetra_SLU.h:92
double * C_
Definition Epetra_SLU.h:97
double * ferr_
Definition Epetra_SLU.h:123
Epetra_MultiVector * X_
Definition Epetra_SLU.h:93
int Solve(bool Verbose=false, bool Equil=true, bool Factor=true, int perm_type=2, double pivot_thresh=-1, bool Refact=true, bool Trans=false)
All computation is performed during the call to Solve()
Epetra_SLU(Epetra_LinearProblem *Problem, int fill_fac=-1, int panel_size=-1, int relax=-1)
Epetra_SLU Constructor.
SLUData * data_
Definition Epetra_SLU.h:101
~Epetra_SLU()
Epetra_SLU Destructor.
double ** TransValues_
Definition Epetra_SLU.h:110
int ** TransIndices_
Definition Epetra_SLU.h:109
int * perm_c_
Definition Epetra_SLU.h:117
SLU::SuperMatrix A
SLU::mem_usage_t mem_usage
SLU::SuperMatrix B
SLU::SuperMatrix L
SLU::factor_param_t iparam
SLU::SuperMatrix U
SLU::SuperMatrix X
int int_t