AbstractLinAlgPack_DirectSparseSolverSuperLU.cpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) 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 Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifdef SPARSE_SOLVER_PACK_USE_SUPERLU
00030 
00031 #include <assert.h>
00032 
00033 #include <fstream>
00034 #include <algorithm>
00035 
00036 #include "AbstractLinAlgPack_DirectSparseSolverSuperLU.hpp"
00037 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00038 #include "DenseLinAlgPack_PermVecMat.hpp"
00039 #include "Teuchos_AbstractFactoryStd.hpp"
00040 #include "Teuchos_TestForException.hpp"
00041 #include "Teuchos_Workspace.hpp"
00042 #include "Teuchos_dyn_cast.hpp"
00043 
00044 namespace {
00045 
00046 // Convert to compressed sparse row
00047 void convet_to_csr(
00048   int                                   n
00049   ,int                                  m
00050   ,int                                  nz
00051   ,const DenseLinAlgPack::value_type    a_val[]
00052   ,const DenseLinAlgPack::index_type    a_row_i[]
00053   ,const DenseLinAlgPack::index_type    a_col_j[]
00054   ,double                               acsr_val[]
00055   ,int                                  acsr_col_j[]
00056   ,int                                  acsr_row_ptr[]
00057   )
00058 {
00059   // Count the number of entries per row and put in acsr_row_ptr[1...m+1]
00060   std::fill_n( &acsr_row_ptr[0], m+1, 0 );
00061   {for( int k = 0; k < nz; ++k ) {
00062     ++acsr_row_ptr[a_row_i[k]]; // a_row_i[] is 1-based so this works out.
00063   }}
00064 
00065   // Transform the counts of entries per row into the start pointers for the rows.
00066   // We will make acsr_row_ptr[0] = 0 and then add form there.  We will then
00067   // shift this data so that acsr_row_ptr[1] = 0. This data
00068   // structure will be used to fill the entries per row.
00069   acsr_row_ptr[0] = 0;
00070   {for( int i = 2; i < m + 1; ++i ) {
00071     acsr_row_ptr[i] += acsr_row_ptr[i-1];
00072   }}
00073   {for( int i = m; i > 0; --i ) {
00074     acsr_row_ptr[i] = acsr_row_ptr[i-1];
00075   }}
00076 
00077   // Now copy into the compressed sparse row data structure
00078   {for( int k = 0; k < nz; ++k ) {
00079     const int row_i   = a_row_i[k];            // one-based
00080     const int row_ptr = acsr_row_ptr[row_i];  // returned value is zero-based
00081     acsr_val[row_ptr]   = a_val[k];
00082     acsr_col_j[row_ptr] = a_col_j[row_ptr] - 1; // from one-based to zero-based
00083     ++acsr_row_ptr[row_i];
00084   }}
00085   assert( acsr_row_ptr[m] == nz );
00086 
00087 }
00088 
00089 } // end namespace
00090 
00091 namespace AbstractLinAlgPack {
00092 
00093 //
00094 // Implementation of DirectSparseSolver(Imp) interface using SuperLU.
00095 //
00096 // Here are some little bits of knowledge about SuperLU that I need
00097 // to record after many hours of hard work.
00098 //
00099 // ToDo: Finish this!
00100 //
00101 
00102 // ToDo:
00103 // a) Add an option for printing the values of the common
00104 //    block parameters to out or to a file.  This can
00105 //    be used to get a feel for the performance of
00106 //    ma28
00107 // b) Add provisions for an external client to change
00108 //    the control options of SuperLU.  Most of these are
00109 //    stored as common block variables.
00110 
00111 // //////////////////////////////////////////////////
00112 // DirectSparseSolverSuperLU::BasisMatrixSuperLU
00113 
00114 // Overridden from BasisMatrixImp
00115 
00116 Teuchos::RefCountPtr<DirectSparseSolverImp::BasisMatrixImp>
00117 DirectSparseSolverSuperLU::BasisMatrixSuperLU::create_matrix() const
00118 {
00119   return Teuchos::rcp(new BasisMatrixSuperLU);
00120 }
00121 
00122 void DirectSparseSolverSuperLU::BasisMatrixSuperLU::V_InvMtV(
00123   VectorMutable* y, BLAS_Cpp::Transp M_trans, const Vector& x
00124   ) const 
00125 {
00126   using Teuchos::dyn_cast;
00127   using Teuchos::Workspace;
00128   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00129   size_type k;
00130 
00131   // Get concrete objects
00132   const FactorizationStructureSuperLU
00133     &fs = dyn_cast<const FactorizationStructureSuperLU>(*this->get_fact_struc());
00134   const FactorizationNonzerosSuperLU
00135     &fn = dyn_cast<const FactorizationNonzerosSuperLU>(*this->get_fact_nonzeros());
00136 
00137   VectorDenseMutableEncap  yd(*y);
00138   VectorDenseEncap         xd(x);
00139 
00140   yd() = xd(); // Copy rhs into lhs for SuperLU
00141 
00142   fs.superlu_solver_->solve(
00143     *fs.fact_struct_
00144     ,*fn.fact_nonzeros_
00145     ,M_trans == BLAS_Cpp::no_trans
00146     ,yd().dim()
00147     ,1
00148     ,&yd()[0]
00149     ,yd().dim()
00150     );
00151 
00152 }
00153 
00154 // //////////////////////////////////////////////////
00155 // DirectSparseSolverSuperLU::FactorizationStructureSuperLU
00156 
00157 DirectSparseSolverSuperLU::FactorizationStructureSuperLU::FactorizationStructureSuperLU()
00158   :superlu_solver_(SuperLUPack::SuperLUSolver::create_solver())
00159 {}
00160 
00161 // //////////////////////////////////////////////////
00162 // DirectSparseSolverSuperLU
00163 
00164 // Constructors/initializers
00165 
00166 DirectSparseSolverSuperLU::DirectSparseSolverSuperLU()
00167 {}
00168 
00169 // Overridden from DirectSparseSolver
00170 
00171 const DirectSparseSolver::basis_matrix_factory_ptr_t
00172 DirectSparseSolverSuperLU::basis_matrix_factory() const
00173 {
00174   namespace mmp = MemMngPack;
00175   return Teuchos::rcp(new Teuchos::AbstractFactoryStd<BasisMatrix,BasisMatrixSuperLU>());
00176 }
00177 
00178 void DirectSparseSolverSuperLU::estimated_fillin_ratio(
00179   value_type estimated_fillin_ratio
00180   )
00181 {
00182   TEST_FOR_EXCEPT(true);
00183 }
00184 
00185 // Overridden from DirectSparseSolverImp
00186 
00187 const Teuchos::RefCountPtr<DirectSparseSolver::FactorizationStructure>
00188 DirectSparseSolverSuperLU::create_fact_struc() const
00189 {
00190   return Teuchos::rcp(new FactorizationStructureSuperLU);
00191 }
00192 
00193 const Teuchos::RefCountPtr<DirectSparseSolverImp::FactorizationNonzeros>
00194 DirectSparseSolverSuperLU::create_fact_nonzeros() const
00195 {
00196   return Teuchos::rcp(new FactorizationNonzerosSuperLU);
00197 }
00198 
00199 void DirectSparseSolverSuperLU::imp_analyze_and_factor(
00200   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00201   ,FactorizationStructure                         *fact_struc
00202   ,FactorizationNonzeros                          *fact_nonzeros
00203   ,DenseLinAlgPack::IVector                            *row_perm
00204   ,DenseLinAlgPack::IVector                            *col_perm
00205   ,size_type                                      *rank
00206   ,std::ostream                                   *out
00207   )
00208 {
00209   namespace mmp = MemMngPack;
00210   using Teuchos::dyn_cast;
00211   typedef MatrixConvertToSparse MCTS;
00212   using Teuchos::Workspace;
00213   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00214 
00215   if(out)
00216     *out << "\nUsing SuperLU to analyze and factor a new matrix ...\n";
00217 
00218   // Get the concrete factorization and nonzeros objects
00219   FactorizationStructureSuperLU
00220     &fs = dyn_cast<FactorizationStructureSuperLU>(*fact_struc);
00221   FactorizationNonzerosSuperLU
00222     &fn = dyn_cast<FactorizationNonzerosSuperLU>(*fact_nonzeros);
00223 
00224   // Allocate new storage if not done so already
00225   if(!fs.fact_struct_.get())
00226     fs.fact_struct_ = SuperLUPack::SuperLUSolver::create_fact_struct();
00227   if(!fn.fact_nonzeros_.get())
00228     fn.fact_nonzeros_ = SuperLUPack::SuperLUSolver::create_fact_nonzeros();
00229 
00230   // Get the dimensions of things.
00231   const index_type
00232     m = A.rows(),
00233     n = A.cols(),
00234     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00235 
00236   // Validate input
00237   TEST_FOR_EXCEPTION(
00238     n <= 0 || m <= 0 || m > n, std::invalid_argument
00239     ,"DirectSparseSolverSuperLU::imp_analyze_and_factor(...) : Error!" );
00240 
00241   // Extract the matrix in coordinate format
00242   Workspace<value_type>   a_val(wss,nz);
00243   Workspace<index_type>   a_row_i(wss,nz);
00244   Workspace<index_type>   a_col_j(wss,nz);
00245   A.coor_extract_nonzeros(
00246     MCTS::EXTRACT_FULL_MATRIX
00247     ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00248     ,nz
00249     ,&a_val[0]
00250     ,nz
00251     ,&a_row_i[0]
00252     ,&a_col_j[0]
00253     );
00254 
00255   //
00256   // Convert to compressed sparse row format (which is compressed sparse
00257   // column of the transposed matrix which will be passed to
00258   // SuperLU for factorization).
00259   //
00260 
00261   Workspace<double>   acsr_val(wss,nz);
00262   Workspace<int>      acsr_col_j(wss,nz);    // Zero-based for SuperLU
00263   Workspace<int>      acsr_row_ptr(wss,m+1);
00264   
00265   convet_to_csr(
00266     n,m,nz,&a_val[0],&a_row_i[0],&a_col_j[0]
00267     ,&acsr_val[0],&acsr_col_j[0],&acsr_row_ptr[0]
00268     );
00269 
00270   //
00271   // Have SuperLU factor this matrix.
00272   //
00273   // SuperLU works with the transpose of the matrix
00274   // That DirectSparseSolver works with.
00275   //
00276   
00277   Workspace<int>    perm_r(wss,m); // Zero-based for SuperLU
00278   Workspace<int>    perm_c(wss,n); // Zero-based for SuperLU
00279   int                    slu_rank = 0;
00280 
00281   fs.superlu_solver_->analyze_and_factor(
00282     n                         // m
00283     ,m                        // n
00284     ,nz                       // nz
00285     ,&acsr_val[0]             // a_val[]
00286     ,&acsr_col_j[0]           // a_row_i[]
00287     ,&acsr_row_ptr[0]         // a_col_ptr[]
00288     ,fs.fact_struct_.get()    // fact_struct
00289     ,fn.fact_nonzeros_.get()  // fact_nonzeros
00290     ,&perm_c[0]               // perm_r[]
00291     ,&perm_r[0]               // perm_c[]
00292     ,&slu_rank                // rank
00293     );
00294 
00295   // Copy the data to the output
00296   row_perm->resize(m);
00297   for( int i = 0; i < m; ++i )
00298     (*row_perm)[i] = perm_r[i] + 1; // Convert from zero based to one based
00299   col_perm->resize(n);
00300   for( int j = 0; j < n; ++j )
00301     (*col_perm)[j] = perm_c[j] + 1; // Convert from zero based to one based
00302   *rank = slu_rank;
00303 
00304   // Sort partitions into assending order (required!)
00305   std::sort(&(*row_perm)[0] , &(*row_perm)[0] + (*rank) );
00306   std::sort(&(*col_perm)[0] , &(*col_perm)[0] + (*rank) );
00307   if( *rank < m )
00308     std::sort(&(*row_perm)[0] + (*rank) , &(*row_perm)[0] + m );
00309   if( *rank < n )
00310     std::sort(&(*col_perm)[0] + (*rank) , &(*col_perm)[0] + n );
00311 
00312 }
00313 
00314 void DirectSparseSolverSuperLU::imp_factor(
00315   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00316   ,const FactorizationStructure                   &fact_struc
00317   ,FactorizationNonzeros                          *fact_nonzeros
00318   ,std::ostream                                   *out
00319   )
00320 {
00321   namespace mmp = MemMngPack;
00322   using Teuchos::dyn_cast;
00323   typedef MatrixConvertToSparse MCTS;
00324   using Teuchos::Workspace;
00325   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00326 
00327   if(out)
00328     *out << "\nUsing SuperLU to refactor the basis matrix ...\n";
00329 
00330   // Get the concrete factorization and nonzeros objects
00331   const FactorizationStructureSuperLU
00332     &fs = dyn_cast<const FactorizationStructureSuperLU>(fact_struc);
00333   FactorizationNonzerosSuperLU
00334     &fn = dyn_cast<FactorizationNonzerosSuperLU>(*fact_nonzeros);
00335 
00336   // Allocate new storage if not done so already
00337   TEST_FOR_EXCEPTION(
00338     !fs.fact_struct_.get(), std::logic_error
00339     ,"DirectSparseSolverSuperLU::imp_factor(...): Error, the factorization sturcture must "
00340     "have already been computed!"
00341     );
00342   if(!fn.fact_nonzeros_.get())
00343     fn.fact_nonzeros_ = SuperLUPack::SuperLUSolver::create_fact_nonzeros();
00344 
00345   // Get the dimensions of things.
00346   const index_type
00347     m = A.rows(),
00348     n = A.cols(),
00349     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00350 
00351   // Validate input
00352   TEST_FOR_EXCEPTION(
00353     n <= 0 || m <= 0 || m > n, std::invalid_argument
00354     ,"DirectSparseSolverSuperLU::imp_factor(...) : Error!" );
00355 
00356   // Extract the matrix in coordinate format
00357   Workspace<value_type>   a_val(wss,nz);
00358   Workspace<index_type>   a_row_i(wss,nz);
00359   Workspace<index_type>   a_col_j(wss,nz);
00360   A.coor_extract_nonzeros(
00361     MCTS::EXTRACT_FULL_MATRIX
00362     ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00363     ,nz
00364     ,&a_val[0]
00365     ,nz
00366     ,&a_row_i[0]
00367     ,&a_col_j[0]
00368     );
00369 
00370   //
00371   // Convert to compressed sparse row format (which is compressed sparse
00372   // column of the transposed matrix which will be passed to
00373   // SuperLU for factorization).
00374   //
00375 
00376   Workspace<double>   acsr_val(wss,nz);
00377   Workspace<int>      acsr_col_j(wss,nz);    // Zero-based for SuperLU
00378   Workspace<int>      acsr_row_ptr(wss,m+1);
00379   
00380   convet_to_csr(
00381     n,m,nz,&a_val[0],&a_row_i[0],&a_col_j[0]
00382     ,&acsr_val[0],&acsr_col_j[0],&acsr_row_ptr[0]
00383     );
00384 
00385   //
00386   // Have SuperLU factor this matrix.
00387   //
00388   // SuperLU works with the transpose of the matrix
00389   // That DirectSparseSolver works with.
00390   //
00391 
00392   fs.superlu_solver_->factor(
00393     n                         // m
00394     ,m                        // n
00395     ,nz                       // nz
00396     ,&acsr_val[0]             // a_val[]
00397     ,&acsr_col_j[0]           // a_row_i[]
00398     ,&acsr_row_ptr[0]         // a_col_ptr[]
00399     ,*fs.fact_struct_         // fact_struct
00400     ,fn.fact_nonzeros_.get()  // fact_nonzeros
00401     );
00402 
00403 }
00404 
00405 // private
00406 
00407 } // end namespace AbstractLinAlgPack 
00408 
00409 #endif // SPARSE_SOLVER_PACK_USE_SUPERLU

Generated on Thu Sep 18 12:33:51 2008 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by doxygen 1.3.9.1