Amesos Package Browser (Single Doxygen Collection) Development
Epetra_SLU.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 "Epetra_SLU.h"
00030 
00031 namespace SLU
00032 {
00033   //FIXME
00034   typedef int int_t;
00035   
00036 extern "C" {
00037 #include "dsp_defs.h"
00038 }
00039 }
00040 
00041 #include "Epetra_CrsGraph.h"
00042 #include "Epetra_LinearProblem.h"
00043 #include "Epetra_MultiVector.h"
00044 #include "Epetra_CrsMatrix.h"
00045 
00046 struct SLUData
00047 {
00048   SLU::SuperMatrix A, B, X, L, U;
00049   SLU::factor_param_t iparam;
00050   SLU::mem_usage_t mem_usage;
00051 };
00052 
00053 Epetra_SLU::Epetra_SLU( Epetra_LinearProblem * Problem,
00054                         int fill_fac,
00055                         int panel_size,
00056                         int relax )
00057 : count_(0)
00058 {
00059   using namespace SLU;
00060 
00061   data_ = new SLUData();
00062 
00063   data_->iparam.panel_size = panel_size;
00064   data_->iparam.relax      = relax;
00065   data_->iparam.diag_pivot_thresh = -1;
00066   data_->iparam.drop_tol = -1;
00067 
00068   A_ = dynamic_cast<Epetra_CrsMatrix *> (Problem->GetOperator());
00069   X_ = Problem->GetLHS();
00070   B_ = Problem->GetRHS();
00071 
00072   AG_ = new Epetra_CrsGraph( A_->Graph() );
00073 
00074   int NumIndices;
00075   int NumMyCols = AG_->NumMyCols();
00076   int NumMyEqs = A_->NumMyRows();
00077   int GlobalMaxNumIndices = AG_->GlobalMaxNumIndices();
00078   int * xIndices;
00079 
00080   TransNumNZ_ = new int[NumMyCols];
00081   for( int i = 0; i < NumMyCols; i++ )
00082     TransNumNZ_[i] = 0;
00083   for( int i = 0; i < NumMyEqs; i++ )
00084   {
00085     assert( AG_->ExtractMyRowView( i, NumIndices, xIndices ) == 0 );
00086     for( int j = 0; j < NumIndices; j++ )
00087       ++TransNumNZ_[ xIndices[j] ];
00088   }
00089   TotNumNZ_ = 0;
00090   for( int i = 0; i < NumMyCols; i++ )
00091    TotNumNZ_ += TransNumNZ_[i];
00092 
00093   TransIndices_ = new int*[ NumMyCols ];
00094   TransValues_ = new double*[ NumMyCols ];
00095 
00096   RowIndices_ = new int[TotNumNZ_];
00097   RowValues_ = new double[TotNumNZ_];
00098   ColPointers_ = new int[NumMyCols+1];
00099 
00100   // Set pointers into the RowIndices and Values arrays, define ColPointers
00101   NumIndices = 0;
00102   for(int i=0; i<NumMyCols; i++) {
00103     ColPointers_[i] = NumIndices;
00104     TransIndices_[i] = RowIndices_ + NumIndices;
00105     TransValues_[i] = RowValues_ + NumIndices;
00106     NumIndices += TransNumNZ_[i];
00107   }
00108   ColPointers_[NumMyCols] = NumIndices;
00109 
00110   Copy();
00111 
00112   int m = A_->NumGlobalRows();
00113   int n = A_->NumGlobalCols();
00114 
00115   /* Create matrix A in the format expected by SuperLU. */
00116   dCreate_CompCol_Matrix( &(data_->A), m, n, TotNumNZ_, RowValues_,
00117   RowIndices_, ColPointers_, SLU_NC, SLU_D, SLU_GE );
00118   
00119   /* Create right-hand side matrix B. */
00120   double * rhs_x;
00121   double * rhs_b;
00122   int LDA_x, LDA_b;
00123   X_->ExtractView( &rhs_x, &LDA_x );
00124   dCreate_Dense_Matrix( &(data_->X), m, 1, rhs_x, m, SLU_DN, SLU_D, SLU_GE);
00125   B_->ExtractView( &rhs_b, &LDA_b );
00126   dCreate_Dense_Matrix( &(data_->B), m, 1, rhs_b, m, SLU_DN, SLU_D, SLU_GE);
00127   
00128   perm_r_ = new int[m];
00129   perm_c_ = new int[n];
00130   
00131   etree_ = new int[n];
00132 
00133   ferr_ = new double[1];
00134   berr_ = new double[1];
00135 
00136   R_ = new double[m];
00137   C_ = new double[n];
00138 }
00139 
00140 Epetra_SLU::~Epetra_SLU()
00141 {
00142   SLU::Destroy_SuperMatrix_Store( &(data_->A) );
00143   SLU::Destroy_SuperMatrix_Store( &(data_->B) );
00144   SLU::Destroy_SuperMatrix_Store( &(data_->X) );
00145 
00146   delete [] TransIndices_;
00147   delete [] TransValues_;
00148 
00149   delete [] TransNumNZ_;
00150 
00151   delete [] RowIndices_;
00152   delete [] RowValues_;
00153   delete [] ColPointers_;
00154 
00155   delete [] perm_r_;
00156   delete [] perm_c_;
00157 
00158   delete [] etree_;
00159 
00160   delete [] ferr_;
00161   delete [] berr_;
00162 
00163   delete [] R_;
00164   delete [] C_;
00165 
00166 //  Destroy_CompCol_Matrix(data_->A);
00167 //  SLU::Destroy_SuperNode_Matrix( &(data_->L) );
00168 //  SLU::Destroy_CompCol_Matrix( &(data_->U) );
00169 
00170   delete data_;
00171 
00172   delete AG_;
00173 }
00174 
00175 void Epetra_SLU::Copy()
00176 {
00177   int NumMyCols = AG_->NumMyCols();
00178   int NumMyEqs = A_->NumMyRows();
00179   int GlobalMaxNumIndices = AG_->GlobalMaxNumIndices();
00180   int NumIndices;
00181   int * xIndices;
00182   double * xValues;
00183 
00184   for ( int i = 0; i < NumMyCols; i++ )
00185     TransNumNZ_[i] = 0;
00186 
00187   for ( int i = 0; i < NumMyEqs; i++ )
00188   {
00189     assert(A_->ExtractMyRowView( i, NumIndices, xValues, xIndices) == 0 );
00190     int ii = A_->GRID(i);
00191     for ( int j = 0; j < NumIndices; j++ )
00192     {
00193       int TransRow = xIndices[j];
00194       int loc = TransNumNZ_[TransRow];
00195       TransIndices_[TransRow][loc] = ii;
00196       TransValues_[TransRow][loc] = xValues[j];
00197       ++TransNumNZ_[TransRow]; // increment counter into current transpose row
00198     }
00199   }
00200 }
00201 
00202 int Epetra_SLU::Solve( bool Verbose,
00203                        bool Equil,
00204                        bool Factor,
00205                        int perm_type,
00206                        double pivot_thresh,
00207                        bool Refact,
00208                        bool Trans )
00209 {
00210   Copy();
00211 
00212   using namespace SLU;
00213 
00214   int m = A_->NumGlobalRows();
00215 
00216   int permt = perm_type;
00217   if( m < 3 ) permt = 0;
00218 
00219   /*
00220    * Get column permutation vector perm_c[], according to permc_spec:
00221    *   permc_spec = 0: use the natural ordering
00222    *   permc_spec = 1: use minimum degree ordering on structure of A'*A
00223    *   permc_spec = 2: use minimum degree ordering on structure of A'+A
00224    */
00225   if( !count_ ) get_perm_c( permt, &(data_->A), perm_c_ );
00226 
00227   if( Verbose ) cout << "MATRIX COPIED!" << endl;
00228 
00229   int info = 0;
00230 
00231   char fact, trans, refact, equed;
00232 
00233   if( Trans ) trans = 'T';
00234   else        trans = 'N';
00235 
00236   if( Equil ) fact = 'E';
00237   else        fact = 'N';
00238 
00239   if( !count_ || !Refact )
00240   {
00241     refact = 'N';
00242     count_ = 0;
00243   }
00244   else
00245   {
00246     refact = 'Y';
00247     if( !Factor ) fact = 'F';
00248   }
00249 
00250   if( Equil ) equed = 'B';
00251   else        equed = 'N';
00252 
00253   data_->iparam.diag_pivot_thresh = pivot_thresh;
00254 
00255   double rpg, rcond;
00256  
00257   if( Verbose ) cout << "TRANS:  " << trans << endl;
00258   if( Verbose ) cout << "REFACT: " << refact << endl;
00259 
00260   dgssvx( &fact, &trans, &refact, &(data_->A), &(data_->iparam), perm_c_,
00261   perm_r_, etree_, &equed, R_, C_, &(data_->L), &(data_->U),
00262   NULL, 0, &(data_->B), &(data_->X), &rpg, &rcond,
00263   ferr_, berr_, &(data_->mem_usage), &info );
00264 
00265   if( info ) cout << "WARNING: SuperLU returned with error code = " << info << endl;
00266 
00267   if( Verbose )
00268   {
00269     cout << "SYSTEM DIRECT SOLVED!" << endl;
00270 
00271     cout << "SuperLU INFO: " << info << "\n\n";
00272     cout << "  ncols = " << A_->NumGlobalCols() << endl;
00273 
00274     cout << "SuperLU Memory Usage\n";
00275     cout << "--------------------\n";
00276     cout << "LU: " << data_->mem_usage.for_lu << endl;
00277     cout << "Total: " << data_->mem_usage.total_needed << endl;
00278     cout << "Expands: " << data_->mem_usage.expansions << endl;
00279     cout << "--------------------\n\n";
00280   
00281     if (m<200) dPrint_CompCol_Matrix("A", &(data_->A));
00282 //    if (m<200) dPrint_CompCol_Matrix("U", &(data_->U));
00283 //    if (m<200) dPrint_SuperNode_Matrix("L", &(data_->L));
00284 //  if (m<200) PrintInt10("\nperm_r", m, perm_r);
00285     if (m<200) dPrint_Dense_Matrix("B", &(data_->B));
00286     if (m<200) dPrint_Dense_Matrix("X", &(data_->X));
00287   }
00288 
00289   count_++;
00290 
00291   return 0;
00292 }
00293 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines