Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_Taucs.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 "Amesos_Taucs.h"
30#include "Epetra_Map.h"
31#include "Epetra_Import.h"
32#include "Epetra_CrsMatrix.h"
33#include "Epetra_Vector.h"
34#include "Epetra_Util.h"
35extern "C" {
36#include "taucs.h"
37}
38
39using namespace Teuchos;
40
41
42//
43// deaalocFunctorHandleDelete requires a fucntion that takes a ** argument, these
44// functions meet that requirement
45//
46void taucs_dccs_free_ptr( taucs_ccs_matrix** taucs_ccs_matrix_in ) {
47 taucs_dccs_free( *taucs_ccs_matrix_in );
48}
49
50void taucs_supernodal_factor_free_ptr( taucs_ccs_matrix** taucs_ccs_matrix_in ) {
51 taucs_supernodal_factor_free( *taucs_ccs_matrix_in );
52}
53
55public:
56
57#define USE_REF_COUNT_PTR_FOR_A_AND_L
58#ifdef USE_REF_COUNT_PTR_FOR_A_AND_L
59 Teuchos::RCP<taucs_ccs_matrix> A_ ;
60 Teuchos::RCP<taucs_ccs_matrix> L_ ;
61#else
62 taucs_ccs_matrix* A_ ;
63 taucs_ccs_matrix* L_ ;
64
66 A_(0),
67 L_(0)
68 {}
69
71 if (A_ != 0)
72 {
73 taucs_dccs_free(A_);
74 A_ = 0;
75 }
76
77 if (L_ != 0)
78 {
79 taucs_supernodal_factor_free(L_);
80 L_ = 0;
81 }
82 }
83
84#endif
85} ;
86
87//=============================================================================
88Amesos_Taucs::Amesos_Taucs(const Epetra_LinearProblem &prob) :
89 PrivateTaucsData_( rcp( new Amesos_Taucs_Pimpl() ) ),
90 Matrix_(0),
91 Problem_(&prob),
92 MtxConvTime_(-1),
93 MtxRedistTime_(-1),
94 VecRedistTime_(-1),
95 SymFactTime_(-1),
96 NumFactTime_(-1),
97 SolveTime_(-1)
98{ }
99
100//=============================================================================
102{
103
104 // print out some information if required by the user
105 if ((verbose_ && PrintTiming_) || verbose_ == 2) PrintTiming();
106 if ((verbose_ && PrintStatus_) || verbose_ == 2) PrintStatus();
107
108#if 0
109 if (A_ != 0)
110 {
111 taucs_dccs_free(A_);
112 A_ = 0;
113 }
114
115 if (L_ != 0)
116 {
117 taucs_supernodal_factor_free(L_);
118 L_ = 0;
119 }
120#endif
121}
122
123//=============================================================================
125{
126 ResetTimer();
127
128 int NumGlobalRows = Matrix_->NumGlobalRows();
129
130 // create a serial map
131 int NumMyRows = 0;
132 if (Comm().MyPID() == 0)
133 NumMyRows = NumGlobalRows;
134
135 SerialMap_ = rcp(new Epetra_Map(-1, NumMyRows, 0, Comm()));
136 if (SerialMap_.get() == 0)
137 AMESOS_CHK_ERR(-1);
138
139 Importer_ = rcp(new Epetra_Import(SerialMap(),Map()));
140 if (Importer_.get() == 0)
141 AMESOS_CHK_ERR(-1);
142
143 SerialCrsMatrix_ = rcp(new Epetra_CrsMatrix(Copy, SerialMap(), 0));
144 if (SerialCrsMatrix_.get() == 0)
145 AMESOS_CHK_ERR(-1);
146
147 AMESOS_CHK_ERR(SerialCrsMatrix().Import(Matrix(), Importer(), Add));
148
149 AMESOS_CHK_ERR(SerialCrsMatrix().FillComplete());
150
151 SerialMatrix_ = rcp(SerialCrsMatrix_.get(), false);
152
153 MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_);
154
155 return 0;
156}
157
158//=============================================================================
160{
161 ResetTimer();
162
163 if (Comm().MyPID() == 0)
164 {
165 int n = SerialMatrix().NumMyRows();
166 int nnz = SerialMatrix().NumMyNonzeros();
167 int nnz_sym = (nnz) / 2 + n;
168
169#ifndef USE_REF_COUNT_PTR_FOR_A_AND_L
170 if (PrivateTaucsData_->A_ != 0)
171 {
172 taucs_dccs_free(PrivateTaucsData_->A_);
173 PrivateTaucsData_->A_ = 0;
174 }
175#endif
176
177#ifdef USE_REF_COUNT_PTR_FOR_A_AND_L
178 PrivateTaucsData_->A_ =
179 rcp( taucs_ccs_create(n, n, nnz_sym, TAUCS_DOUBLE),
180 deallocFunctorHandleDelete<taucs_ccs_matrix>(taucs_dccs_free_ptr), true );
181#else
182 PrivateTaucsData_->A_ = taucs_ccs_create(n, n, nnz_sym, TAUCS_DOUBLE) ;
183#endif
184 PrivateTaucsData_->A_->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
185
186 int count = 0;
187 int MaxNumEntries = SerialMatrix().MaxNumEntries();
188 std::vector<int> Indices(MaxNumEntries);
189 std::vector<double> Values(MaxNumEntries);
190
191 PrivateTaucsData_->A_->colptr[0] = 0;
192
193 for (int i = 0 ; i < n ; ++i)
194 {
195 int count2 = 0;
196 int ierr, NumEntriesThisRow;
197 ierr = SerialMatrix().ExtractMyRowCopy(i, MaxNumEntries,
198 NumEntriesThisRow,
199 &Values[0], &Indices[0]);
200 if (ierr < 0)
201 AMESOS_CHK_ERR(ierr);
202
203 for (int j = 0 ; j < NumEntriesThisRow ; ++j)
204 {
205 if (Indices[j] < i)
206 continue;
207
208 if (Indices[j] == i)
209 Values[j] += AddToDiag_;
210
211 PrivateTaucsData_->A_->rowind[count] = Indices[j];
212 PrivateTaucsData_->A_->values.d[count] = Values[j];
213 ++count;
214 ++count2;
215 }
216
217 PrivateTaucsData_->A_->colptr[i + 1] = PrivateTaucsData_->A_->colptr[i] + count2;
218 }
219
220 if (count > SerialMatrix().NumMyNonzeros())
221 AMESOS_CHK_ERR(-1); // something wrong here
222 }
223
224 MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
225
226 return 0;
227}
228
229//=============================================================================
230int Amesos_Taucs::SetParameters( Teuchos::ParameterList &ParameterList)
231{
232 // retrive general parameters
233
234 SetStatusParameters( ParameterList );
235
236 SetControlParameters( ParameterList );
237
238 return 0;
239}
240
241//=============================================================================
243{
244 ResetTimer();
245
246 if (Comm().MyPID() == 0)
247 {
248
249#ifdef USE_REF_COUNT_PTR_FOR_A_AND_L
250 PrivateTaucsData_->L_ =
251 rcp( ((taucs_ccs_matrix*)taucs_ccs_factor_llt_symbolic(&*PrivateTaucsData_->A_)),
252 deallocFunctorHandleDelete<taucs_ccs_matrix>(taucs_supernodal_factor_free_ptr), true );
253#else
254 if (PrivateTaucsData_->L_ != 0)
255 {
256 taucs_supernodal_factor_free(PrivateTaucsData_->L_);
257 PrivateTaucsData_->L_ = 0;
258 }
259 PrivateTaucsData_->L_ =
260 (taucs_ccs_matrix*)taucs_ccs_factor_llt_symbolic(&*PrivateTaucsData_->A_) ;
261#endif
262 }
263
264 SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
265
266 return 0;
267}
268
269//=============================================================================
271{
272 ResetTimer();
273
274 if (Comm().MyPID() == 0)
275 {
276 taucs_supernodal_factor_free_numeric( &*PrivateTaucsData_->L_);
277
278 int ierr = taucs_ccs_factor_llt_numeric(&*PrivateTaucsData_->A_, &*PrivateTaucsData_->L_);
279
280 if (ierr != 0)
281 {
282 std::cerr << "Amesos_Taucs: error during numeric factorization ("
283 << ierr << ")" << std::endl;
284 AMESOS_CHK_ERR(-1);
285 }
286 }
287
288 NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
289
290 return 0;
291}
292
293//=============================================================================
295{
296 bool OK = true;
297
298 if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
299 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() )
300 {
301 OK = false;
302 }
303 return OK;
304}
305
306//=============================================================================
308{
309 if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__ << " Entering SymbolicFactorization()" << std::endl ;
312
313 CreateTimer(Comm());
314
316
317 Matrix_ = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
318 Map_ = &(Matrix_->RowMatrixRowMap());
319
320 // =========================================================== //
321 // redistribute and create all import/export objects only //
322 // if more than one processor is used. Otherwise simply set //
323 // dummy pointers to Matrix() and Map(), without giving the //
324 // ownership to the smart pointer. //
325 // =========================================================== //
326
327 if (Comm().NumProc() != 1)
329 else
330 {
331 SerialMap_ = rcp(const_cast<Epetra_Map*>(&Map()), false);
332 SerialMatrix_ = rcp(const_cast<Epetra_RowMatrix*>(&Matrix()), false);
333 }
334
335 // =========================================================== //
336 // Only on processor zero, convert the matrix into CSR format, //
337 // as required by TAUCS. //
338 // =========================================================== //
339
341
343
345
346 if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__ << " Leaving SymbolicFactorization()" << std::endl ;
347 return(0);
348}
349
350//=============================================================================
352{
353 if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__ << " Entering NumericFactorization()" << std::endl ;
355
356 if (IsSymbolicFactorizationOK_ == false)
358
360
361 // FIXME: this must be checked, now all the matrix is shipped twice here
364
366
368
369 if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__
370 << " Leaving NumericFactorization()" << std::endl ;
371 return(0);
372}
373
374//=============================================================================
376{
377 if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__
378 << " Entering Solve()" << std::endl ;
379 if (IsNumericFactorizationOK_ == false)
381
382 ++NumSolve_;
383
384 Epetra_MultiVector* X = Problem_->GetLHS();
385 Epetra_MultiVector* B = Problem_->GetRHS();
386
387 if ((X == 0) || (B == 0))
388 AMESOS_CHK_ERR(-1);
389
390 int NumVectors = X->NumVectors();
391 if (NumVectors != B->NumVectors())
392 AMESOS_CHK_ERR(-1);
393
394 // vectors with SerialMap_
395 RCP<Epetra_MultiVector> SerialB;
396 RCP<Epetra_MultiVector> SerialX;
397
398 ResetTimer();
399
400 if (Comm().NumProc() == 1)
401 {
402 SerialB = rcp(B,false);
403 SerialX = rcp(X,false);
404 }
405 else
406 {
407 SerialX = rcp(new Epetra_MultiVector(SerialMap(),NumVectors));
408 SerialB = rcp(new Epetra_MultiVector(SerialMap(),NumVectors));
409
410 SerialB->Import(*B,Importer(),Insert);
411 }
412
413 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
414
415 ResetTimer();
416
417 if (Comm().MyPID() == 0)
418 {
419 double* SerialXValues;
420 double* SerialBValues;
421 int LDA;
422
423 AMESOS_CHK_ERR(SerialX->ExtractView(&SerialXValues,&LDA));
424
425 // FIXME: check LDA
426 AMESOS_CHK_ERR(SerialB->ExtractView(&SerialBValues,&LDA));
427
428 for (int i = 0 ; i < NumVectors ; ++i)
429 {
430 int ierr = taucs_supernodal_solve_llt(&*PrivateTaucsData_->L_,
431 SerialXValues + i * LDA,
432 SerialBValues + i * LDA);
433
434 if (ierr != TAUCS_SUCCESS)
435 {
436 std::cerr << "Error occurred in taucs_ccs_solve()" << std::endl;
437 AMESOS_CHK_ERR(-1);
438 }
439 }
440 }
441
442 SolveTime_ = AddTime("Total solve time", SolveTime_);
443
444 ResetTimer();
445
446 if (Comm().NumProc() != 1)
447 X->Export(*SerialX, Importer(), Insert);
448
449 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
450
452 ComputeTrueResidual(Matrix(), *X, *B, false, "Amesos_Taucs");
453
455 ComputeVectorNorms(*X, *B, "Amesos_Taucs");
456
457 if ( debug_ > 0 ) std::cout << __FILE__ << "::" << __LINE__
458 << " Leaving Solve()" << std::endl ;
459 return(0) ;
460}
461
462// ======================================================================
464{
465 if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
466 return;
467
468 std::string p = "Amesos_Taucs : ";
469
470 if (Matrix_ == 0) return;
471
472
473 PrintLine();
474 int n = Matrix().NumGlobalRows();
475 int nnz = Matrix().NumGlobalNonzeros();
476
477 std::cout << p << "Matrix has " << n << " rows"
478 << " and " << nnz << " nonzeros" << std::endl;
479 if (n > 0)
480 {
481 std::cout << p << "Nonzero elements per row = "
482 << 1.0 * nnz / n << std::endl;
483 std::cout << p << "Percentage of nonzero elements = "
484 << 100.0 * nnz /(pow(n,2.0)) << std::endl;
485 }
486 PrintLine();
487
488 return;
489}
490
491// ======================================================================
493{
494 if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
495 return;
496
497 double ConTime = GetTime(MtxConvTime_);
498 double MatTime = GetTime(MtxRedistTime_);
499 double VecTime = GetTime(VecRedistTime_);
500 double SymTime = GetTime(SymFactTime_);
501 double NumTime = GetTime(NumFactTime_);
502 double SolTime = GetTime(SolveTime_);
503
505 SymTime /= NumSymbolicFact_;
506
507 if (NumNumericFact_)
508 NumTime /= NumNumericFact_;
509
510 if (NumSolve_)
511 SolTime /= NumSolve_;
512
513 std::string p = "Amesos_Taucs : ";
514 PrintLine();
515
516 std::cout << p << "Time to convert matrix to Taucs format = "
517 << ConTime << " (s)" << std::endl;
518 std::cout << p << "Time to redistribute matrix = "
519 << MatTime << " (s)" << std::endl;
520 std::cout << p << "Time to redistribute vectors = "
521 << VecTime << " (s)" << std::endl;
522 std::cout << p << "Number of symbolic factorizations = "
523 << NumSymbolicFact_ << std::endl;
524 std::cout << p << "Time for sym fact = "
525 << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
526 std::cout << p << "Number of numeric factorizations = "
527 << NumNumericFact_ << std::endl;
528 std::cout << p << "Time for num fact = "
529 << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
530 std::cout << p << "Number of solve phases = "
531 << NumSolve_ << std::endl;
532 std::cout << p << "Time for solve = "
533 << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
534
535 PrintLine();
536
537 return;
538}
#define AMESOS_CHK_ERR(a)
void taucs_supernodal_factor_free_ptr(taucs_ccs_matrix **taucs_ccs_matrix_in)
void taucs_dccs_free_ptr(taucs_ccs_matrix **taucs_ccs_matrix_in)
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
double AddToDiag_
Add this value to the diagonal.
int debug_
Sets the level of debug_ output.
bool PrintTiming_
If true, prints timing information in the destructor.
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
int verbose_
Toggles the output level.
int NumSymbolicFact_
Number of symbolic factorization phases.
bool IsNumericFactorizationOK_
If true, NumericFactorization() has been successfully called.
int NumSolve_
Number of solves.
bool ComputeTrueResidual_
If true, computes the true residual in Solve().
bool PrintStatus_
If true, print additional information in the destructor.
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
bool IsSymbolicFactorizationOK_
If true, SymbolicFactorization() has been successfully called.
int NumNumericFact_
Number of numeric factorization phases.
Teuchos::RCP< taucs_ccs_matrix > A_
Teuchos::RCP< taucs_ccs_matrix > L_
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
int ConvertToTaucs()
Converts the Epetra_RowMatrix into TAUCS format.
const Epetra_RowMatrix & Matrix() const
Returns a reference to the linear system matrix.
int NumericFactorization()
Performs NumericFactorization on the matrix A.
Epetra_CrsMatrix & SerialCrsMatrix()
Returns a reference to the already SerialMatrix as Crs (if allocated).
int ConvertToSerial()
Constructs a matrix with all rows on processor 0.
Teuchos::RCP< Amesos_Taucs_Pimpl > PrivateTaucsData_
Amesos_Taucs(const Epetra_LinearProblem &LinearProblem)
Default constructor.
bool MatrixShapeOK() const
Returns true if the solver can handle this matrix shape.
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Epetra_Import & Importer()
Returns a reference to the already allocated Importer.
Teuchos::RCP< Epetra_Import > Importer_
Epetra_Map & SerialMap()
Returns a reference to the already allocated SerialMap.
Teuchos::RCP< Epetra_RowMatrix > SerialMatrix_
int MtxConvTime_
Quick accessor pointer to internal timing data.
Epetra_RowMatrix & SerialMatrix()
Returns a reference to the SerialMatrix.
int PerformSymbolicFactorization()
Performs the symbolic factorization.
int PerformNumericFactorization()
Performs the numeric factorization.
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
const Epetra_LinearProblem * GetProblem() const
Returns the Epetra_LinearProblem.
Teuchos::RCP< Epetra_CrsMatrix > SerialCrsMatrix_
const Epetra_RowMatrix * Matrix_
Teuchos::RCP< Epetra_Map > SerialMap_
void PrintTiming() const
Prints timing information.
const Epetra_LinearProblem * Problem_
Pointer to the linear system problem.
void PrintStatus() const
Prints status information.
~Amesos_Taucs(void)
Default destructor.
int Solve()
Solves A X = B (or AT x = B)
const Epetra_Map & Map() const
Returns a reference to the RowMatrixRowMap().
const Epetra_Map * Map_
double GetTime(const std::string what) const
Gets the cumulative time using the string.
void ResetTimer(const int timerID=0)
Resets the internally stored time object.
Definition Amesos_Time.h:74
int AddTime(const std::string what, int dataID, const int timerID=0)
Adds to field what the time elapsed since last call to ResetTimer().
Definition Amesos_Time.h:80
void CreateTimer(const Epetra_Comm &Comm, int size=1)
Initializes the Time object.
Definition Amesos_Time.h:64
void ComputeTrueResidual(const Epetra_RowMatrix &Matrix, const Epetra_MultiVector &X, const Epetra_MultiVector &B, const bool UseTranspose, const std::string prefix) const
Computes the true residual, B - Matrix * X, and prints the results.
void ComputeVectorNorms(const Epetra_MultiVector &X, const Epetra_MultiVector &B, const std::string prefix) const
Computes the norms of X and B and print the results.
void PrintLine() const
Prints line on std::cout.