Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Amesos_CSparse.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#ifdef HAVE_AMESOS_CSPARSE
30#include "Amesos_CSparse.h"
31#include "Epetra_Map.h"
32#include "Epetra_Import.h"
33#include "Epetra_CrsMatrix.h"
34#include "Epetra_Vector.h"
35#include "Epetra_Util.h"
36
37
38using namespace Teuchos;
39
40//=============================================================================
41Amesos_CSparse::Amesos_CSparse(const Epetra_LinearProblem &prob) :
42 UseTranspose_(false),
43 Problem_(&prob),
44 MtxConvTime_(-1),
45 MtxRedistTime_(-1),
46 VecRedistTime_(-1),
47 SymFactTime_(-1),
48 NumFactTime_(-1),
49 SolveTime_(-1)
50{
51 // Initialize data structures for CSparse
52}
53
54//=============================================================================
55Amesos_CSparse::~Amesos_CSparse()
56{
57 // Clean up
58
59 //AMESOS_CHK_ERRV(CheckError(error));
60 if ((verbose_ && PrintTiming_) || verbose_ == 2) PrintTiming();
61 if ((verbose_ && PrintStatus_) || verbose_ == 2) PrintStatus();
62}
63
64//=============================================================================
65int Amesos_CSparse::ConvertToSerial()
66{
67 ResetTimer();
68 const Epetra_RowMatrix* StdIndexMatrix_ ;
69 Epetra_CrsMatrix* CrsMatrixA_;
70
71 if ( Reindex_ ) {
72 if ( Matrix_ == 0 ) AMESOS_CHK_ERR(-7);
73#ifdef HAVE_AMESOS_EPETRAEXT
74 const Epetra_Map& OriginalMap = Matrix_->RowMatrixRowMap();
75
76 const Epetra_Map &OriginalDomainMap =
77 UseTranspose()?GetProblem()->GetOperator()->OperatorRangeMap():
78 GetProblem()->GetOperator()->OperatorDomainMap();
79 const Epetra_Map &OriginalRangeMap =
80 UseTranspose()?GetProblem()->GetOperator()->OperatorDomainMap():
81 GetProblem()->GetOperator()->OperatorRangeMap();
82
83 StdIndex_ = rcp( new Amesos_StandardIndex( OriginalMap ) );
84 StdIndexDomain_ = rcp( new Amesos_StandardIndex( OriginalDomainMap ) );
85 StdIndexRange_ = rcp( new Amesos_StandardIndex( OriginalRangeMap ) );
86
87 CrsMatrixA_ = dynamic_cast<Epetra_CrsMatrix *>(Problem_->GetOperator());
88 StdIndexMatrix_ = StdIndex_->StandardizeIndex( CrsMatrixA_ );
89#else
90 std::cerr << "Amesos_CSparse requires EpetraExt to reindex matrices." << std::endl
92#endif
93 } else {
94 StdIndexMatrix_ = Matrix_ ;
95 }
96
97 int NumGlobalRows = Matrix_->NumGlobalRows();
98
99 // create a serial map
100 int NumMyRows = 0;
101 if (Comm().MyPID() == 0)
102 NumMyRows = NumGlobalRows;
103
104 SerialMap_ = rcp(new Epetra_Map(-1, NumMyRows, 0, Comm()));
105 if (SerialMap_.get() == 0)
106 AMESOS_CHK_ERR(-1);
107
108 Importer_ = rcp(new Epetra_Import(SerialMap(),StdIndexMatrix_->RowMatrixRowMap()));
109 if (Importer_.get() == 0)
110 AMESOS_CHK_ERR(-1);
111
112 SerialCrsMatrix_ = rcp(new Epetra_CrsMatrix(Copy, SerialMap(), 0));
113 if (SerialCrsMatrix_.get() == 0)
114 AMESOS_CHK_ERR(-1);
115
116 AMESOS_CHK_ERR(SerialCrsMatrix().Import(*StdIndexMatrix_, Importer(), Add));
117
118 AMESOS_CHK_ERR(SerialCrsMatrix().FillComplete());
119
120 SerialMatrix_ = rcp(SerialCrsMatrix_.get(), false);
121
122 MtxRedistTime_ = AddTime("Total matrix redistribution time", MtxRedistTime_);
123
124 return 0;
125}
126
127//=============================================================================
128int Amesos_CSparse::ConvertToCSparse()
129{
130
131#ifdef HAVE_AMESOS_CSPARSE
132
133 ResetTimer();
134
135 if (Comm().MyPID() == 0)
136 {
137 csMatrix.p = (ptrdiff_t *) malloc((SerialMatrix().NumMyRows()+1)*sizeof(ptrdiff_t));
138 csMatrix.i = (ptrdiff_t *) malloc(SerialMatrix().NumMyNonzeros()*sizeof(ptrdiff_t));
139 csMatrix.x = (double *) malloc(SerialMatrix().NumMyNonzeros()*sizeof(double));
140 csMatrix.nzmax = SerialMatrix().NumMyNonzeros();
141 csMatrix.m = SerialMatrix().NumMyRows();
142 csMatrix.n = SerialMatrix().NumMyRows();
143 csMatrix.nz = -1;
144
145 int MaxNumEntries = SerialMatrix().MaxNumEntries();
146 std::vector<int> Indices(MaxNumEntries);
147 std::vector<double> Values(MaxNumEntries);
148
149 csMatrix.p[0] = 0;
150 int count = 0;
151
152 for (int i = 0 ; i < SerialMatrix().NumMyRows() ; ++i)
153 {
154 int ierr, NumEntriesThisRow;
155 ierr = SerialMatrix().ExtractMyRowCopy(i, MaxNumEntries,
156 NumEntriesThisRow,
157 &Values[0], &Indices[0]);
158 if (ierr < 0)
159 AMESOS_CHK_ERR(ierr);
160
161 csMatrix.p[i + 1] = csMatrix.p[i] + NumEntriesThisRow;
162
163 for (int j = 0 ; j < NumEntriesThisRow ; ++j)
164 {
165 if (Indices[j] == i)
166 Values[j] += AddToDiag_;
167
168 csMatrix.i[count] = Indices[j];
169 csMatrix.x[count] = Values[j];
170 ++count;
171 }
172 }
173
174 if (count != SerialMatrix().NumMyNonzeros())
175 AMESOS_CHK_ERR(-1); // something wrong here*/
176
177 // TODO: Avoid this transpose once we have cs_sptsolve()
178 csTranMatrix = cs_transpose(&csMatrix, 1);
179
180 free(csMatrix.p);
181 free(csMatrix.i);
182 free(csMatrix.x);
183 }
184
185 MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
186
187 return 0;
188#else
189 AMESOS_CHK_ERR(-1); // Don't have CSPARSE
190 return 1;
191#endif
192}
193
194//=============================================================================
195int Amesos_CSparse::SetParameters( Teuchos::ParameterList &ParameterList)
196{
197 // retrive general parameters
198
199 SetStatusParameters( ParameterList );
200
201 SetControlParameters( ParameterList );
202
203 // retrive CSparse's specific parameters
204
205 if (ParameterList.isSublist("CSparse"))
206 {
207 // set solver specific parameters here.
208
209 }
210
211 return 0;
212}
213
214//=============================================================================
215int Amesos_CSparse::PerformSymbolicFactorization()
216{
217#ifdef HAVE_AMESOS_CSPARSE
218
219 ResetTimer();
220
221 if (Comm().MyPID() == 0)
222 {
223 // ============================================================== //
224 // Setup CSparse parameteres and call symbolic factorization.
225 // ============================================================== //
226 int error = 0;
227
228 // Call Csparse here.
229 csSymbolic = cs_sqr(2, csTranMatrix, 0);
230
232 }
233
234 SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_);
235
236 return 0;
237#else
238 AMESOS_CHK_ERR(-1); // Don't have CSPARSE
239 return 1;
240#endif
241}
242
243//=============================================================================
244int Amesos_CSparse::PerformNumericFactorization( )
245{
246#ifdef HAVE_AMESOS_CSPARSE
247 ResetTimer();
248
249 if (Comm().MyPID() == 0)
250 {
251 int error = 0;
252 // Call Csparse here.
253 csNumeric = cs_lu(csTranMatrix, csSymbolic, 1e-15);
254
256 }
257
258 NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_);
259
260 return 0;
261#else
262 AMESOS_CHK_ERR(-1); // Don't have CSPARSE
263 return 1;
264#endif
265}
266
267//=============================================================================
268bool Amesos_CSparse::MatrixShapeOK() const
269{
270 bool OK = true;
271
272 if (GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
273 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() )
274 {
275 OK = false;
276 }
277 return OK;
278}
279
280//=============================================================================
281int Amesos_CSparse::SymbolicFactorization()
282{
283 IsSymbolicFactorizationOK_ = false;
284 IsNumericFactorizationOK_ = false;
285
286 CreateTimer(Comm());
287
288 ++NumSymbolicFact_;
289
290 Matrix_ = dynamic_cast<Epetra_RowMatrix*>(Problem_->GetOperator());
291 Map_ = &(Matrix_->RowMatrixRowMap());
292
293 // =========================================================== //
294 // redistribute and create all import/export objects only //
295 // if more than one processor is used. Otherwise simply set //
296 // dummy pointers to Matrix() and Map(), without giving the //
297 // ownership to the smart pointer. //
298 // =========================================================== //
299
300 if (Comm().NumProc() != 1)
301 ConvertToSerial();
302 else
303 {
304 SerialMap_ = rcp(const_cast<Epetra_Map*>(&Map()), false);
305 SerialMatrix_ = rcp(const_cast<Epetra_RowMatrix*>(&Matrix()), false);
306 }
307
308 // =========================================================== //
309 // Only on processor zero, convert the matrix into CSR format, //
310 // as required by CSparse. //
311 // =========================================================== //
312
313 ConvertToCSparse();
314
315 PerformSymbolicFactorization();
316
317 IsSymbolicFactorizationOK_ = true;
318
319 return(0);
320}
321
322//=============================================================================
323int Amesos_CSparse::NumericFactorization()
324{
325 IsNumericFactorizationOK_ = false;
326
327 if (IsSymbolicFactorizationOK_ == false)
328 AMESOS_CHK_ERR(SymbolicFactorization());
329
330 ++NumNumericFact_;
331
332 // FIXME: this must be checked, now all the matrix is shipped twice here
333 ConvertToSerial();
334 ConvertToCSparse();
335
336 PerformNumericFactorization();
337
338 IsNumericFactorizationOK_ = true;
339
340 return(0);
341}
342
343//=============================================================================
344int Amesos_CSparse::Solve()
345{
346#ifdef HAVE_AMESOS_CSPARSE
347 Epetra_MultiVector* vecX = 0 ;
348 Epetra_MultiVector* vecB = 0 ;
349
350#ifdef HAVE_AMESOS_EPETRAEXT
351 Teuchos::RCP<Epetra_MultiVector> vecX_rcp;
352 Teuchos::RCP<Epetra_MultiVector> vecB_rcp;
353#endif
354
355 if (IsNumericFactorizationOK_ == false)
356 AMESOS_CHK_ERR(NumericFactorization());
357
358 Epetra_MultiVector* X = Problem_->GetLHS();
359 Epetra_MultiVector* B = Problem_->GetRHS();
360
361 if ((X == 0) || (B == 0))
362 AMESOS_CHK_ERR(-1);
363
364 int NumVectors = X->NumVectors();
365 if (NumVectors != B->NumVectors())
366 AMESOS_CHK_ERR(-1);
367
368 if ( Reindex_ ) {
369#ifdef HAVE_AMESOS_EPETRAEXT
370 vecX_rcp = StdIndexDomain_->StandardizeIndex( *X ) ;
371 vecB_rcp = StdIndexRange_->StandardizeIndex( *B ) ;
372
373 vecX = &*vecX_rcp;
374 vecB = &*vecB_rcp;
375#else
376 AMESOS_CHK_ERR( -13 ) ; // Amesos_CSparse can't handle non-standard indexing without EpetraExt
377#endif
378 } else {
379 vecX = X ;
380 vecB = B ;
381 }
382
383 // vectors with SerialMap_
384 Epetra_MultiVector* SerialB;
385 Epetra_MultiVector* SerialX;
386
387 ResetTimer();
388
389 SerialX = new Epetra_MultiVector(SerialMap(),NumVectors);
390 SerialB = new Epetra_MultiVector(SerialMap(),NumVectors);
391
392 SerialB->Import(*vecB,Importer(),Insert);
393
394 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
395
396 ResetTimer();
397
398 if (Comm().MyPID() == 0)
399 {
400 double* SerialXValues;
401 double* SerialBValues;
402 int LDA;
403
404 AMESOS_CHK_ERR(SerialX->ExtractView(&SerialXValues,&LDA));
405
406 // FIXME: check LDA
407 AMESOS_CHK_ERR(SerialB->ExtractView(&SerialBValues,&LDA));
408
409 int error = 0;
410 int n = SerialMatrix().NumMyRows();
411
412 double *x = (double *) malloc(n * sizeof(double));
413 // TODO: OMP here ??
414 for (int i = 0 ; i < NumVectors ; ++i)
415 {
416 // Call Csparse here
417 cs_ipvec(csNumeric->pinv, SerialBValues+i*n, x, n);
418 cs_lsolve(csNumeric->L, x);
419 cs_usolve(csNumeric->U, x);
420 cs_ipvec(csSymbolic->q, x, SerialXValues+i*n, n);
421 }
422 free(x);
423
425 }
426
427 SolveTime_ = AddTime("Total solve time", SolveTime_);
428
429 // Copy X back to the original vector
430
431 ResetTimer();
432
433 vecX->Export(*SerialX, Importer(), Insert);
434 delete SerialB;
435 delete SerialX;
436
437 VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_);
438
439 if (ComputeTrueResidual_)
440 ComputeTrueResidual(Matrix(), *X, *B, UseTranspose(), "Amesos_CSparse");
441
442 if (ComputeVectorNorms_)
443 ComputeVectorNorms(*X, *B, "Amesos_CSparse");
444
445 ++NumSolve_;
446
447 return(0) ;
448#else
449 AMESOS_CHK_ERR(-1); // Don't have CSPARSE
450 return 1;
451#endif
452}
453
454// ======================================================================
455void Amesos_CSparse::PrintStatus() const
456{
457 if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
458 return;
459
460 std::string p = "Amesos_CSparse : ";
461 PrintLine();
462
463 PrintLine();
464
465 return;
466}
467
468// ======================================================================
469void Amesos_CSparse::PrintTiming() const
470{
471 if (Problem_->GetOperator() == 0 || Comm().MyPID() != 0)
472 return;
473
474 double ConTime = GetTime(MtxConvTime_);
475 double MatTime = GetTime(MtxRedistTime_);
476 double VecTime = GetTime(VecRedistTime_);
477 double SymTime = GetTime(SymFactTime_);
478 double NumTime = GetTime(NumFactTime_);
479 double SolTime = GetTime(SolveTime_);
480
481 if (NumSymbolicFact_)
482 SymTime /= NumSymbolicFact_;
483
484 if (NumNumericFact_)
485 NumTime /= NumNumericFact_;
486
487 if (NumSolve_)
488 SolTime /= NumSolve_;
489
490 std::string p = "Amesos_CSparse : ";
491 PrintLine();
492
493 std::cout << p << "Time to convert matrix to Csparse format = "
494 << ConTime << " (s)" << std::endl;
495 std::cout << p << "Time to redistribute matrix = "
496 << MatTime << " (s)" << std::endl;
497 std::cout << p << "Time to redistribute vectors = "
498 << VecTime << " (s)" << std::endl;
499 std::cout << p << "Number of symbolic factorizations = "
500 << NumSymbolicFact_ << std::endl;
501 std::cout << p << "Time for sym fact = "
502 << SymTime << " (s), avg = " << SymTime << " (s)" << std::endl;
503 std::cout << p << "Number of numeric factorizations = "
504 << NumNumericFact_ << std::endl;
505 std::cout << p << "Time for num fact = "
506 << NumTime << " (s), avg = " << NumTime << " (s)" << std::endl;
507 std::cout << p << "Number of solve phases = "
508 << NumSolve_ << std::endl;
509 std::cout << p << "Time for solve = "
510 << SolTime << " (s), avg = " << SolTime << " (s)" << std::endl;
511
512 PrintLine();
513
514 return;
515}
516
517// ======================================================================
518int Amesos_CSparse::CheckError(const int error) const
519{
520 if (!error)
521 return 0;
522
523 std::cerr << "Amesos: CSparse returned error code " << error << std::endl;
524
525 AMESOS_RETURN(error);
526}
527
528#endif
#define AMESOS_CHK_ERR(a)
#define AMESOS_RETURN(amesos_err)
bool CheckError(const std::string SolverType, const std::string Descriptor, const Epetra_RowMatrix &A, const Epetra_MultiVector &x, const Epetra_MultiVector &b, const Epetra_MultiVector &x_exact)