Amesos Package Browser (Single Doxygen Collection) Development
Amesos_Scalapack.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: Direct Sparse Solver Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Amesos_Scalapack.h"
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Import.h"
00032 #include "Epetra_Export.h"
00033 #include "Epetra_CrsMatrix.h"
00034 #include "Epetra_Vector.h"
00035 #include "Epetra_Util.h"
00036 #include "CrsMatrixTranspose.h"
00037 #include "Amesos_SCALAPACK_wrappers.h"
00038 //  #include "Epetra_LAPACK_wrappers.h"
00039 
00040 
00041 //
00042 //  pcolnum computes the process column which a given column belongs to
00043 //  in the ScaLAPACK 2D grid.
00044 //
00045 inline int pcolnum( int j, int nb, int npcol ) {
00046   return ((j/nb)%npcol) ; }
00047 
00048 
00049 //=============================================================================
00050 Amesos_Scalapack::Amesos_Scalapack(const Epetra_LinearProblem &prob ):
00051   ictxt_(-1313),
00052   ScaLAPACK1DMap_(0), 
00053   ScaLAPACK1DMatrix_(0), 
00054   VectorMap_(0),
00055   UseTranspose_(false),   // Overwritten by call to SetParameters below
00056   Problem_(&prob), 
00057   ConTime_(0.0),
00058   SymTime_(0.0),
00059   NumTime_(0.0),
00060   SolTime_(0.0),
00061   VecTime_(0.0),
00062   MatTime_(0.0),
00063   FatOut_(0),
00064   Time_(0)
00065 {
00066     Teuchos::ParameterList ParamList ;
00067     SetParameters( ParamList ) ; 
00068 }
00069 
00070 //=============================================================================
00071 Amesos_Scalapack::~Amesos_Scalapack(void) {
00072 
00073   if ( ScaLAPACK1DMap_ ) delete ScaLAPACK1DMap_ ; 
00074   if ( ScaLAPACK1DMatrix_ ) delete ScaLAPACK1DMatrix_ ; 
00075   // print out some information if required by the user
00076   if( (verbose_ && PrintTiming_) || verbose_ == 2 ) PrintTiming();
00077   if( (verbose_ && PrintStatus_) || verbose_ == 2 ) PrintStatus();
00078 
00079   if( Time_ ) delete Time_;
00080   if( FatOut_ ) delete FatOut_;
00081   if( VectorMap_ ) delete VectorMap_;
00082 }
00083 //  See  pre and post conditions in Amesos_Scalapack.h
00084 
00085 //
00086 //
00087 //  Distribution of the matrix:
00088 //    Amesos_Scalapack uses one of two data distributions:  
00089 //    1)  A 1D blocked data distribution in which each row is assigned
00090 //        to a single process.  The first (n/p) rows are assigned to 
00091 //        the first process and the i^th (n/p) rows are assigned to 
00092 //        process i.   A^T X = B is computed.
00093 //    2)  A 2D data distribution (A X = B is computed)
00094 //
00095 //    The 1D data distribution should be reasonably efficient for the matrices
00096 //    that we expect to deal with, i.e. n on the order of 4,000 
00097 //    with up to 16 processes.  For n >= 10,000 we should switch to 
00098 //    a 2D block-cyclis data distribution.  
00099 //
00100 //  Distribution of the vector(s)
00101 //    1)  1D data distribution
00102 //      ScaLAPACK requires that the vectors be distributed in the same manner
00103 //      as the matrices.  Since PDGETRF factors the transpose of the matrix, 
00104 //      using NPROW=1, the vectors must also be distributed with NPROW=1, i.e.
00105 //      each vector fully owned by a particular process.  And, the first nb
00106 //      vectors belong on the first process. 
00107 //
00108 //      The easiest way to deal with this is to limit ourselves to nb right hand 
00109 //      sides at a time (this is not a significant limitation as nb >= n/p ) 
00110 //      and using our basic heuristic for the number of processes to use 
00111 //      nb >= min(200,n) as well. 
00112 //   
00113 //      Limiting the number of right hand sides to <= nb means that all right hand 
00114 //      sides are stored on process 0.  Hence, they can be treated as a serial 
00115 //      matrix of vectors.
00116 //
00117 //    2)  2D data distribution
00118 //      If we restrict ourselves to nb (typically 32) right hand sides,
00119 //      we can use a simple epetra exporter to export the vector a single
00120 //      ScaLAPACK process column (i.e. to nprow processes)
00121 //
00122 
00123 //
00124 int Amesos_Scalapack::RedistributeA( ) {
00125 
00126   if( debug_ == 1 ) std::cout << "Entering `RedistributeA()'" << std::endl;
00127 
00128   Time_->ResetStartTime();
00129   
00130   Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
00131   EPETRA_CHK_ERR( RowMatrixA == 0 ) ; 
00132 
00133   const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
00134   int NumberOfProcesses = Comm().NumProc() ; 
00135 
00136   //
00137   //  Compute a uniform distribution as ScaLAPACK would want it
00138   //    MyFirstElement - The first element which this processor would have
00139   //    NumExpectedElemetns - The number of elements which this processor would have
00140   //
00141 
00142   int NumRows_ = RowMatrixA->NumGlobalRows() ; 
00143   int NumColumns_  = RowMatrixA->NumGlobalCols() ; 
00144   if ( MaxProcesses_ > 0 ) {
00145     NumberOfProcesses = EPETRA_MIN( NumberOfProcesses, MaxProcesses_ ) ; 
00146   }
00147   else {
00148     int ProcessNumHeuristic = (1+NumRows_/200)*(1+NumRows_/200);
00149     NumberOfProcesses = EPETRA_MIN( NumberOfProcesses,  ProcessNumHeuristic );
00150   }
00151   
00152   if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:171" << std::endl;
00153   //
00154   // Create the ScaLAPACK data distribution.
00155   // The TwoD data distribution is created in a completely different
00156   // manner and is not transposed (whereas the SaLAPACK 1D data
00157   // distribution was transposed) 
00158   //
00159   if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:163" << std::endl;
00160   Comm().Barrier(); 
00161   if ( TwoD_distribution_ ) { 
00162     if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:166" << std::endl;
00163     Comm().Barrier(); 
00164     npcol_ = EPETRA_MIN( NumberOfProcesses, 
00165        EPETRA_MAX ( 2, (int) sqrt( NumberOfProcesses * 0.5 ) ) ) ; 
00166     nprow_ = NumberOfProcesses / npcol_ ;
00167 
00168     //
00169     //  Create the map for FatA - our first intermediate matrix
00170     //
00171     int NumMyElements = RowMatrixA->RowMatrixRowMap().NumMyElements() ;
00172     std::vector<int> MyGlobalElements( NumMyElements );
00173     RowMatrixA->RowMatrixRowMap().MyGlobalElements( &MyGlobalElements[0] ) ;
00174 
00175     int NumMyColumns = RowMatrixA->RowMatrixColMap().NumMyElements() ;
00176     std::vector<int> MyGlobalColumns( NumMyColumns );
00177     RowMatrixA->RowMatrixColMap().MyGlobalElements( &MyGlobalColumns[0] ) ;
00178 
00179     if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:194" << std::endl;
00180 
00181     std::vector<int> MyFatElements( NumMyElements * npcol_ ); 
00182     
00183     for( int LocalRow=0; LocalRow<NumMyElements; LocalRow++ ) {
00184       for (int i = 0 ; i < npcol_; i++ ){
00185   MyFatElements[LocalRow*npcol_+i] = MyGlobalElements[LocalRow]*npcol_+i;
00186       } 
00187     }
00188     
00189     Epetra_Map FatInMap( npcol_*NumRows_, NumMyElements*npcol_, 
00190        &MyFatElements[0], 0, Comm() ); 
00191     
00192     //
00193     //  Create FatIn, our first intermediate matrix
00194     //
00195     Epetra_CrsMatrix FatIn( Copy, FatInMap, 0 );
00196     
00197     
00198     std::vector<std::vector<int> > FatColumnIndices(npcol_,std::vector<int>(1));
00199     std::vector<std::vector<double> > FatMatrixValues(npcol_,std::vector<double>(1));
00200     std::vector<int> FatRowPtrs(npcol_);  // A FatRowPtrs[i] = the number 
00201     // of entries in local row LocalRow*npcol_ + i 
00202     
00203     if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:219" << std::endl;
00204     //
00205     mypcol_ = iam_%npcol_;
00206     myprow_ = (iam_/npcol_)%nprow_;
00207     if ( iam_ >= nprow_ * npcol_ ) {
00208       myprow_ = nprow_; 
00209       mypcol_ = npcol_; 
00210     }
00211     //  Each row is split into npcol_ rows, with each of the 
00212     //  new rows containing only those elements belonging to 
00213     //  its process column (in the ScaLAPACK 2D process grid)
00214     //
00215     int MaxNumIndices = RowMatrixA->MaxNumEntries();
00216     int NumIndices;
00217     std::vector<int> ColumnIndices(MaxNumIndices);
00218     std::vector<double> MatrixValues(MaxNumIndices); 
00219     
00220     if ( debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:232 NumMyElements = " 
00221           << NumMyElements 
00222           << std::endl;
00223     
00224     nb_ = grid_nb_;
00225     
00226     for( int LocalRow=0; LocalRow<NumMyElements; ++LocalRow ) {
00227       
00228       RowMatrixA->ExtractMyRowCopy( LocalRow, 
00229             MaxNumIndices,
00230             NumIndices, 
00231             &MatrixValues[0],
00232             &ColumnIndices[0] );
00233       
00234       for (int i=0; i<npcol_; i++ )  FatRowPtrs[i] = 0 ; 
00235 
00236       //
00237       //  Deal the individual matrix entries out to the row owned by
00238       //  the process to which this matrix entry will belong.
00239       //
00240       for( int i=0 ; i<NumIndices ; ++i ) {
00241   int GlobalCol = MyGlobalColumns[ ColumnIndices[i] ];
00242   int pcol_i = pcolnum( GlobalCol, nb_, npcol_ ) ;
00243   if ( FatRowPtrs[ pcol_i ]+1 >= FatColumnIndices[ pcol_i ].size() ) {
00244     FatColumnIndices[ pcol_i ]. resize( 2 * FatRowPtrs[ pcol_i ]+1 );
00245     FatMatrixValues[ pcol_i ]. resize( 2 * FatRowPtrs[ pcol_i ]+1 );
00246   }
00247   FatColumnIndices[pcol_i][FatRowPtrs[pcol_i]] = GlobalCol ;
00248   FatMatrixValues[pcol_i][FatRowPtrs[pcol_i]] = MatrixValues[i];
00249   
00250   FatRowPtrs[ pcol_i ]++;
00251       }
00252       
00253       //
00254       //  Insert each of the npcol_ rows individually
00255       //
00256       for ( int pcol_i = 0 ; pcol_i < npcol_ ; pcol_i++ ) { 
00257   FatIn.InsertGlobalValues( MyGlobalElements[LocalRow]*npcol_ + pcol_i, 
00258           FatRowPtrs[ pcol_i ],
00259           &FatMatrixValues[ pcol_i ][0], 
00260           &FatColumnIndices[ pcol_i ][0] );
00261       }
00262     }
00263     FatIn.FillComplete( false );
00264     
00265     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:260" << std::endl;
00266     if (  debug_ == 1) std::cout  << "Amesos_Scalapack.cpp:265B" 
00267            << " iam_ = " << iam_ 
00268            << " nb_ = " << nb_ 
00269            << " nprow_ = " << nprow_ 
00270            << " npcol_ = " << npcol_ 
00271            << std::endl;
00272     
00273     //
00274     //  Compute the map for our second intermediate matrix, FatOut
00275     //
00276     //  Compute directly
00277     int UniformRows =  ( NumRows_ / ( nprow_ * nb_ ) ) * nb_ ; 
00278     int AllExcessRows = NumRows_ - UniformRows * nprow_ ; 
00279     int OurExcessRows = EPETRA_MIN( nb_, AllExcessRows - ( myprow_ * nb_ ) ) ; 
00280     OurExcessRows = EPETRA_MAX( 0, OurExcessRows );
00281     NumOurRows_ = UniformRows + OurExcessRows ; 
00282     
00283     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:277" << std::endl;
00284     int UniformColumns =  ( NumColumns_ / ( npcol_ * nb_ ) ) * nb_ ; 
00285     int AllExcessColumns = NumColumns_ - UniformColumns * npcol_ ; 
00286     int OurExcessColumns = EPETRA_MIN( nb_, AllExcessColumns - ( mypcol_ * nb_ ) ) ; 
00287     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:281" << std::endl;
00288     OurExcessColumns = EPETRA_MAX( 0, OurExcessColumns );
00289     NumOurColumns_ = UniformColumns + OurExcessColumns ; 
00290     
00291     if ( iam_ >= nprow_ * npcol_ ) {
00292       UniformRows = 0;
00293       NumOurRows_ = 0;
00294       NumOurColumns_ = 0;
00295     }
00296     
00297     Comm().Barrier(); 
00298     
00299     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:295" << std::endl;
00300 #if 0
00301     //  Compute using ScaLAPACK's numroc routine, assert agreement  
00302     int izero = 0; // All matrices start at process 0
00303     int NumRocSays = numroc_( &NumRows_, &nb_, &myprow_, &izero, &nprow_ );
00304     assert( NumOurRows_ == NumRocSays );
00305 #endif
00306     //
00307     //  Compute the rows which this process row owns in the ScaLAPACK 2D
00308     //  process grid.
00309     //
00310     std::vector<int> AllOurRows(NumOurRows_);
00311     
00312     int RowIndex = 0 ; 
00313     int BlockRow = 0 ;
00314     for ( ; BlockRow < UniformRows / nb_ ; BlockRow++ ) {
00315       for ( int RowOffset = 0; RowOffset < nb_ ; RowOffset++ ) {
00316   AllOurRows[RowIndex++] = BlockRow*nb_*nprow_  + myprow_*nb_ + RowOffset ;
00317       } 
00318     }
00319     Comm().Barrier(); 
00320     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:315" << std::endl;
00321     assert ( BlockRow == UniformRows / nb_ ) ; 
00322     for ( int RowOffset = 0; RowOffset < OurExcessRows ; RowOffset++ ) {
00323       AllOurRows[RowIndex++] = BlockRow*nb_*nprow_ + myprow_*nb_ + RowOffset ;
00324     } 
00325     assert( RowIndex == NumOurRows_ );
00326     //
00327     //  Distribute those rows amongst all the processes in that process row
00328     //  This is an artificial distribution with the following properties:
00329     //  1)  It is a 1D data distribution (each row belogs entirely to 
00330     //      a single process
00331     //  2)  All data which will eventually belong to a given process row, 
00332     //      is entirely contained within the processes in that row.
00333     //
00334     
00335     Comm().Barrier(); 
00336     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:312" << std::endl;
00337     //
00338     //  Compute MyRows directly
00339     //
00340     std::vector<int>MyRows(NumOurRows_);
00341     RowIndex = 0 ; 
00342     BlockRow = 0 ;
00343     for ( ; BlockRow < UniformRows / nb_ ; BlockRow++ ) {
00344       for ( int RowOffset = 0; RowOffset < nb_ ; RowOffset++ ) {
00345   MyRows[RowIndex++] = BlockRow*nb_*nprow_*npcol_  + 
00346     myprow_*nb_*npcol_ + RowOffset*npcol_  + mypcol_ ;
00347       } 
00348     }
00349     
00350     Comm().Barrier(); 
00351     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:326" << std::endl;
00352     
00353     assert ( BlockRow == UniformRows / nb_ ) ; 
00354     for ( int RowOffset = 0; RowOffset < OurExcessRows ; RowOffset++ ) {
00355       MyRows[RowIndex++] = BlockRow*nb_*nprow_*npcol_  + 
00356   myprow_*nb_*npcol_ + RowOffset*npcol_  + mypcol_ ;
00357     } 
00358     
00359     Comm().Barrier(); 
00360     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:334" << std::endl;
00361     Comm().Barrier(); 
00362     
00363     for (int i=0; i < NumOurRows_; i++ ) { 
00364       assert( MyRows[i] == AllOurRows[i]*npcol_+mypcol_ );
00365     } 
00366     
00367     Comm().Barrier(); 
00368     if (  debug_ == 1) std::cout  << "Amesos_Scalapack.cpp:340" 
00369            << " iam_ = " << iam_ 
00370            << " myprow_ = " << myprow_ 
00371            << " mypcol_ = " << mypcol_ 
00372            << " NumRows_ = " << NumRows_ 
00373            << " NumOurRows_ = " << NumOurRows_ 
00374            << std::endl;
00375     
00376     Comm().Barrier(); 
00377     Epetra_Map FatOutMap( npcol_*NumRows_, NumOurRows_, &MyRows[0], 0, Comm() ); 
00378     
00379     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:344" << std::endl;
00380     Comm().Barrier(); 
00381     
00382     if ( FatOut_ ) delete FatOut_ ; 
00383     FatOut_ = new Epetra_CrsMatrix( Copy, FatOutMap, 0 ) ;
00384     
00385     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:348" << std::endl;
00386     
00387     Epetra_Export ExportToFatOut( FatInMap, FatOutMap ) ;
00388     
00389     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:360" << std::endl;
00390     
00391     FatOut_->Export( FatIn, ExportToFatOut, Add );
00392     FatOut_->FillComplete( false );
00393     
00394     //
00395     //  Create a map to allow us to redistribute the vectors X and B 
00396     //
00397     Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
00398     const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
00399     assert( NumGlobalElements_ == OriginalMap.NumGlobalElements() ) ;
00400     int NumMyVecElements = 0 ;
00401     if ( mypcol_ == 0 ) { 
00402       NumMyVecElements = NumOurRows_;
00403     }
00404     
00405     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:385" << std::endl;
00406     
00407     if (VectorMap_) { delete VectorMap_ ; VectorMap_ = 0 ; } 
00408     VectorMap_ = new Epetra_Map( NumGlobalElements_, 
00409          NumMyVecElements, 
00410          &AllOurRows[0], 
00411          0, 
00412          Comm() );
00413     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << " Amesos_Scalapack.cpp:393 debug_ = "
00414            << debug_ << std::endl;
00415     
00416   } else {
00417     nprow_ = 1 ;
00418     npcol_ = NumberOfProcesses / nprow_ ;
00419     assert ( nprow_ * npcol_ == NumberOfProcesses ) ; 
00420     
00421     m_per_p_ = ( NumRows_ + NumberOfProcesses - 1 ) / NumberOfProcesses ;
00422     int MyFirstElement = EPETRA_MIN( iam_ * m_per_p_, NumRows_ ) ;
00423     int MyFirstNonElement = EPETRA_MIN( (iam_+1) * m_per_p_, NumRows_ ) ;
00424     int NumExpectedElements = MyFirstNonElement - MyFirstElement ; 
00425     
00426     assert( NumRows_ ==  RowMatrixA->NumGlobalRows() ) ; 
00427     if ( ScaLAPACK1DMap_ ) delete( ScaLAPACK1DMap_ ) ; 
00428     ScaLAPACK1DMap_ = new Epetra_Map( NumRows_, NumExpectedElements, 0, Comm() );
00429     if ( ScaLAPACK1DMatrix_ ) delete( ScaLAPACK1DMatrix_ ) ; 
00430     ScaLAPACK1DMatrix_ = new Epetra_CrsMatrix(Copy, *ScaLAPACK1DMap_, 0);
00431     Epetra_Export ExportToScaLAPACK1D_( OriginalMap, *ScaLAPACK1DMap_);
00432     
00433     ScaLAPACK1DMatrix_->Export( *RowMatrixA, ExportToScaLAPACK1D_, Add ); 
00434     
00435     ScaLAPACK1DMatrix_->FillComplete( false ) ; 
00436   }
00437   if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << " Amesos_Scalapack.cpp:417 debug_ = "
00438          << debug_ << std::endl;
00439   if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:402"
00440          << " nprow_ = " << nprow_
00441          << " npcol_ = " << npcol_ << std::endl ; 
00442   int info; 
00443   const int zero = 0 ; 
00444   if ( ictxt_ == -1313 ) {
00445     ictxt_ = 0 ; 
00446     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:408" << std::endl;
00447     SL_INIT_F77(&ictxt_, &nprow_, &npcol_) ; 
00448   }
00449   if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:410A" << std::endl;
00450   
00451   int nprow;
00452   int npcol;
00453   int myrow;
00454   int mycol;
00455   BLACS_GRIDINFO_F77(&ictxt_, &nprow, &npcol, &myrow, &mycol) ; 
00456   if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "iam_ = " << iam_ << " Amesos_Scalapack.cpp:410" << std::endl;
00457   if ( iam_ < nprow_ * npcol_ ) { 
00458     assert( nprow == nprow_ ) ; 
00459     if ( npcol != npcol_ ) std::cout << "Amesos_Scalapack.cpp:430 npcol = " << 
00460       npcol << " npcol_ = " << npcol_ << std::endl ; 
00461     assert( npcol == npcol_ ) ; 
00462     if ( TwoD_distribution_ ) {
00463       assert( myrow == myprow_ ) ; 
00464       assert( mycol == mypcol_ ) ; 
00465       lda_ = EPETRA_MAX(1,NumOurRows_) ;
00466     } else { 
00467       assert( myrow == 0 ) ; 
00468       assert( mycol == iam_ ) ; 
00469       nb_ = m_per_p_;
00470       lda_ = EPETRA_MAX(1,NumGlobalElements_);
00471     }
00472     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  
00473            << "Amesos_Scalapack.cpp: " << __LINE__ 
00474            << " TwoD_distribution_ = "  << TwoD_distribution_ 
00475            << " NumGlobalElements_ = "  << NumGlobalElements_ 
00476            << " debug_ = "  << debug_ 
00477            << " nb_ = "  << nb_ 
00478            << " lda_ = "  << lda_ 
00479            << " nprow_ = "  << nprow_ 
00480            << " npcol_ = "  << npcol_ 
00481            << " myprow_ = "  << myprow_ 
00482            << " mypcol_ = "  << mypcol_ 
00483            << " iam_ = "  << iam_ << std::endl ;
00484     AMESOS_PRINT( myprow_ );
00485     DESCINIT_F77(DescA_, 
00486      &NumGlobalElements_, 
00487      &NumGlobalElements_, 
00488      &nb_,
00489      &nb_,
00490      &zero,
00491      &zero,
00492      &ictxt_,
00493      &lda_,
00494      &info) ;
00495     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:441" << std::endl;
00496     assert( info == 0 ) ; 
00497   } else {
00498     DescA_[0] = -13;
00499     if (  debug_ == 1) std::cout  << "iam_ = " << iam_  << "Amesos_Scalapack.cpp:458 nprow = " << nprow << std::endl;
00500     assert( nprow == -1 ) ; 
00501   }
00502   
00503   if (  debug_ == 1) std::cout  << "Amesos_Scalapack.cpp:446" << std::endl;
00504   MatTime_ += Time_->ElapsedTime();
00505   
00506   return 0;
00507 }
00508 
00509 
00510 int Amesos_Scalapack::ConvertToScalapack(){
00511   
00512   //
00513   //  Convert matrix and vector to the form that Scalapack expects
00514   //  ScaLAPACK accepts the matrix to be in any 2D block-cyclic form
00515   //
00516   //  Amesos_ScaLAPACK uses one of two 2D data distributions: 
00517   //  a simple 1D non-cyclic data distribution with npcol= 1, or a 
00518   //  full  2D block-cyclic data distribution.
00519   //
00520   //  2D data distribvution:
00521   //    Because the Epetra export operation is oriented toward a 1D 
00522   //    data distribution in which each row is entirely stored on 
00523   //    a single process, we create two intermediate matrices: FatIn and
00524   //    FatOut, both of which have dimension:  
00525   //      NumGlobalElements * nprow by NumGlobalElements
00526   //    This allows each row of FatOut to be owned by a single process.
00527   //    The larger dimension does not significantly increase the 
00528   //    storage requirements and allows the export operation to be 
00529   //    efficient.  
00530   //
00531   //  1D data distribution:
00532   //  We have chosen the simplest 2D block-cyclic form, a 1D blocked (not-cyclic)
00533   //  data distribution, for the matrix A.
00534   //  We use the same distribution for the multivectors X and B.  However, 
00535   //  except for very large numbers of right hand sides, this places all of X and B
00536   //  on process 0, making it effectively a serial matrix.  
00537   //  
00538   //  For now, we simply treat X and B as serial matrices (as viewed from epetra)
00539   //  though ScaLAPACK treats them as distributed matrices. 
00540   //
00541   
00542   if( debug_ == 1 ) std::cout << "Entering `ConvertToScalapack()'" << std::endl;
00543   
00544   Time_->ResetStartTime();
00545   
00546   if ( iam_ < nprow_ * npcol_ ) { 
00547     if ( TwoD_distribution_ ) { 
00548       
00549       DenseA_.resize( NumOurRows_ * NumOurColumns_ ); 
00550       for ( int i = 0 ; i < (int)DenseA_.size() ; i++ ) DenseA_[i] = 0 ; 
00551       assert( lda_ == EPETRA_MAX(1,NumOurRows_) ) ;
00552       assert( DescA_[8] == lda_ ) ;
00553       
00554       int NzThisRow ;
00555       int MyRow;
00556       
00557       double *RowValues;
00558       int *ColIndices;
00559       int MaxNumEntries = FatOut_->MaxNumEntries();
00560       
00561       std::vector<int>ColIndicesV(MaxNumEntries);
00562       std::vector<double>RowValuesV(MaxNumEntries);
00563       
00564       int NumMyElements = FatOut_->NumMyRows() ; 
00565       for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
00566   EPETRA_CHK_ERR( FatOut_->
00567       ExtractMyRowView( MyRow, NzThisRow, RowValues, 
00568             ColIndices ) != 0 ) ;
00569   //
00570   //  The following eight lines are just a sanity check on MyRow:
00571   //
00572   int MyGlobalRow =  FatOut_->GRID( MyRow );
00573   assert( MyGlobalRow%npcol_ == mypcol_ ) ;   // I should only own rows belonging to my processor column
00574   int MyTrueRow = MyGlobalRow/npcol_ ;  // This is the original row
00575   int UniformRows =  ( MyTrueRow / ( nprow_ * nb_ ) ) * nb_ ; 
00576   int AllExcessRows = MyTrueRow - UniformRows * nprow_ ; 
00577   int OurExcessRows =  AllExcessRows - ( myprow_ * nb_ ) ; 
00578   
00579   if (  MyRow != UniformRows + OurExcessRows ) {
00580     std::cout << " iam _ = " << iam_ 
00581          << " MyGlobalRow = " << MyGlobalRow 
00582          << " MyTrueRow = " << MyTrueRow 
00583          << " UniformRows = " << UniformRows 
00584          << " AllExcessRows = " << AllExcessRows 
00585          << " OurExcessRows = " << OurExcessRows 
00586          << " MyRow = " << MyRow  << std::endl ;  
00587   }
00588   
00589   assert( OurExcessRows >= 0 &&  OurExcessRows < nb_ );
00590   assert( MyRow == UniformRows + OurExcessRows ) ; 
00591   
00592   for ( int j = 0; j < NzThisRow; j++ ) { 
00593     assert(  FatOut_->RowMatrixColMap().GID( ColIndices[j] ) ==
00594        FatOut_->GCID( ColIndices[j] ) );
00595     
00596     int MyGlobalCol =  FatOut_->GCID( ColIndices[j] );
00597     assert( (MyGlobalCol/nb_)%npcol_ == mypcol_ ) ; 
00598     int UniformCols =  ( MyGlobalCol / ( npcol_ * nb_ ) ) * nb_ ; 
00599     int AllExcessCols = MyGlobalCol - UniformCols * npcol_ ; 
00600     int OurExcessCols =  AllExcessCols - ( mypcol_ * nb_ ) ; 
00601     assert( OurExcessCols >= 0 &&  OurExcessCols < nb_ );
00602     int MyCol = UniformCols + OurExcessCols ; 
00603     
00604     DenseA_[ MyCol * lda_ + MyRow ] = RowValues[j] ; 
00605   }
00606       }
00607       
00608     } else { 
00609       
00610       int NumMyElements = ScaLAPACK1DMatrix_->NumMyRows() ; 
00611       
00612       assert( NumGlobalElements_ ==ScaLAPACK1DMatrix_->NumGlobalRows());
00613       assert( NumGlobalElements_ ==ScaLAPACK1DMatrix_->NumGlobalCols());
00614       DenseA_.resize( NumGlobalElements_ * NumMyElements ) ;
00615       for ( int i = 0 ; i < (int)DenseA_.size() ; i++ ) DenseA_[i] = 0 ; 
00616       
00617       int NzThisRow ;
00618       int MyRow;
00619       
00620       double *RowValues;
00621       int *ColIndices;
00622       int MaxNumEntries = ScaLAPACK1DMatrix_->MaxNumEntries();
00623       
00624       assert( DescA_[8] == lda_ ) ; //  Double check Lda
00625       std::vector<int>ColIndicesV(MaxNumEntries);
00626       std::vector<double>RowValuesV(MaxNumEntries);
00627       
00628       for ( MyRow = 0; MyRow < NumMyElements ; MyRow++ ) {
00629   EPETRA_CHK_ERR( ScaLAPACK1DMatrix_->
00630       ExtractMyRowView( MyRow, NzThisRow, RowValues, 
00631             ColIndices ) != 0 ) ;
00632   
00633   for ( int j = 0; j < NzThisRow; j++ ) { 
00634     DenseA_[ ( ScaLAPACK1DMatrix_->RowMatrixColMap().GID( ColIndices[j] ) ) 
00635        + MyRow * NumGlobalElements_ ] = RowValues[j] ; 
00636   }
00637       }
00638       //
00639       //  Create a map to allow us to redistribute the vectors X and B 
00640       //
00641       Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
00642       const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
00643       assert( NumGlobalElements_ == OriginalMap.NumGlobalElements() ) ;
00644       int NumMyElements_ = 0 ;
00645       if (iam_==0) NumMyElements_ = NumGlobalElements_;
00646       
00647       if (VectorMap_) { delete VectorMap_ ; VectorMap_ = 0 ; } 
00648       VectorMap_ = new Epetra_Map( NumGlobalElements_, NumMyElements_, 0, Comm() );
00649     }
00650   }
00651   ConTime_ += Time_->ElapsedTime();
00652   
00653   return 0;
00654 }   
00655 
00656 
00657 int Amesos_Scalapack::SetParameters( Teuchos::ParameterList &ParameterList ) {
00658   
00659   if( debug_ == 1 ) std::cout << "Entering `SetParameters()'" << std::endl;
00660   
00661   // retrive general parameters
00662   SetStatusParameters( ParameterList );
00663   
00664   SetControlParameters( ParameterList );
00665   
00666   //
00667   //  We have to set these to their defaults here because user codes 
00668   //  are not guaranteed to have a "Scalapack" parameter list.
00669   //
00670   TwoD_distribution_ = true; 
00671   grid_nb_ = 32; 
00672   
00673   //  Some compilers reject the following cast:
00674   //  if( &ParameterList == 0 ) return 0;
00675   
00676   // ========================================= //
00677   // retrive ScaLAPACK's parameters from list. //
00678   // ========================================= //
00679   
00680   // retrive general parameters
00681   //  check to see if they exist before trying to retrieve them
00682   if (ParameterList.isSublist("Scalapack") ) {
00683     const Teuchos::ParameterList ScalapackParams = ParameterList.sublist("Scalapack") ;
00684     //  Fix Bug #3251 
00685     if ( ScalapackParams.isParameter("2D distribution") )
00686       TwoD_distribution_ = ScalapackParams.get<bool>("2D distribution");
00687     if ( ScalapackParams.isParameter("grid_nb") )
00688       grid_nb_ = ScalapackParams.get<int>("grid_nb");
00689   }  
00690   
00691   return 0;
00692 }
00693 
00694 int Amesos_Scalapack::PerformNumericFactorization( ) {
00695   
00696   if( debug_ == 1 ) std::cout << "Entering `PerformNumericFactorization()'" << std::endl;
00697   
00698   Time_->ResetStartTime();  
00699   
00700   Ipiv_.resize(NumGlobalElements_) ;
00701   for (int i=0; i <NumGlobalElements_ ; i++) 
00702     Ipiv_[i] = 0 ; // kludge - just to see if this makes valgrind happy
00703   
00704   if ( false) std::cout  << " Amesos_Scalapack.cpp: 711 iam_ = " << iam_ << " DescA = " 
00705         << DescA_[0] << " " 
00706         << DescA_[1] << " " 
00707         << DescA_[2] << " " 
00708         << DescA_[3] << " " 
00709         << DescA_[4] << " " 
00710         << DescA_[5] << " " 
00711         << DescA_[6] << " " 
00712         << DescA_[7] << " " 
00713         << DescA_[8] << " " 
00714         << std::endl ; 
00715   
00716 #if 1
00717   if( NumGlobalElements_ < 10 && nprow_ == 1 && npcol_ == 1 && debug_ == 1 ) {
00718     assert( lda_ == NumGlobalElements_ ) ; 
00719     std::cout << " DenseA = " << std::endl ; 
00720     for (int i=0 ; i < NumGlobalElements_; i++ ) {
00721       for (int j=0 ; j < NumGlobalElements_; j++ ) {
00722   if ( DenseA_[ i+j*lda_ ] < 0 ) {
00723     DenseA_[ i+j*lda_ ] *= (1+1e-5) ;   // kludge fixme debugxx - just to let vaglrind check to be sure that DenseA is initialized
00724   }
00725   //  std::cout << DenseA_[ i+j*lda_ ] << "\t"; 
00726       }
00727       //      std::cout << std::endl ; 
00728     }
00729   }
00730 #endif
00731   
00732   int Ierr[1] ; 
00733   Ierr[0] = 0 ; 
00734   const int one = 1 ; 
00735   if ( iam_ < nprow_ * npcol_ ) {
00736     if ( nprow_ * npcol_ == 1 ) { 
00737       DGETRF_F77(&NumGlobalElements_,  
00738      &NumGlobalElements_, 
00739      &DenseA_[0],
00740      &lda_,
00741      &Ipiv_[0],
00742      Ierr) ;
00743     } else { 
00744       PDGETRF_F77(&NumGlobalElements_,  
00745       &NumGlobalElements_, 
00746       &DenseA_[0],
00747       &one,
00748       &one, 
00749       DescA_,
00750       &Ipiv_[0],
00751       Ierr) ;
00752     }
00753   }
00754   
00755   if ( debug_ == 1  && Ierr[0] != 0 ) 
00756     std::cout << " Amesos_Scalapack.cpp:738 iam_ = " << iam_ 
00757    << " Ierr = " << Ierr[0]  << std::endl ; 
00758   
00759   //  All processes should return the same error code
00760   if ( nprow_ * npcol_ < Comm().NumProc() ) 
00761     Comm().Broadcast( Ierr, 1, 0 ) ; 
00762   
00763   NumTime_ += Time_->ElapsedTime();
00764   
00765   return Ierr[0];
00766 }
00767 
00768 
00769 bool Amesos_Scalapack::MatrixShapeOK() const { 
00770   bool OK ;
00771   
00772   if ( GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() != 
00773        GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) OK = false;
00774   return OK; 
00775 }
00776 
00777 
00778 int Amesos_Scalapack::SymbolicFactorization() {
00779   
00780   if( debug_ == 1 ) std::cout << "Entering `PerformSymbolicFactorization()'" << std::endl;
00781   
00782   if( Time_ == 0 ) Time_ = new Epetra_Time( Comm() );
00783   
00784   NumSymbolicFact_++;
00785   
00786   return 0;
00787 }
00788 
00789 int Amesos_Scalapack::NumericFactorization() {
00790   
00791   if( debug_ == 1 ) std::cout << "Entering `NumericFactorization()'" << std::endl;
00792   
00793   NumNumericFact_++;
00794   
00795   iam_ = Comm().MyPID();
00796   
00797   Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
00798   const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
00799   NumGlobalElements_ = OriginalMap.NumGlobalElements();
00800 
00801   NumGlobalNonzeros_ = RowMatrixA->NumGlobalNonzeros();
00802   
00803   RedistributeA();
00804   ConvertToScalapack();
00805   
00806   return PerformNumericFactorization( );
00807 }
00808 
00809 
00810 int Amesos_Scalapack::Solve() { 
00811   
00812   if( debug_ == 1 ) std::cout << "Entering `Solve()'" << std::endl;
00813   
00814   NumSolve_++;
00815   
00816   Epetra_MultiVector   *vecX = Problem_->GetLHS() ; 
00817   Epetra_MultiVector   *vecB = Problem_->GetRHS() ; 
00818   
00819   //
00820   //  Compute the number of right hands sides 
00821   //  (and check that X and B have the same shape) 
00822   //
00823   int nrhs; 
00824   if ( vecX == 0 ) { 
00825     nrhs = 0 ;
00826     EPETRA_CHK_ERR( vecB != 0 ) ; 
00827   } else { 
00828     nrhs = vecX->NumVectors() ; 
00829     EPETRA_CHK_ERR( vecB->NumVectors() != nrhs ) ; 
00830   }
00831   
00832   Epetra_MultiVector *ScalapackB =0;
00833   Epetra_MultiVector *ScalapackX =0;
00834   //
00835   //  Extract Scalapack versions of X and B 
00836   //
00837   double *ScalapackXvalues ;
00838   
00839   Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
00840   Time_->ResetStartTime(); // track time to broadcast vectors
00841   //
00842   //  Copy B to the scalapack version of B
00843   //
00844   const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap();
00845   Epetra_MultiVector *ScalapackXextract = new Epetra_MultiVector( *VectorMap_, nrhs ) ; 
00846   Epetra_MultiVector *ScalapackBextract = new Epetra_MultiVector( *VectorMap_, nrhs ) ; 
00847   
00848   Epetra_Import ImportToScalapack( *VectorMap_, OriginalMap );
00849   ScalapackBextract->Import( *vecB, ImportToScalapack, Insert ) ;
00850   ScalapackB = ScalapackBextract ; 
00851   ScalapackX = ScalapackXextract ; 
00852   
00853   VecTime_ += Time_->ElapsedTime();
00854   
00855   //
00856   //  Call SCALAPACKs PDGETRS to perform the solve
00857   //
00858   
00859   int DescX[10];  
00860   
00861   ScalapackX->Scale(1.0, *ScalapackB) ;  
00862   
00863   int ScalapackXlda ; 
00864   
00865   Time_->ResetStartTime(); // tract time to solve
00866   
00867   //
00868   //  Setup DescX 
00869   //
00870   
00871   if( nrhs > nb_ ) {
00872     EPETRA_CHK_ERR( -2 );  
00873   }
00874   
00875   int Ierr[1] ; 
00876   Ierr[0] = 0 ; 
00877   const int zero = 0 ; 
00878   const int one = 1 ; 
00879   if ( iam_ < nprow_ * npcol_ ) {
00880     assert( ScalapackX->ExtractView( &ScalapackXvalues, &ScalapackXlda ) == 0 ) ; 
00881     
00882     if ( false ) std::cout << "Amesos_Scalapack.cpp: " << __LINE__ << " ScalapackXlda = "  <<  ScalapackXlda 
00883           << " lda_ = "  << lda_ 
00884           << " nprow_ = "  << nprow_ 
00885           << " npcol_ = "  << npcol_ 
00886           << " myprow_ = "  << myprow_ 
00887           << " mypcol_ = "  << mypcol_ 
00888           << " iam_ = "  << iam_ << std::endl ;
00889     if (  TwoD_distribution_ )    assert( mypcol_ >0 || EPETRA_MAX(ScalapackXlda,1) == lda_ ) ; 
00890     
00891     DESCINIT_F77(DescX, 
00892      &NumGlobalElements_, 
00893      &nrhs, 
00894      &nb_,
00895      &nb_,
00896      &zero,
00897      &zero,
00898      &ictxt_,
00899      &lda_,
00900      Ierr ) ;
00901     assert( Ierr[0] == 0 ) ; 
00902     
00903     //
00904     //  For the 1D data distribution, we factor the transposed 
00905     //  matrix, hence we must invert the sense of the transposition
00906     //
00907     char trans = 'N';
00908     if ( TwoD_distribution_ ) {
00909       if ( UseTranspose() ) trans = 'T' ;
00910     } else {
00911       
00912       if ( ! UseTranspose() ) trans = 'T' ;
00913     }
00914     
00915     if ( nprow_ * npcol_ == 1 ) { 
00916       DGETRS_F77(&trans,
00917      &NumGlobalElements_,  
00918      &nrhs, 
00919      &DenseA_[0],
00920      &lda_,
00921      &Ipiv_[0],
00922      ScalapackXvalues,
00923      &lda_,
00924      Ierr ) ;
00925     } else { 
00926       PDGETRS_F77(&trans,
00927       &NumGlobalElements_,  
00928       &nrhs, 
00929       &DenseA_[0],
00930       &one,
00931       &one, 
00932       DescA_,
00933       &Ipiv_[0],
00934       ScalapackXvalues,
00935       &one,
00936       &one, 
00937       DescX,
00938       Ierr ) ;
00939     }
00940   }
00941   
00942   SolTime_ += Time_->ElapsedTime();
00943   
00944   Time_->ResetStartTime();  // track time to broadcast vectors
00945   //
00946   //  Copy X back to the original vector
00947   // 
00948   Epetra_Import ImportFromScalapack( OriginalMap, *VectorMap_ );
00949   vecX->Import( *ScalapackX, ImportFromScalapack, Insert ) ;
00950   delete ScalapackBextract ;
00951   delete ScalapackXextract ;
00952   
00953   VecTime_ += Time_->ElapsedTime();
00954   
00955   //  All processes should return the same error code
00956   if ( nprow_ * npcol_ < Comm().NumProc() ) 
00957     Comm().Broadcast( Ierr, 1, 0 ) ; 
00958   
00959   // MS // compute vector norms
00960   if( ComputeVectorNorms_ == true || verbose_ == 2 ) {
00961     double NormLHS, NormRHS;
00962     for( int i=0 ; i<nrhs ; ++i ) {
00963       assert((*vecX)(i)->Norm2(&NormLHS)==0);
00964       assert((*vecB)(i)->Norm2(&NormRHS)==0);
00965       if( verbose_ && Comm().MyPID() == 0 ) {
00966   std::cout << "Amesos_Scalapack : vector " << i << ", ||x|| = " << NormLHS
00967        << ", ||b|| = " << NormRHS << std::endl;
00968       }
00969     }
00970   }
00971   
00972   // MS // compute true residual
00973   if( ComputeTrueResidual_ == true || verbose_ == 2  ) {
00974     double Norm;
00975     Epetra_MultiVector Ax(vecB->Map(),nrhs);
00976     for( int i=0 ; i<nrhs ; ++i ) {
00977       (Problem_->GetMatrix()->Multiply(UseTranspose(), *((*vecX)(i)), Ax));
00978       (Ax.Update(1.0, *((*vecB)(i)), -1.0));
00979       (Ax.Norm2(&Norm));
00980       
00981       if( verbose_ && Comm().MyPID() == 0 ) {
00982   std::cout << "Amesos_Scalapack : vector " << i << ", ||Ax - b|| = " << Norm << std::endl;
00983       }
00984     }
00985   }
00986   
00987   return Ierr[0];
00988   
00989 }
00990 
00991 
00992 // ================================================ ====== ==== ==== == =
00993 
00994 void Amesos_Scalapack::PrintStatus() const
00995 {
00996   
00997   if( iam_ != 0  ) return;
00998   
00999   std::cout << "----------------------------------------------------------------------------" << std::endl;
01000   std::cout << "Amesos_Scalapack : Matrix has " << NumGlobalElements_ << " rows"
01001        << " and " << NumGlobalNonzeros_ << " nonzeros" << std::endl;
01002   std::cout << "Amesos_Scalapack : Nonzero elements per row = "
01003        << 1.0*NumGlobalNonzeros_/NumGlobalElements_ << std::endl;
01004   std::cout << "Amesos_Scalapack : Percentage of nonzero elements = "
01005        << 100.0*NumGlobalNonzeros_/(pow(NumGlobalElements_,2.0)) << std::endl;
01006   std::cout << "Amesos_Scalapack : Use transpose = " << UseTranspose_ << std::endl;
01007   std::cout << "----------------------------------------------------------------------------" << std::endl;
01008   
01009   return;
01010   
01011 }
01012 
01013 // ================================================ ====== ==== ==== == =
01014 
01015 void Amesos_Scalapack::PrintTiming() const
01016 {
01017   if( iam_ ) return;
01018   
01019   std::cout << "----------------------------------------------------------------------------" << std::endl;
01020   std::cout << "Amesos_Scalapack : Time to convert matrix to ScaLAPACK format = "
01021        << ConTime_ << " (s)" << std::endl;
01022   std::cout << "Amesos_Scalapack : Time to redistribute matrix = "
01023        << MatTime_ << " (s)" << std::endl;
01024   std::cout << "Amesos_Scalapack : Time to redistribute vectors = "
01025        << VecTime_ << " (s)" << std::endl;
01026   std::cout << "Amesos_Scalapack : Number of symbolic factorizations = "
01027        << NumSymbolicFact_ << std::endl;
01028   std::cout << "Amesos_Scalapack : Time for sym fact = "
01029        << SymTime_ << " (s), avg = " << SymTime_/NumSymbolicFact_
01030        << " (s)" << std::endl;
01031   std::cout << "Amesos_Scalapack : Number of numeric factorizations = "
01032        << NumNumericFact_ << std::endl;
01033   std::cout << "Amesos_Scalapack : Time for num fact = "
01034        << NumTime_ << " (s), avg = " << NumTime_/NumNumericFact_
01035        << " (s)" << std::endl;
01036   std::cout << "Amesos_Scalapack : Number of solve phases = "
01037        << NumSolve_ << std::endl;
01038   std::cout << "Amesos_Scalapack : Time for solve = "
01039        << SolTime_ << " (s), avg = " << SolTime_/NumSolve_
01040        << " (s)" << std::endl;
01041   std::cout << "----------------------------------------------------------------------------" << std::endl;
01042   
01043   return;
01044 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines