Amesos Package Browser (Single Doxygen Collection) Development
Amesos_Dscpack.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_Dscpack.h"
00030 #include "Epetra_Map.h"
00031 #include "Epetra_Import.h"
00032 #include "Epetra_CrsMatrix.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_Util.h"
00035 #include "Teuchos_RCP.hpp"
00036 #include <algorithm>
00037 
00038 #define DBL_R_NUM
00039 extern "C" {
00040 #include "dscmain.h"
00041 }
00042 
00043 class Amesos_Dscpack_Pimpl {
00044 public:
00045   DSC_Solver  MyDSCObject_;
00046 } ;
00047 //=============================================================================
00048 Amesos_Dscpack::Amesos_Dscpack(const Epetra_LinearProblem &prob ) : 
00049   PrivateDscpackData_( rcp( new Amesos_Dscpack_Pimpl() ) ),
00050   DscNumProcs(-1), // will be set later
00051   MaxProcs_(-1),
00052   MtxRedistTime_(-1),
00053   MtxConvTime_(-1),
00054   VecRedistTime_(-1),
00055   SymFactTime_(-1),
00056   NumFactTime_(-1),
00057   SolveTime_(-1),
00058   OverheadTime_(-1)
00059 {  
00060   Problem_ = &prob ; 
00061   A_and_LU_built = false ; 
00062   FirstCallToSolve_ = true ; 
00063   PrivateDscpackData_->MyDSCObject_ = DSC_Begin() ; 
00064 
00065   MyDscRank = -1 ; 
00066 }
00067 
00068 //=============================================================================
00069 Amesos_Dscpack::~Amesos_Dscpack(void) 
00070 {
00071   // MS // print out some information if required by the user
00072   if( (verbose_ && PrintTiming_) || verbose_ == 2 ) PrintTiming();
00073   if( (verbose_ && PrintStatus_) || verbose_ == 2 ) PrintStatus();
00074 
00075   if ( MyDscRank>=0 && A_and_LU_built ) { 
00076     DSC_FreeAll( PrivateDscpackData_->MyDSCObject_ ) ; 
00077     DSC_Close0( PrivateDscpackData_->MyDSCObject_ ) ; 
00078   }
00079   DSC_End( PrivateDscpackData_->MyDSCObject_ ) ; 
00080 }
00081 
00082 //=============================================================================
00083 int Amesos_Dscpack::SetParameters( Teuchos::ParameterList &ParameterList) 
00084 {
00085   // ========================================= //
00086   // retrive DSCPACK's parameters from list.   //
00087   // default values defined in the constructor //
00088   // ========================================= //
00089 
00090   // retrive general parameters
00091 
00092   SetStatusParameters( ParameterList );
00093 
00094   SetControlParameters( ParameterList );
00095   
00096   // MS // NO DSCPACK-specify parameters at this point, uncomment
00097   // MS // as necessary
00098   /*
00099   if (ParameterList.isSublist("Dscpack") ) {
00100     Teuchos::ParameterList DscpackParams = ParameterList.sublist("Dscpack") ;
00101   }
00102   */
00103   
00104   return 0;
00105 }
00106 
00107 //=============================================================================
00108 int Amesos_Dscpack::PerformSymbolicFactorization()
00109 {
00110   ResetTimer(0);
00111   ResetTimer(1);
00112 
00113   MyPID_    = Comm().MyPID();
00114   NumProcs_ = Comm().NumProc();
00115   
00116   Epetra_RowMatrix *RowMatrixA = Problem_->GetMatrix();
00117   if (RowMatrixA == 0)
00118     AMESOS_CHK_ERR(-1);
00119 
00120   const Epetra_Map& OriginalMap = RowMatrixA->RowMatrixRowMap() ;
00121   const Epetra_MpiComm& comm1   = dynamic_cast<const Epetra_MpiComm &> (Comm());
00122   int numrows                   = RowMatrixA->NumGlobalRows();
00123   int numentries                = RowMatrixA->NumGlobalNonzeros();
00124 
00125   Teuchos::RCP<Epetra_CrsGraph> Graph;
00126 
00127   Epetra_CrsMatrix* CastCrsMatrixA = 
00128     dynamic_cast<Epetra_CrsMatrix*>(RowMatrixA); 
00129 
00130   if (CastCrsMatrixA)
00131   {
00132     Graph = Teuchos::rcp(const_cast<Epetra_CrsGraph*>(&(CastCrsMatrixA->Graph())), false);
00133   }
00134   else
00135   {
00136     int MaxNumEntries = RowMatrixA->MaxNumEntries();
00137     Graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, OriginalMap, MaxNumEntries));
00138 
00139     std::vector<int>    Indices(MaxNumEntries);
00140     std::vector<double> Values(MaxNumEntries);
00141 
00142     for (int i = 0 ; i < RowMatrixA->NumMyRows() ; ++i)
00143     {
00144       int NumEntries;
00145       RowMatrixA->ExtractMyRowCopy(i, MaxNumEntries, NumEntries,
00146                                    &Values[0], &Indices[0]);
00147 
00148       for (int j = 0 ; j < NumEntries ; ++j)
00149         Indices[j] = RowMatrixA->RowMatrixColMap().GID(Indices[j]);
00150 
00151       int GlobalRow = RowMatrixA->RowMatrixRowMap().GID(i);
00152       Graph->InsertGlobalIndices(GlobalRow, NumEntries, &Indices[0]);
00153     }
00154 
00155     Graph->FillComplete();
00156   }
00157 
00158   //
00159   //  Create a replicated map and graph 
00160   //
00161   std::vector<int> AllIDs( numrows ) ; 
00162   for ( int i = 0; i < numrows ; i++ ) AllIDs[i] = i ; 
00163 
00164   Epetra_Map      ReplicatedMap( -1, numrows, &AllIDs[0], 0, Comm());
00165   Epetra_Import   ReplicatedImporter(ReplicatedMap, OriginalMap);
00166   Epetra_CrsGraph ReplicatedGraph( Copy, ReplicatedMap, 0 ); 
00167 
00168   AMESOS_CHK_ERR(ReplicatedGraph.Import(*Graph, ReplicatedImporter, Insert));
00169   AMESOS_CHK_ERR(ReplicatedGraph.FillComplete());
00170 
00171   //
00172   //  Convert the matrix to Ap, Ai
00173   //
00174   std::vector <int> Replicates(numrows);
00175   std::vector <int> Ap(numrows + 1);
00176   std::vector <int> Ai(EPETRA_MAX(numrows, numentries));
00177 
00178   for( int i = 0 ; i < numrows; i++ ) Replicates[i] = 1; 
00179   
00180   int NumEntriesPerRow ;
00181   int *ColIndices = 0 ;
00182   int Ai_index = 0 ; 
00183   for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
00184     AMESOS_CHK_ERR( ReplicatedGraph.ExtractMyRowView( MyRow, NumEntriesPerRow, ColIndices ) );
00185     Ap[MyRow] = Ai_index ; 
00186     for ( int j = 0; j < NumEntriesPerRow; j++ ) { 
00187       Ai[Ai_index] = ColIndices[j] ; 
00188       Ai_index++;
00189     }
00190   }
00191   assert( Ai_index == numentries ) ; 
00192   Ap[ numrows ] = Ai_index ; 
00193   
00194   MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);
00195 
00196   ResetTimer(0);
00197 
00198   //
00199   //  Call Dscpack Symbolic Factorization
00200   //  
00201   int OrderCode = 2;
00202   std::vector<double> MyANonZ;
00203   
00204   NumLocalNonz = 0 ; 
00205   GlobalStructNewColNum = 0 ; 
00206   GlobalStructNewNum = 0 ;  
00207   GlobalStructOwner = 0 ; 
00208   LocalStructOldNum = 0 ; 
00209   
00210   NumGlobalCols = 0 ; 
00211   
00212   // MS // Have to define the maximum number of processes to be used
00213   // MS // This is only a suggestion as Dscpack uses a number of processes that is a power of 2  
00214 
00215   int NumGlobalNonzeros = GetProblem()->GetMatrix()->NumGlobalNonzeros();
00216   int NumRows = GetProblem()->GetMatrix()->NumGlobalRows(); 
00217 
00218   // optimal value for MaxProcs == -1
00219   
00220   int OptNumProcs1 = 1+EPETRA_MAX( NumRows/10000, NumGlobalNonzeros/1000000 );
00221   OptNumProcs1 = EPETRA_MIN(NumProcs_,OptNumProcs1 );
00222 
00223   // optimal value for MaxProcs == -2
00224 
00225   int OptNumProcs2 = (int)sqrt(1.0 * NumProcs_);
00226   if( OptNumProcs2 < 1 ) OptNumProcs2 = 1;
00227 
00228   // fix the value of MaxProcs
00229 
00230   switch (MaxProcs_) 
00231   {
00232   case -1:
00233     MaxProcs_ = EPETRA_MIN(OptNumProcs1, NumProcs_);
00234     break;
00235   case -2:
00236     MaxProcs_ = EPETRA_MIN(OptNumProcs2, NumProcs_);
00237     break;
00238   case -3:
00239     MaxProcs_ = NumProcs_;
00240     break;
00241   }
00242 
00243 #if 0
00244   if (MyDscRank>=0 && A_and_LU_built) { 
00245     DSC_ReFactorInitialize(PrivateDscpackData_->MyDSCObject);
00246   }
00247 #endif
00248   //  if ( ! A_and_LU_built ) { 
00249   //    DSC_End( PrivateDscpackData_->MyDSCObject ) ; 
00250   //    PrivateDscpackData_->MyDSCObject = DSC_Begin() ;
00251   //  } 
00252 
00253   // MS // here I continue with the old code...
00254   
00255   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00256 
00257   DscNumProcs = 1 ; 
00258   int DscMax = DSC_Analyze( numrows, &Ap[0], &Ai[0], &Replicates[0] );
00259 
00260   while ( DscNumProcs * 2 <=EPETRA_MIN( MaxProcs_, DscMax ) )  DscNumProcs *= 2 ;
00261   
00262   MyDscRank = -1; 
00263   DSC_Open0( PrivateDscpackData_->MyDSCObject_, DscNumProcs, &MyDscRank, comm1.Comm()) ; 
00264   
00265   NumLocalCols = 0 ; // This is for those processes not in the Dsc grid
00266   if ( MyDscRank >= 0 ) { 
00267     assert( MyPID_ == MyDscRank ) ; 
00268     AMESOS_CHK_ERR( DSC_Order ( PrivateDscpackData_->MyDSCObject_, OrderCode, numrows, &Ap[0], &Ai[0], 
00269         &Replicates[0], &NumGlobalCols, &NumLocalStructs, 
00270         &NumLocalCols, &NumLocalNonz, 
00271         &GlobalStructNewColNum, &GlobalStructNewNum, 
00272         &GlobalStructOwner, &LocalStructOldNum ) ) ; 
00273     assert( NumGlobalCols == numrows ) ; 
00274     assert( NumLocalCols == NumLocalStructs ) ; 
00275   }
00276 
00277   if ( MyDscRank >= 0 ) { 
00278     int MaxSingleBlock; 
00279     
00280     const int Limit = 5000000 ;  //  Memory Limit set to 5 Terabytes 
00281     AMESOS_CHK_ERR( DSC_SFactor ( PrivateDscpackData_->MyDSCObject_, &TotalMemory_, 
00282           &MaxSingleBlock, Limit, DSC_LBLAS3, DSC_DBLAS2 ) ) ; 
00283     
00284   }
00285   
00286   //  A_and_LU_built = true;   // If you uncomment this, TestOptions fails
00287   
00288   SymFactTime_ = AddTime("Total symbolic factorization time", SymFactTime_, 0);
00289 
00290   return(0);
00291 }
00292 
00293 //=============================================================================
00294 int Amesos_Dscpack::PerformNumericFactorization()
00295 {
00296   ResetTimer(0);
00297   ResetTimer(1);
00298 
00299   Epetra_RowMatrix* RowMatrixA = Problem_->GetMatrix();
00300   if (RowMatrixA == 0)
00301     AMESOS_CHK_ERR(-1);
00302 
00303   const Epetra_Map& OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
00304 
00305   int numrows = RowMatrixA->NumGlobalRows();
00306   assert( numrows == RowMatrixA->NumGlobalCols() );
00307   
00308   //
00309   //  Call Dscpack to perform Numeric Factorization
00310   //  
00311   std::vector<double> MyANonZ;
00312 #if 0
00313     if ( IsNumericFactorizationOK_ ) { 
00314       DSC_ReFactorInitialize(PrivateDscpackData_->MyDSCObject);
00315     }
00316 #endif
00317 
00318   DscRowMap_ = Teuchos::rcp(new Epetra_Map(numrows, NumLocalCols, 
00319                                            LocalStructOldNum, 0, Comm()));
00320 
00321   if (DscRowMap_.get() == 0)
00322     AMESOS_CHK_ERR(-1);
00323   
00324   Importer_ = rcp(new Epetra_Import(DscRowMap(), OriginalMap));
00325     
00326   //
00327   //  Import from the CrsMatrix
00328   //
00329   Epetra_CrsMatrix DscMat(Copy, DscRowMap(), 0);
00330   AMESOS_CHK_ERR(DscMat.Import(*RowMatrixA, Importer(), Insert));
00331   AMESOS_CHK_ERR(DscMat.FillComplete()); 
00332 
00333   DscColMap_ = Teuchos::rcp(new Epetra_Map(DscMat.RowMatrixColMap()));
00334 
00335   assert( MyDscRank >= 0 || NumLocalNonz == 0 ) ;
00336   assert( MyDscRank >= 0 || NumLocalCols == 0 ) ;
00337   assert( MyDscRank >= 0 || NumGlobalCols == 0  ) ; 
00338   MyANonZ.resize( NumLocalNonz ) ; 
00339   int NonZIndex = 0 ;
00340 
00341   int max_num_entries = DscMat.MaxNumEntries() ; 
00342   std::vector<int> col_indices( max_num_entries ) ; 
00343   std::vector<double> mat_values( max_num_entries ) ; 
00344   assert( NumLocalCols == DscRowMap().NumMyElements() ) ;
00345   std::vector<int> my_global_elements( NumLocalCols ) ; 
00346   AMESOS_CHK_ERR(DscRowMap().MyGlobalElements( &my_global_elements[0] ) ) ;
00347   
00348   std::vector<int> GlobalStructOldColNum( NumGlobalCols ) ; 
00349   
00350   typedef std::pair<int, double> Data; 
00351   std::vector<Data> sort_array(max_num_entries); 
00352   std::vector<int>  sort_indices(max_num_entries);
00353   
00354   for ( int i = 0; i < NumLocalCols ; i++ ) { 
00355     assert( my_global_elements[i] == LocalStructOldNum[i] ) ; 
00356     int num_entries_this_row; 
00357 #ifdef USE_LOCAL
00358     AMESOS_CHK_ERR( DscMat.ExtractMyRowCopy( i, max_num_entries, num_entries_this_row, 
00359                &mat_values[0], &col_indices[0] ) ) ; 
00360 #else
00361     AMESOS_CHK_ERR( DscMat.ExtractGlobalRowCopy( DscMat.GRID(i), max_num_entries, num_entries_this_row, 
00362              &mat_values[0], &col_indices[0] ) ) ; 
00363 #endif
00364     int OldRowNumber =  LocalStructOldNum[i] ;
00365     if (GlobalStructOwner[ OldRowNumber ] == -1)
00366       AMESOS_CHK_ERR(-1);
00367     
00368     int NewRowNumber = GlobalStructNewColNum[ my_global_elements[ i ] ] ; 
00369     
00370     //
00371     //  Sort the column elements 
00372     //
00373     
00374     for ( int j = 0; j < num_entries_this_row; j++ ) { 
00375 #ifdef USE_LOCAL
00376       sort_array[j].first = GlobalStructNewColNum[ DscMat.GCID( col_indices[j])] ; 
00377       sort_indices[j] =  GlobalStructNewColNum[ DscMat.GCID( col_indices[j])] ; 
00378 #else
00379       sort_array[j].first = GlobalStructNewColNum[ col_indices[j] ]; 
00380 #endif
00381       sort_array[j].second = mat_values[j] ; 
00382     }
00383     sort(&sort_array[0], &sort_array[num_entries_this_row]);
00384     
00385     for ( int j = 0; j < num_entries_this_row; j++ ) { 
00386       int NewColNumber = sort_array[j].first ; 
00387       if ( NewRowNumber <= NewColNumber ) MyANonZ[ NonZIndex++ ] = sort_array[j].second ; 
00388 
00389 #ifndef USE_LOCAL
00390       assert( NonZIndex <= NumLocalNonz );  // This assert can fail on non-symmetric matrices
00391 #endif
00392     }
00393   }
00394   
00395   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00396 
00397   if ( MyDscRank >= 0 ) { 
00398     const int SchemeCode = 1; 
00399 #ifndef USE_LOCAL
00400     assert( NonZIndex == NumLocalNonz );
00401 #endif
00402     
00403     AMESOS_CHK_ERR( DSC_NFactor ( PrivateDscpackData_->MyDSCObject_, SchemeCode, &MyANonZ[0], 
00404           DSC_LLT,  DSC_LBLAS3, DSC_DBLAS2 ) ) ;
00405     
00406   }        //     if ( MyDscRank >= 0 ) 
00407   
00408   IsNumericFactorizationOK_ = true ; 
00409 
00410   NumFactTime_ = AddTime("Total numeric factorization time", NumFactTime_, 0);
00411   
00412   return(0);
00413 }
00414 
00415 bool Amesos_Dscpack::MatrixShapeOK() const 
00416 {
00417   bool OK =  GetProblem()->IsOperatorSymmetric() ;
00418 
00419   //
00420   //  The following test is redundant.  I have left it here in case the 
00421   //  IsOperatorSymmetric test turns out not to be reliable.
00422   //
00423   if ( GetProblem()->GetOperator()->OperatorRangeMap().NumGlobalPoints() != 
00424        GetProblem()->GetOperator()->OperatorDomainMap().NumGlobalPoints() ) OK = false;
00425   return OK; 
00426 }
00427 
00428 //=============================================================================
00429 int Amesos_Dscpack::SymbolicFactorization()
00430 {
00431   IsSymbolicFactorizationOK_ = false;
00432   IsNumericFactorizationOK_ = false;
00433   
00434   CreateTimer(Comm(), 2);
00435 
00436   AMESOS_CHK_ERR(PerformSymbolicFactorization());
00437 
00438   IsSymbolicFactorizationOK_ = true; 
00439   NumSymbolicFact_++;
00440 
00441   return(0);
00442 }
00443 
00444 //=============================================================================
00445 int Amesos_Dscpack::NumericFactorization()
00446 {
00447   IsNumericFactorizationOK_ = false;
00448 
00449   if (!IsSymbolicFactorizationOK_) 
00450     AMESOS_CHK_ERR(SymbolicFactorization());
00451 
00452   AMESOS_CHK_ERR(PerformNumericFactorization());
00453 
00454   IsNumericFactorizationOK_ = true;
00455   NumNumericFact_++;
00456   
00457   return(0);
00458 }
00459 
00460 //=============================================================================
00461 int Amesos_Dscpack::Solve()
00462 {
00463   if (IsNumericFactorizationOK_ == false) 
00464     AMESOS_CHK_ERR(NumericFactorization());
00465 
00466   ResetTimer(0);
00467   ResetTimer(1);
00468   
00469   Epetra_RowMatrix *RowMatrixA = Problem_->GetMatrix();
00470   if (RowMatrixA == 0)
00471     AMESOS_CHK_ERR(-1);
00472 
00473   // MS // some checks on matrix size
00474   if (RowMatrixA->NumGlobalRows() != RowMatrixA->NumGlobalCols())
00475     AMESOS_CHK_ERR(-1);
00476 
00477   //  Convert vector b to a vector in the form that DSCPACK needs it
00478   //
00479   Epetra_MultiVector* vecX = Problem_->GetLHS(); 
00480   Epetra_MultiVector* vecB = Problem_->GetRHS(); 
00481 
00482   if ((vecX == 0) || (vecB == 0))
00483     AMESOS_CHK_ERR(-1); // something wrong with input
00484 
00485   int NumVectors = vecX->NumVectors(); 
00486   if (NumVectors != vecB->NumVectors())
00487     AMESOS_CHK_ERR(-2);
00488 
00489   double *dscmapXvalues ;
00490   int dscmapXlda ;
00491   Epetra_MultiVector dscmapX(DscRowMap(),NumVectors) ; 
00492   int ierr;
00493   AMESOS_CHK_ERR(dscmapX.ExtractView(&dscmapXvalues,&dscmapXlda));
00494   assert (dscmapXlda == NumLocalCols); 
00495 
00496   double *dscmapBvalues ;
00497   int dscmapBlda ;
00498   Epetra_MultiVector dscmapB(DscRowMap(), NumVectors ) ; 
00499   ierr = dscmapB.ExtractView( &dscmapBvalues, &dscmapBlda );
00500   AMESOS_CHK_ERR(ierr);
00501   assert( dscmapBlda == NumLocalCols ) ; 
00502 
00503   AMESOS_CHK_ERR(dscmapB.Import(*vecB, Importer(), Insert));
00504 
00505   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
00506   ResetTimer(0);
00507   
00508   // MS // now solve the problem
00509   
00510   std::vector<double> ValuesInNewOrder( NumLocalCols ) ;
00511 
00512   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00513 
00514   if ( MyDscRank >= 0 ) {
00515     for ( int j =0 ; j < NumVectors; j++ ) { 
00516       for ( int i = 0; i < NumLocalCols; i++ ) { 
00517   ValuesInNewOrder[i] = dscmapBvalues[DscColMap().LID( LocalStructOldNum[i] ) +j*dscmapBlda ] ;
00518       }
00519       AMESOS_CHK_ERR( DSC_InputRhsLocalVec ( PrivateDscpackData_->MyDSCObject_, &ValuesInNewOrder[0], NumLocalCols ) ) ;
00520       AMESOS_CHK_ERR( DSC_Solve ( PrivateDscpackData_->MyDSCObject_ ) ) ; 
00521       AMESOS_CHK_ERR( DSC_GetLocalSolution ( PrivateDscpackData_->MyDSCObject_, &ValuesInNewOrder[0], NumLocalCols ) ) ; 
00522       for ( int i = 0; i < NumLocalCols; i++ ) { 
00523   dscmapXvalues[DscColMap().LID( LocalStructOldNum[i] ) +j*dscmapXlda ] = ValuesInNewOrder[i];
00524       }
00525     }
00526   }
00527 
00528   SolveTime_ = AddTime("Total solve time", SolveTime_, 0);
00529   ResetTimer(0);
00530   ResetTimer(1);
00531 
00532   vecX->Export( dscmapX, Importer(), Insert ) ;
00533 
00534   VecRedistTime_ = AddTime("Total vector redistribution time", VecRedistTime_, 0);
00535 
00536   if (ComputeTrueResidual_)
00537     ComputeTrueResidual(*(GetProblem()->GetMatrix()), *vecX, *vecB, 
00538                         false, "Amesos_Dscpack");
00539 
00540   if (ComputeVectorNorms_)
00541     ComputeVectorNorms(*vecX, *vecB, "Amesos_Dscpack");
00542   OverheadTime_ = AddTime("Total Amesos overhead time", OverheadTime_, 1);
00543   
00544   NumSolve_++;
00545 
00546   return(0) ; 
00547 }
00548 
00549 // ======================================================================
00550 void Amesos_Dscpack::PrintStatus() const
00551 {
00552   if (Problem_->GetOperator() != 0 && MyPID_ != 0)
00553   {
00554     std::string p = "Amesos_Dscpack : ";
00555     PrintLine();
00556 
00557     int n = GetProblem()->GetMatrix()->NumGlobalRows();
00558     int nnz = GetProblem()->GetMatrix()->NumGlobalNonzeros();
00559 
00560     std::cout << p << "Matrix has " << n << " rows"
00561          << " and " << nnz << " nonzeros" << std::endl;
00562     std::cout << p << "Nonzero elements per row = "
00563          << 1.0 *  nnz / n << std::endl;
00564     std::cout << p << "Percentage of nonzero elements = "
00565          << 100.0 * nnz /(pow(n,2.0)) << std::endl;
00566     std::cout << p << "Available process(es) = " << NumProcs_ << std::endl;
00567     std::cout << p << "Process(es) used = " << DscNumProcs
00568          << ", idle = " << NumProcs_ - DscNumProcs << std::endl;
00569     std::cout << p << "Estimated total memory for factorization =  " 
00570          << TotalMemory_ << " Mbytes" << std::endl; 
00571   }
00572 
00573   if ( MyDscRank >= 0 ) 
00574     DSC_DoStats( PrivateDscpackData_->MyDSCObject_ );
00575 
00576   if (!MyPID_)
00577     PrintLine();
00578 
00579   return;
00580 }
00581 
00582 // ====================================================================== 
00583 void Amesos_Dscpack::PrintTiming() const
00584 {
00585   if (Problem_->GetOperator() == 0 || MyPID_ != 0)
00586     return;
00587 
00588   double ConTime = GetTime(MtxConvTime_);
00589   double MatTime = GetTime(MtxRedistTime_);
00590   double VecTime = GetTime(VecRedistTime_);
00591   double SymTime = GetTime(SymFactTime_);
00592   double NumTime = GetTime(NumFactTime_);
00593   double SolTime = GetTime(SolveTime_);
00594   double OveTime = GetTime(OverheadTime_);
00595 
00596   if (NumSymbolicFact_)
00597     SymTime /= NumSymbolicFact_;
00598 
00599   if (NumNumericFact_)
00600     NumTime /= NumNumericFact_;
00601 
00602   if (NumSolve_)
00603     SolTime /= NumSolve_;
00604 
00605   std::string p = "Amesos_Dscpack : ";
00606   PrintLine();
00607 
00608   std::cout << p << "Time to convert matrix to DSCPACK format = "
00609        << ConTime << " (s)" << std::endl;
00610   std::cout << p << "Time to redistribute matrix = "
00611        << MatTime << " (s)" << std::endl;
00612   std::cout << p << "Time to redistribute vectors = "
00613        << VecTime << " (s)" << std::endl;
00614   std::cout << p << "Number of symbolic factorizations = "
00615        << NumSymbolicFact_ << std::endl;
00616   std::cout << p << "Time for sym fact = "
00617        << SymTime * NumSymbolicFact_ << " (s), avg = " << SymTime << " (s)" << std::endl;
00618   std::cout << p << "Number of numeric factorizations = "
00619        << NumNumericFact_ << std::endl;
00620   std::cout << p << "Time for num fact = "
00621        << NumTime * NumNumericFact_ << " (s), avg = " << NumTime << " (s)" << std::endl;
00622   std::cout << p << "Number of solve phases = "
00623        << NumSolve_ << std::endl;
00624   std::cout << p << "Time for solve = "
00625        << SolTime * NumSolve_ << " (s), avg = " << SolTime << " (s)" << std::endl;
00626 
00627   double tt = SymTime * NumSymbolicFact_ + NumTime * NumNumericFact_ + SolTime * NumSolve_;
00628   if (tt != 0)
00629   {
00630     std::cout << p << "Total time spent in Amesos = " << tt << " (s) " << std::endl;
00631     std::cout << p << "Total time spent in the Amesos interface = " << OveTime << " (s)" << std::endl;
00632     std::cout << p << "(the above time does not include DSCPACK time)" << std::endl;
00633     std::cout << p << "Amesos interface time / total time = " << OveTime / tt << std::endl;
00634   }
00635 
00636   PrintLine();
00637 
00638   return;
00639 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines