45 inline int pcolnum(
int j,
int nb,
int npcol ) {
46 return ((j/nb)%npcol) ; }
53 ScaLAPACK1DMatrix_(0),
126 if(
debug_ == 1 ) std::cout <<
"Entering `RedistributeA()'" << std::endl;
148 int ProcessNumHeuristic = (1+NumRows_/200)*(1+NumRows_/200);
149 NumberOfProcesses =
EPETRA_MIN( NumberOfProcesses, ProcessNumHeuristic );
152 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:171" << std::endl;
159 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:163" << std::endl;
162 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:166" << std::endl;
165 EPETRA_MAX ( 2, (
int) sqrt( NumberOfProcesses * 0.5 ) ) ) ;
172 std::vector<int> MyGlobalElements( NumMyElements );
176 std::vector<int> MyGlobalColumns( NumMyColumns );
179 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:194" << std::endl;
181 std::vector<int> MyFatElements( NumMyElements *
npcol_ );
183 for(
int LocalRow=0; LocalRow<NumMyElements; LocalRow++ ) {
184 for (
int i = 0 ; i <
npcol_; i++ ){
185 MyFatElements[LocalRow*
npcol_+i] = MyGlobalElements[LocalRow]*
npcol_+i;
190 &MyFatElements[0], 0,
Comm() );
198 std::vector<std::vector<int> > FatColumnIndices(
npcol_,std::vector<int>(1));
199 std::vector<std::vector<double> > FatMatrixValues(
npcol_,std::vector<double>(1));
200 std::vector<int> FatRowPtrs(
npcol_);
203 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:219" << std::endl;
217 std::vector<int> ColumnIndices(MaxNumIndices);
218 std::vector<double> MatrixValues(MaxNumIndices);
220 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:232 NumMyElements = " 226 for(
int LocalRow=0; LocalRow<NumMyElements; ++LocalRow ) {
234 for (
int i=0; i<
npcol_; i++ ) FatRowPtrs[i] = 0 ;
240 for(
int i=0 ; i<NumIndices ; ++i ) {
241 int GlobalCol = MyGlobalColumns[ ColumnIndices[i] ];
243 if ( FatRowPtrs[ pcol_i ]+1 >= FatColumnIndices[ pcol_i ].size() ) {
244 FatColumnIndices[ pcol_i ]. resize( 2 * FatRowPtrs[ pcol_i ]+1 );
245 FatMatrixValues[ pcol_i ]. resize( 2 * FatRowPtrs[ pcol_i ]+1 );
247 FatColumnIndices[pcol_i][FatRowPtrs[pcol_i]] = GlobalCol ;
248 FatMatrixValues[pcol_i][FatRowPtrs[pcol_i]] = MatrixValues[i];
250 FatRowPtrs[ pcol_i ]++;
256 for (
int pcol_i = 0 ; pcol_i <
npcol_ ; pcol_i++ ) {
258 FatRowPtrs[ pcol_i ],
259 &FatMatrixValues[ pcol_i ][0],
260 &FatColumnIndices[ pcol_i ][0] );
265 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:260" << std::endl;
266 if (
debug_ == 1) std::cout <<
"Amesos_Scalapack.cpp:265B" 267 <<
" iam_ = " <<
iam_ 277 int UniformRows = ( NumRows_ / (
nprow_ *
nb_ ) ) *
nb_ ;
278 int AllExcessRows = NumRows_ - UniformRows *
nprow_ ;
280 OurExcessRows =
EPETRA_MAX( 0, OurExcessRows );
283 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:277" << std::endl;
284 int UniformColumns = ( NumColumns_ / (
npcol_ *
nb_ ) ) *
nb_ ;
285 int AllExcessColumns = NumColumns_ - UniformColumns *
npcol_ ;
287 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:281" << std::endl;
288 OurExcessColumns =
EPETRA_MAX( 0, OurExcessColumns );
299 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:295" << std::endl;
314 for ( ; BlockRow < UniformRows /
nb_ ; BlockRow++ ) {
315 for (
int RowOffset = 0; RowOffset <
nb_ ; RowOffset++ ) {
320 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:315" << std::endl;
321 assert ( BlockRow == UniformRows /
nb_ ) ;
322 for (
int RowOffset = 0; RowOffset < OurExcessRows ; RowOffset++ ) {
336 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:312" << std::endl;
343 for ( ; BlockRow < UniformRows /
nb_ ; BlockRow++ ) {
344 for (
int RowOffset = 0; RowOffset <
nb_ ; RowOffset++ ) {
351 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:326" << std::endl;
353 assert ( BlockRow == UniformRows /
nb_ ) ;
354 for (
int RowOffset = 0; RowOffset < OurExcessRows ; RowOffset++ ) {
360 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:334" << std::endl;
368 if (
debug_ == 1) std::cout <<
"Amesos_Scalapack.cpp:340" 369 <<
" iam_ = " <<
iam_ 372 <<
" NumRows_ = " << NumRows_
379 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:344" << std::endl;
385 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:348" << std::endl;
389 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:360" << std::endl;
400 int NumMyVecElements = 0 ;
405 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:385" << std::endl;
413 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
" Amesos_Scalapack.cpp:393 debug_ = " 421 m_per_p_ = ( NumRows_ + NumberOfProcesses - 1 ) / NumberOfProcesses ;
424 int NumExpectedElements = MyFirstNonElement - MyFirstElement ;
437 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
" Amesos_Scalapack.cpp:417 debug_ = " 439 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:402" 441 <<
" npcol_ = " <<
npcol_ << std::endl ;
446 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:408" << std::endl;
449 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:410A" << std::endl;
456 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"iam_ = " <<
iam_ <<
" Amesos_Scalapack.cpp:410" << std::endl;
458 assert( nprow ==
nprow_ ) ;
459 if ( npcol !=
npcol_ ) std::cout <<
"Amesos_Scalapack.cpp:430 npcol = " <<
460 npcol <<
" npcol_ = " <<
npcol_ << std::endl ;
461 assert( npcol ==
npcol_ ) ;
467 assert( myrow == 0 ) ;
468 assert( mycol ==
iam_ ) ;
472 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ 473 <<
"Amesos_Scalapack.cpp: " << __LINE__
478 <<
" lda_ = " <<
lda_ 483 <<
" iam_ = " <<
iam_ << std::endl ;
495 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:441" << std::endl;
496 assert( info == 0 ) ;
499 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:458 nprow = " << nprow << std::endl;
500 assert( nprow == -1 ) ;
503 if (
debug_ == 1) std::cout <<
"Amesos_Scalapack.cpp:446" << std::endl;
542 if(
debug_ == 1 ) std::cout <<
"Entering `ConvertToScalapack()'" << std::endl;
550 for (
int i = 0 ; i < (int)
DenseA_.size() ; i++ )
DenseA_[i] = 0 ;
561 std::vector<int>ColIndicesV(MaxNumEntries);
562 std::vector<double>RowValuesV(MaxNumEntries);
565 for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
567 ExtractMyRowView( MyRow, NzThisRow, RowValues,
568 ColIndices ) != 0 ) ;
574 int MyTrueRow = MyGlobalRow/
npcol_ ;
575 int UniformRows = ( MyTrueRow / (
nprow_ *
nb_ ) ) *
nb_ ;
576 int AllExcessRows = MyTrueRow - UniformRows *
nprow_ ;
577 int OurExcessRows = AllExcessRows - (
myprow_ *
nb_ ) ;
579 if ( MyRow != UniformRows + OurExcessRows ) {
580 std::cout <<
" iam _ = " <<
iam_ 581 <<
" MyGlobalRow = " << MyGlobalRow
582 <<
" MyTrueRow = " << MyTrueRow
583 <<
" UniformRows = " << UniformRows
584 <<
" AllExcessRows = " << AllExcessRows
585 <<
" OurExcessRows = " << OurExcessRows
586 <<
" MyRow = " << MyRow << std::endl ;
589 assert( OurExcessRows >= 0 && OurExcessRows <
nb_ );
590 assert( MyRow == UniformRows + OurExcessRows ) ;
592 for (
int j = 0; j < NzThisRow; j++ ) {
598 int UniformCols = ( MyGlobalCol / (
npcol_ *
nb_ ) ) *
nb_ ;
599 int AllExcessCols = MyGlobalCol - UniformCols *
npcol_ ;
600 int OurExcessCols = AllExcessCols - (
mypcol_ *
nb_ ) ;
601 assert( OurExcessCols >= 0 && OurExcessCols <
nb_ );
602 int MyCol = UniformCols + OurExcessCols ;
615 for (
int i = 0 ; i < (int)
DenseA_.size() ; i++ )
DenseA_[i] = 0 ;
625 std::vector<int>ColIndicesV(MaxNumEntries);
626 std::vector<double>RowValuesV(MaxNumEntries);
628 for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
630 ExtractMyRowView( MyRow, NzThisRow, RowValues,
631 ColIndices ) != 0 ) ;
633 for (
int j = 0; j < NzThisRow; j++ ) {
644 int NumMyElements_ = 0 ;
659 if(
debug_ == 1 ) std::cout <<
"Entering `SetParameters()'" << std::endl;
685 if ( ScalapackParams.
isParameter(
"2D distribution") )
696 if(
debug_ == 1 ) std::cout <<
"Entering `PerformNumericFactorization()'" << std::endl;
704 if (
false) std::cout <<
" Amesos_Scalapack.cpp: 711 iam_ = " <<
iam_ <<
" DescA = " 719 std::cout <<
" DenseA = " << std::endl ;
755 if (
debug_ == 1 && Ierr[0] != 0 )
756 std::cout <<
" Amesos_Scalapack.cpp:738 iam_ = " <<
iam_ 757 <<
" Ierr = " << Ierr[0] << std::endl ;
772 if (
GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() !=
773 GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) OK =
false;
780 if(
debug_ == 1 ) std::cout <<
"Entering `PerformSymbolicFactorization()'" << std::endl;
791 if(
debug_ == 1 ) std::cout <<
"Entering `NumericFactorization()'" << std::endl;
812 if(
debug_ == 1 ) std::cout <<
"Entering `Solve()'" << std::endl;
837 double *ScalapackXvalues ;
849 ScalapackBextract->
Import( *vecB, ImportToScalapack,
Insert ) ;
850 ScalapackB = ScalapackBextract ;
851 ScalapackX = ScalapackXextract ;
861 ScalapackX->
Scale(1.0, *ScalapackB) ;
880 assert( ScalapackX->
ExtractView( &ScalapackXvalues, &ScalapackXlda ) == 0 ) ;
882 if (
false ) std::cout <<
"Amesos_Scalapack.cpp: " << __LINE__ <<
" ScalapackXlda = " << ScalapackXlda
883 <<
" lda_ = " <<
lda_ 888 <<
" iam_ = " <<
iam_ << std::endl ;
901 assert( Ierr[0] == 0 ) ;
949 vecX->
Import( *ScalapackX, ImportFromScalapack,
Insert ) ;
950 delete ScalapackBextract ;
951 delete ScalapackXextract ;
961 double NormLHS, NormRHS;
962 for(
int i=0 ; i<nrhs ; ++i ) {
963 assert((*vecX)(i)->Norm2(&NormLHS)==0);
964 assert((*vecB)(i)->Norm2(&NormRHS)==0);
966 std::cout <<
"Amesos_Scalapack : vector " << i <<
", ||x|| = " << NormLHS
967 <<
", ||b|| = " << NormRHS << std::endl;
976 for(
int i=0 ; i<nrhs ; ++i ) {
978 (Ax.Update(1.0, *((*vecB)(i)), -1.0));
982 std::cout <<
"Amesos_Scalapack : vector " << i <<
", ||Ax - b|| = " << Norm << std::endl;
997 if(
iam_ != 0 )
return;
999 std::cout <<
"----------------------------------------------------------------------------" << std::endl;
1002 std::cout <<
"Amesos_Scalapack : Nonzero elements per row = " 1004 std::cout <<
"Amesos_Scalapack : Percentage of nonzero elements = " 1006 std::cout <<
"Amesos_Scalapack : Use transpose = " <<
UseTranspose_ << std::endl;
1007 std::cout <<
"----------------------------------------------------------------------------" << std::endl;
1019 std::cout <<
"----------------------------------------------------------------------------" << std::endl;
1020 std::cout <<
"Amesos_Scalapack : Time to convert matrix to ScaLAPACK format = " 1021 <<
ConTime_ <<
" (s)" << std::endl;
1022 std::cout <<
"Amesos_Scalapack : Time to redistribute matrix = " 1023 <<
MatTime_ <<
" (s)" << std::endl;
1024 std::cout <<
"Amesos_Scalapack : Time to redistribute vectors = " 1025 <<
VecTime_ <<
" (s)" << std::endl;
1026 std::cout <<
"Amesos_Scalapack : Number of symbolic factorizations = " 1028 std::cout <<
"Amesos_Scalapack : Time for sym fact = " 1030 <<
" (s)" << std::endl;
1031 std::cout <<
"Amesos_Scalapack : Number of numeric factorizations = " 1033 std::cout <<
"Amesos_Scalapack : Time for num fact = " 1035 <<
" (s)" << std::endl;
1036 std::cout <<
"Amesos_Scalapack : Number of solve phases = " 1038 std::cout <<
"Amesos_Scalapack : Time for solve = " 1040 <<
" (s)" << std::endl;
1041 std::cout <<
"----------------------------------------------------------------------------" << std::endl;
int NumSymbolicFact_
Number of symbolic factorization phases.
void PREFIX BLACS_GRIDINFO_F77(int *blacs_context, const int *nprow, const int *npcol, const int *myrow, const int *mycol)
const Epetra_BlockMap & Map() const
int MaxNumEntries() const
int MyGlobalElements(int *MyGlobalElementList) const
void PREFIX PDGETRS_F77(Epetra_fcd, const int *n, const int *nrhs, const double *A, const int *Ai, const int *Aj, const int *DescA, const int *ipiv, double *X, const int *Xi, const int *Xj, const int *DescX, int *info)
virtual const Epetra_Map & RowMatrixRowMap() const=0
int SymbolicFactorization()
Performs SymbolicFactorization on the matrix A.
void PrintTiming() const
Print timing information.
T & get(ParameterList &l, const std::string &name)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
bool isSublist(const std::string &name) const
void PREFIX DESCINIT_F77(int *DescA, const int *m, const int *n, const int *mblock, const int *nblock, const int *rsrc, const int *csrc, const int *blacs_context, const int *Lda, int *ierr)
RCP< ParameterList > sublist(const RCP< ParameterList > ¶mList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
int NumGlobalCols() const
virtual void Barrier() const=0
const Epetra_LinearProblem * GetProblem() const
Get a pointer to the Problem.
int Solve()
Solves A X = B (or AT X = B)
int NumNumericFact_
Number of numeric factorization phases.
virtual int NumGlobalNonzeros() const=0
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
int NumSolve_
Number of solves.
virtual int MaxNumEntries() const=0
int Scale(double ScalarValue)
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
void SetStatusParameters(const Teuchos::ParameterList &ParameterList)
bool ComputeVectorNorms_
If true, prints the norms of X and B in Solve().
bool UseTranspose() const
Returns the current UseTranspose setting.
void SetControlParameters(const Teuchos::ParameterList &ParameterList)
int GRID(int LRID_in) const
int NumMyElements() const
void PREFIX SL_INIT_F77(int *blacs_context, const int *nprow, const int *npcol)
bool PrintTiming_
If true, prints timing information in the destructor.
bool PrintStatus_
If true, print additional information in the destructor.
virtual int NumGlobalCols() const=0
#define EPETRA_CHK_ERR(xxx)
const Epetra_LinearProblem * Problem_
int GCID(int LCID_in) const
virtual int Broadcast(double *MyVals, int Count, int Root) const=0
const Epetra_Map & RowMatrixColMap() const
void PREFIX DGETRS_F77(Teuchos_fcd, const int *n, const int *nrhs, const double *a, const int *lda, const int *ipiv, double *x, const int *ldx, int *info)
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
int NumericFactorization()
Performs NumericFactorization on the matrix A.
#define AMESOS_PRINT(variable)
void PREFIX PDGETRF_F77(const int *m, const int *n, double *A, const int *Ai, const int *Aj, const int *DescA, int *ipiv, int *info)
Epetra_CrsMatrix * FatOut_
Epetra_MultiVector * GetRHS() const
virtual int NumProc() const=0
int NumGlobalElements() const
int SetParameters(Teuchos::ParameterList &ParameterList)
Updates internal variables.
bool MatrixShapeOK() const
Returns true if SCALAPACK can handle this matrix shape.
std::vector< double > DenseA_
bool isParameter(const std::string &name) const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
void PrintStatus() const
Print information about the factorization and solution phases.
Epetra_RowMatrix * GetMatrix() const
Epetra_MultiVector * GetLHS() const
virtual const Epetra_Map & RowMatrixColMap() const=0
int ExtractView(double **A, int *MyLDA) const
int verbose_
Toggles the output level.
int NumGlobalRows() const
Amesos_Scalapack(const Epetra_LinearProblem &LinearProblem)
Amesos_Scalapack Constructor.
int debug_
Sets the level of debug_ output.
int PerformNumericFactorization()
Epetra_Map * ScaLAPACK1DMap_
void ResetStartTime(void)
Epetra_Operator * GetOperator() const
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
int pcolnum(int j, int nb, int npcol)
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Epetra_CrsMatrix * ScaLAPACK1DMatrix_
void PREFIX DGETRF_F77(const int *m, const int *n, double *a, const int *lda, int *ipiv, int *info)
virtual int NumGlobalRows() const=0
~Amesos_Scalapack(void)
Amesos_Scalapack Destructor.
double ElapsedTime(void) const
bool ComputeTrueResidual_
If true, computes the true residual in Solve().