126 if(
debug_ == 1 ) std::cout <<
"Entering `RedistributeA()'" << std::endl;
128 Time_->ResetStartTime();
130 Epetra_RowMatrix *RowMatrixA =
dynamic_cast<Epetra_RowMatrix *
>(
Problem_->GetOperator());
133 const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ;
134 int NumberOfProcesses =
Comm().NumProc() ;
142 int NumRows_ = RowMatrixA->NumGlobalRows() ;
143 int NumColumns_ = RowMatrixA->NumGlobalCols() ;
145 NumberOfProcesses = EPETRA_MIN( NumberOfProcesses,
MaxProcesses_ ) ;
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;
164 npcol_ = EPETRA_MIN( NumberOfProcesses,
165 EPETRA_MAX ( 2, (
int) sqrt( NumberOfProcesses * 0.5 ) ) ) ;
171 int NumMyElements = RowMatrixA->RowMatrixRowMap().NumMyElements() ;
172 std::vector<int> MyGlobalElements( NumMyElements );
173 RowMatrixA->RowMatrixRowMap().MyGlobalElements( &MyGlobalElements[0] ) ;
175 int NumMyColumns = RowMatrixA->RowMatrixColMap().NumMyElements() ;
176 std::vector<int> MyGlobalColumns( NumMyColumns );
177 RowMatrixA->RowMatrixColMap().MyGlobalElements( &MyGlobalColumns[0] ) ;
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;
189 Epetra_Map FatInMap(
npcol_*NumRows_, NumMyElements*
npcol_,
190 &MyFatElements[0], 0,
Comm() );
195 Epetra_CrsMatrix FatIn( Copy, FatInMap, 0 );
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;
215 int MaxNumIndices = RowMatrixA->MaxNumEntries();
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 ) {
228 RowMatrixA->ExtractMyRowCopy( 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++ ) {
257 FatIn.InsertGlobalValues( MyGlobalElements[LocalRow]*
npcol_ + pcol_i,
258 FatRowPtrs[ pcol_i ],
259 &FatMatrixValues[ pcol_i ][0],
260 &FatColumnIndices[ pcol_i ][0] );
263 FatIn.FillComplete(
false );
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_ ;
279 int OurExcessRows = EPETRA_MIN(
nb_, AllExcessRows - (
myprow_ *
nb_ ) ) ;
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_ ;
286 int OurExcessColumns = EPETRA_MIN(
nb_, AllExcessColumns - (
mypcol_ *
nb_ ) ) ;
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;
383 FatOut_ =
new Epetra_CrsMatrix( Copy, FatOutMap, 0 ) ;
385 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:348" << std::endl;
387 Epetra_Export ExportToFatOut( FatInMap, FatOutMap ) ;
389 if (
debug_ == 1) std::cout <<
"iam_ = " <<
iam_ <<
"Amesos_Scalapack.cpp:360" << std::endl;
391 FatOut_->Export( FatIn, ExportToFatOut, Add );
392 FatOut_->FillComplete(
false );
397 Epetra_RowMatrix *RowMatrixA =
dynamic_cast<Epetra_RowMatrix *
>(
Problem_->GetOperator());
398 const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ;
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 ;
422 int MyFirstElement = EPETRA_MIN(
iam_ *
m_per_p_, NumRows_ ) ;
423 int MyFirstNonElement = EPETRA_MIN( (
iam_+1) *
m_per_p_, NumRows_ ) ;
424 int NumExpectedElements = MyFirstNonElement - MyFirstElement ;
426 assert( NumRows_ == RowMatrixA->NumGlobalRows() ) ;
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;
544 Time_->ResetStartTime();
550 for (
int i = 0 ; i < (int)
DenseA_.size() ; i++ )
DenseA_[i] = 0 ;
559 int MaxNumEntries =
FatOut_->MaxNumEntries();
561 std::vector<int>ColIndicesV(MaxNumEntries);
562 std::vector<double>RowValuesV(MaxNumEntries);
564 int NumMyElements =
FatOut_->NumMyRows() ;
565 for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
567 ExtractMyRowView( MyRow, NzThisRow, RowValues,
568 ColIndices ) != 0 ) ;
572 int MyGlobalRow =
FatOut_->GRID( MyRow );
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++ ) {
593 assert(
FatOut_->RowMatrixColMap().GID( ColIndices[j] ) ==
594 FatOut_->GCID( ColIndices[j] ) );
596 int MyGlobalCol =
FatOut_->GCID( ColIndices[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++ ) {
641 Epetra_RowMatrix *RowMatrixA =
dynamic_cast<Epetra_RowMatrix *
>(
Problem_->GetOperator());
642 const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ;
644 int NumMyElements_ = 0 ;
812 if(
debug_ == 1 ) std::cout <<
"Entering `Solve()'" << std::endl;
816 Epetra_MultiVector *vecX =
Problem_->GetLHS() ;
817 Epetra_MultiVector *vecB =
Problem_->GetRHS() ;
828 nrhs = vecX->NumVectors() ;
832 Epetra_MultiVector *ScalapackB =0;
833 Epetra_MultiVector *ScalapackX =0;
837 double *ScalapackXvalues ;
839 Epetra_RowMatrix *RowMatrixA =
dynamic_cast<Epetra_RowMatrix *
>(
Problem_->GetOperator());
840 Time_->ResetStartTime();
844 const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap();
845 Epetra_MultiVector *ScalapackXextract =
new Epetra_MultiVector( *
VectorMap_, nrhs ) ;
846 Epetra_MultiVector *ScalapackBextract =
new Epetra_MultiVector( *
VectorMap_, nrhs ) ;
848 Epetra_Import ImportToScalapack( *
VectorMap_, OriginalMap );
849 ScalapackBextract->Import( *vecB, ImportToScalapack, Insert ) ;
850 ScalapackB = ScalapackBextract ;
851 ScalapackX = ScalapackXextract ;
861 ScalapackX->Scale(1.0, *ScalapackB) ;
865 Time_->ResetStartTime();
880 EPETRA_CHK_ERR( ScalapackX->ExtractView( &ScalapackXvalues, &ScalapackXlda ) ) ;
882 if (
false ) std::cout <<
"Amesos_Scalapack.cpp: " << __LINE__ <<
" ScalapackXlda = " << ScalapackXlda
883 <<
" lda_ = " <<
lda_
888 <<
" iam_ = " <<
iam_ << std::endl ;
901 assert( Ierr[0] == 0 ) ;
944 Time_->ResetStartTime();
948 Epetra_Import ImportFromScalapack( OriginalMap, *
VectorMap_ );
949 vecX->Import( *ScalapackX, ImportFromScalapack, Insert ) ;
950 delete ScalapackBextract ;
951 delete ScalapackXextract ;
957 Comm().Broadcast( Ierr, 1, 0 ) ;
961 double NormLHS, NormRHS;
962 for(
int i=0 ; i<nrhs ; ++i ) {
966 std::cout <<
"Amesos_Scalapack : vector " << i <<
", ||x|| = " << NormLHS
967 <<
", ||b|| = " << NormRHS << std::endl;
975 Epetra_MultiVector Ax(vecB->Map(),nrhs);
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;