MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_DirectSparseSolverSuperLU.cpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifdef SPARSE_SOLVER_PACK_USE_SUPERLU
00043 
00044 #include <assert.h>
00045 
00046 #include <fstream>
00047 #include <algorithm>
00048 
00049 #include "AbstractLinAlgPack_DirectSparseSolverSuperLU.hpp"
00050 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00051 #include "DenseLinAlgPack_PermVecMat.hpp"
00052 #include "Teuchos_AbstractFactoryStd.hpp"
00053 #include "Teuchos_Assert.hpp"
00054 #include "Teuchos_Workspace.hpp"
00055 #include "Teuchos_dyn_cast.hpp"
00056 
00057 namespace {
00058 
00059 // Convert to compressed sparse row
00060 void convet_to_csr(
00061   int                                   n
00062   ,int                                  m
00063   ,int                                  nz
00064   ,const DenseLinAlgPack::value_type    a_val[]
00065   ,const DenseLinAlgPack::index_type    a_row_i[]
00066   ,const DenseLinAlgPack::index_type    a_col_j[]
00067   ,double                               acsr_val[]
00068   ,int                                  acsr_col_j[]
00069   ,int                                  acsr_row_ptr[]
00070   )
00071 {
00072   // Count the number of entries per row and put in acsr_row_ptr[1...m+1]
00073   std::fill_n( &acsr_row_ptr[0], m+1, 0 );
00074   {for( int k = 0; k < nz; ++k ) {
00075     ++acsr_row_ptr[a_row_i[k]]; // a_row_i[] is 1-based so this works out.
00076   }}
00077 
00078   // Transform the counts of entries per row into the start pointers for the rows.
00079   // We will make acsr_row_ptr[0] = 0 and then add form there.  We will then
00080   // shift this data so that acsr_row_ptr[1] = 0. This data
00081   // structure will be used to fill the entries per row.
00082   acsr_row_ptr[0] = 0;
00083   {for( int i = 2; i < m + 1; ++i ) {
00084     acsr_row_ptr[i] += acsr_row_ptr[i-1];
00085   }}
00086   {for( int i = m; i > 0; --i ) {
00087     acsr_row_ptr[i] = acsr_row_ptr[i-1];
00088   }}
00089 
00090   // Now copy into the compressed sparse row data structure
00091   {for( int k = 0; k < nz; ++k ) {
00092     const int row_i   = a_row_i[k];            // one-based
00093     const int row_ptr = acsr_row_ptr[row_i];  // returned value is zero-based
00094     acsr_val[row_ptr]   = a_val[k];
00095     acsr_col_j[row_ptr] = a_col_j[row_ptr] - 1; // from one-based to zero-based
00096     ++acsr_row_ptr[row_i];
00097   }}
00098   TEUCHOS_TEST_FOR_EXCEPT( !(  acsr_row_ptr[m] == nz  ) );
00099 
00100 }
00101 
00102 } // end namespace
00103 
00104 namespace AbstractLinAlgPack {
00105 
00106 //
00107 // Implementation of DirectSparseSolver(Imp) interface using SuperLU.
00108 //
00109 // Here are some little bits of knowledge about SuperLU that I need
00110 // to record after many hours of hard work.
00111 //
00112 // ToDo: Finish this!
00113 //
00114 
00115 // ToDo:
00116 // a) Add an option for printing the values of the common
00117 //    block parameters to out or to a file.  This can
00118 //    be used to get a feel for the performance of
00119 //    ma28
00120 // b) Add provisions for an external client to change
00121 //    the control options of SuperLU.  Most of these are
00122 //    stored as common block variables.
00123 
00124 // //////////////////////////////////////////////////
00125 // DirectSparseSolverSuperLU::BasisMatrixSuperLU
00126 
00127 // Overridden from BasisMatrixImp
00128 
00129 Teuchos::RCP<DirectSparseSolverImp::BasisMatrixImp>
00130 DirectSparseSolverSuperLU::BasisMatrixSuperLU::create_matrix() const
00131 {
00132   return Teuchos::rcp(new BasisMatrixSuperLU);
00133 }
00134 
00135 void DirectSparseSolverSuperLU::BasisMatrixSuperLU::V_InvMtV(
00136   VectorMutable* y, BLAS_Cpp::Transp M_trans, const Vector& x
00137   ) const 
00138 {
00139   using Teuchos::dyn_cast;
00140   using Teuchos::Workspace;
00141   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00142   size_type k;
00143 
00144   // Get concrete objects
00145   const FactorizationStructureSuperLU
00146     &fs = dyn_cast<const FactorizationStructureSuperLU>(*this->get_fact_struc());
00147   const FactorizationNonzerosSuperLU
00148     &fn = dyn_cast<const FactorizationNonzerosSuperLU>(*this->get_fact_nonzeros());
00149 
00150   VectorDenseMutableEncap  yd(*y);
00151   VectorDenseEncap         xd(x);
00152 
00153   yd() = xd(); // Copy rhs into lhs for SuperLU
00154 
00155   fs.superlu_solver_->solve(
00156     *fs.fact_struct_
00157     ,*fn.fact_nonzeros_
00158     ,M_trans == BLAS_Cpp::no_trans
00159     ,yd().dim()
00160     ,1
00161     ,&yd()[0]
00162     ,yd().dim()
00163     );
00164 
00165 }
00166 
00167 // //////////////////////////////////////////////////
00168 // DirectSparseSolverSuperLU::FactorizationStructureSuperLU
00169 
00170 DirectSparseSolverSuperLU::FactorizationStructureSuperLU::FactorizationStructureSuperLU()
00171   :superlu_solver_(SuperLUPack::SuperLUSolver::create_solver())
00172 {}
00173 
00174 // //////////////////////////////////////////////////
00175 // DirectSparseSolverSuperLU
00176 
00177 // Constructors/initializers
00178 
00179 DirectSparseSolverSuperLU::DirectSparseSolverSuperLU()
00180 {}
00181 
00182 // Overridden from DirectSparseSolver
00183 
00184 const DirectSparseSolver::basis_matrix_factory_ptr_t
00185 DirectSparseSolverSuperLU::basis_matrix_factory() const
00186 {
00187   namespace mmp = MemMngPack;
00188   return Teuchos::rcp(new Teuchos::AbstractFactoryStd<BasisMatrix,BasisMatrixSuperLU>());
00189 }
00190 
00191 void DirectSparseSolverSuperLU::estimated_fillin_ratio(
00192   value_type estimated_fillin_ratio
00193   )
00194 {
00195   TEUCHOS_TEST_FOR_EXCEPT(true);
00196 }
00197 
00198 // Overridden from DirectSparseSolverImp
00199 
00200 const Teuchos::RCP<DirectSparseSolver::FactorizationStructure>
00201 DirectSparseSolverSuperLU::create_fact_struc() const
00202 {
00203   return Teuchos::rcp(new FactorizationStructureSuperLU);
00204 }
00205 
00206 const Teuchos::RCP<DirectSparseSolverImp::FactorizationNonzeros>
00207 DirectSparseSolverSuperLU::create_fact_nonzeros() const
00208 {
00209   return Teuchos::rcp(new FactorizationNonzerosSuperLU);
00210 }
00211 
00212 void DirectSparseSolverSuperLU::imp_analyze_and_factor(
00213   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00214   ,FactorizationStructure                         *fact_struc
00215   ,FactorizationNonzeros                          *fact_nonzeros
00216   ,DenseLinAlgPack::IVector                            *row_perm
00217   ,DenseLinAlgPack::IVector                            *col_perm
00218   ,size_type                                      *rank
00219   ,std::ostream                                   *out
00220   )
00221 {
00222   namespace mmp = MemMngPack;
00223   using Teuchos::dyn_cast;
00224   typedef MatrixConvertToSparse MCTS;
00225   using Teuchos::Workspace;
00226   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00227 
00228   if(out)
00229     *out << "\nUsing SuperLU to analyze and factor a new matrix ...\n";
00230 
00231   // Get the concrete factorization and nonzeros objects
00232   FactorizationStructureSuperLU
00233     &fs = dyn_cast<FactorizationStructureSuperLU>(*fact_struc);
00234   FactorizationNonzerosSuperLU
00235     &fn = dyn_cast<FactorizationNonzerosSuperLU>(*fact_nonzeros);
00236 
00237   // Allocate new storage if not done so already
00238   if(!fs.fact_struct_.get())
00239     fs.fact_struct_ = SuperLUPack::SuperLUSolver::create_fact_struct();
00240   if(!fn.fact_nonzeros_.get())
00241     fn.fact_nonzeros_ = SuperLUPack::SuperLUSolver::create_fact_nonzeros();
00242 
00243   // Get the dimensions of things.
00244   const index_type
00245     m = A.rows(),
00246     n = A.cols(),
00247     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00248 
00249   // Validate input
00250   TEUCHOS_TEST_FOR_EXCEPTION(
00251     n <= 0 || m <= 0 || m > n, std::invalid_argument
00252     ,"DirectSparseSolverSuperLU::imp_analyze_and_factor(...) : Error!" );
00253 
00254   // Extract the matrix in coordinate format
00255   Workspace<value_type>   a_val(wss,nz);
00256   Workspace<index_type>   a_row_i(wss,nz);
00257   Workspace<index_type>   a_col_j(wss,nz);
00258   A.coor_extract_nonzeros(
00259     MCTS::EXTRACT_FULL_MATRIX
00260     ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00261     ,nz
00262     ,&a_val[0]
00263     ,nz
00264     ,&a_row_i[0]
00265     ,&a_col_j[0]
00266     );
00267 
00268   //
00269   // Convert to compressed sparse row format (which is compressed sparse
00270   // column of the transposed matrix which will be passed to
00271   // SuperLU for factorization).
00272   //
00273 
00274   Workspace<double>   acsr_val(wss,nz);
00275   Workspace<int>      acsr_col_j(wss,nz);    // Zero-based for SuperLU
00276   Workspace<int>      acsr_row_ptr(wss,m+1);
00277   
00278   convet_to_csr(
00279     n,m,nz,&a_val[0],&a_row_i[0],&a_col_j[0]
00280     ,&acsr_val[0],&acsr_col_j[0],&acsr_row_ptr[0]
00281     );
00282 
00283   //
00284   // Have SuperLU factor this matrix.
00285   //
00286   // SuperLU works with the transpose of the matrix
00287   // That DirectSparseSolver works with.
00288   //
00289   
00290   Workspace<int>    perm_r(wss,m); // Zero-based for SuperLU
00291   Workspace<int>    perm_c(wss,n); // Zero-based for SuperLU
00292   int                    slu_rank = 0;
00293 
00294   fs.superlu_solver_->analyze_and_factor(
00295     n                         // m
00296     ,m                        // n
00297     ,nz                       // nz
00298     ,&acsr_val[0]             // a_val[]
00299     ,&acsr_col_j[0]           // a_row_i[]
00300     ,&acsr_row_ptr[0]         // a_col_ptr[]
00301     ,fs.fact_struct_.get()    // fact_struct
00302     ,fn.fact_nonzeros_.get()  // fact_nonzeros
00303     ,&perm_c[0]               // perm_r[]
00304     ,&perm_r[0]               // perm_c[]
00305     ,&slu_rank                // rank
00306     );
00307 
00308   // Copy the data to the output
00309   row_perm->resize(m);
00310   for( int i = 0; i < m; ++i )
00311     (*row_perm)[i] = perm_r[i] + 1; // Convert from zero based to one based
00312   col_perm->resize(n);
00313   for( int j = 0; j < n; ++j )
00314     (*col_perm)[j] = perm_c[j] + 1; // Convert from zero based to one based
00315   *rank = slu_rank;
00316 
00317   // Sort partitions into assending order (required!)
00318   std::sort(&(*row_perm)[0] , &(*row_perm)[0] + (*rank) );
00319   std::sort(&(*col_perm)[0] , &(*col_perm)[0] + (*rank) );
00320   if( *rank < m )
00321     std::sort(&(*row_perm)[0] + (*rank) , &(*row_perm)[0] + m );
00322   if( *rank < n )
00323     std::sort(&(*col_perm)[0] + (*rank) , &(*col_perm)[0] + n );
00324 
00325 }
00326 
00327 void DirectSparseSolverSuperLU::imp_factor(
00328   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00329   ,const FactorizationStructure                   &fact_struc
00330   ,FactorizationNonzeros                          *fact_nonzeros
00331   ,std::ostream                                   *out
00332   )
00333 {
00334   namespace mmp = MemMngPack;
00335   using Teuchos::dyn_cast;
00336   typedef MatrixConvertToSparse MCTS;
00337   using Teuchos::Workspace;
00338   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00339 
00340   if(out)
00341     *out << "\nUsing SuperLU to refactor the basis matrix ...\n";
00342 
00343   // Get the concrete factorization and nonzeros objects
00344   const FactorizationStructureSuperLU
00345     &fs = dyn_cast<const FactorizationStructureSuperLU>(fact_struc);
00346   FactorizationNonzerosSuperLU
00347     &fn = dyn_cast<FactorizationNonzerosSuperLU>(*fact_nonzeros);
00348 
00349   // Allocate new storage if not done so already
00350   TEUCHOS_TEST_FOR_EXCEPTION(
00351     !fs.fact_struct_.get(), std::logic_error
00352     ,"DirectSparseSolverSuperLU::imp_factor(...): Error, the factorization sturcture must "
00353     "have already been computed!"
00354     );
00355   if(!fn.fact_nonzeros_.get())
00356     fn.fact_nonzeros_ = SuperLUPack::SuperLUSolver::create_fact_nonzeros();
00357 
00358   // Get the dimensions of things.
00359   const index_type
00360     m = A.rows(),
00361     n = A.cols(),
00362     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00363 
00364   // Validate input
00365   TEUCHOS_TEST_FOR_EXCEPTION(
00366     n <= 0 || m <= 0 || m > n, std::invalid_argument
00367     ,"DirectSparseSolverSuperLU::imp_factor(...) : Error!" );
00368 
00369   // Extract the matrix in coordinate format
00370   Workspace<value_type>   a_val(wss,nz);
00371   Workspace<index_type>   a_row_i(wss,nz);
00372   Workspace<index_type>   a_col_j(wss,nz);
00373   A.coor_extract_nonzeros(
00374     MCTS::EXTRACT_FULL_MATRIX
00375     ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00376     ,nz
00377     ,&a_val[0]
00378     ,nz
00379     ,&a_row_i[0]
00380     ,&a_col_j[0]
00381     );
00382 
00383   //
00384   // Convert to compressed sparse row format (which is compressed sparse
00385   // column of the transposed matrix which will be passed to
00386   // SuperLU for factorization).
00387   //
00388 
00389   Workspace<double>   acsr_val(wss,nz);
00390   Workspace<int>      acsr_col_j(wss,nz);    // Zero-based for SuperLU
00391   Workspace<int>      acsr_row_ptr(wss,m+1);
00392   
00393   convet_to_csr(
00394     n,m,nz,&a_val[0],&a_row_i[0],&a_col_j[0]
00395     ,&acsr_val[0],&acsr_col_j[0],&acsr_row_ptr[0]
00396     );
00397 
00398   //
00399   // Have SuperLU factor this matrix.
00400   //
00401   // SuperLU works with the transpose of the matrix
00402   // That DirectSparseSolver works with.
00403   //
00404 
00405   fs.superlu_solver_->factor(
00406     n                         // m
00407     ,m                        // n
00408     ,nz                       // nz
00409     ,&acsr_val[0]             // a_val[]
00410     ,&acsr_col_j[0]           // a_row_i[]
00411     ,&acsr_row_ptr[0]         // a_col_ptr[]
00412     ,*fs.fact_struct_         // fact_struct
00413     ,fn.fact_nonzeros_.get()  // fact_nonzeros
00414     );
00415 
00416 }
00417 
00418 // private
00419 
00420 } // end namespace AbstractLinAlgPack 
00421 
00422 #endif // SPARSE_SOLVER_PACK_USE_SUPERLU
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines