MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_SuperLUSolver.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 #include <valarray>
00046 #include <stdexcept>
00047 
00048 #include "AbstractLinAlgPack_SuperLUSolver.hpp"
00049 #include "Teuchos_dyn_cast.hpp"
00050 #include "Teuchos_Workspace.hpp"
00051 #include "Teuchos_Assert.hpp"
00052 
00053 // SuperLU
00054 #include "dsp_defs.h"
00055 #include "util.h"
00056 
00057 namespace {
00058 
00059 // Static SuperLU stuff
00060 
00061 int local_panel_size  = 0;
00062 int local_relax       = 0;
00063 
00064 class StaticSuperLUInit {
00065 public:
00066   StaticSuperLUInit()
00067     {
00068       local_panel_size = sp_ienv(1);
00069       local_relax      = sp_ienv(2);
00070       StatInit(local_panel_size,local_relax);
00071     }
00072   ~StaticSuperLUInit()
00073     {
00074       StatFree();
00075     }
00076 };
00077 
00078 StaticSuperLUInit static_super_lu_init; // Will be created early and destroyed late!
00079 
00080 // ToDo: RAB: 2002/10/14: We must find a better way to work with
00081 // SuperLU than this.  It should not be too hard
00082 // to do better in the future.
00083 
00084 // A cast to const is needed because the standard does not return a reference from
00085 // valarray<>::operator[]() const.
00086 template <class T>
00087 std::valarray<T>& cva(const std::valarray<T>& va )
00088 {
00089   return const_cast<std::valarray<T>&>(va);
00090 }
00091 
00092 } // end namespace
00093 
00094 namespace SuperLUPack {
00095 
00096 class SuperLUSolverImpl;
00097 
00102 class SuperLUSolverImpl : public SuperLUSolver {
00103 public:
00104 
00107 
00109   class FactorizationStructureImpl : public FactorizationStructure {
00110   public:
00111     friend class SuperLUSolverImpl;
00112   private:
00113     int                   rank_;        // For square basis
00114     int                   nz_;          // ...
00115     std::valarray<int>    perm_r_;      // ...
00116     std::valarray<int>    perm_c_;      // ...
00117     std::valarray<int>    etree_;       // ...
00118     int                   m_orig_;      // For original rectangular matrix
00119     int                   n_orig_;      // ...
00120     int                   nz_orig_;     // ...
00121     std::valarray<int>    perm_r_orig_; // ...
00122     std::valarray<int>    perm_c_orig_; // ...
00123   };
00124 
00126   class FactorizationNonzerosImpl : public FactorizationNonzeros {
00127   public:
00128     friend class SuperLUSolverImpl;
00129   private:
00130     SuperMatrix   L_;
00131     SuperMatrix   U_;
00132   };
00133 
00135 
00138 
00140   void analyze_and_factor(
00141     int                         m
00142     ,int                        n
00143     ,int                        nz
00144     ,const double               a_val[]
00145     ,const int                  a_row_i[]
00146     ,const int                  a_col_ptr[]
00147     ,FactorizationStructure     *fact_struct
00148     ,FactorizationNonzeros      *fact_nonzeros
00149     ,int                        row_perm[]
00150     ,int                        col_perm[]
00151     ,int                        *rank
00152     );
00154   void factor(
00155     int                             m
00156     ,int                            n
00157     ,int                            nz
00158     ,const double                   a_val[]
00159     ,const int                      a_row_i[]
00160     ,const int                      a_col_ptr[]
00161     ,const FactorizationStructure   &fact_struct
00162     ,FactorizationNonzeros          *fact_nonzeros
00163     );
00165   void solve(
00166     const FactorizationStructure    &fact_struct
00167     ,const FactorizationNonzeros    &fact_nonzeros
00168     ,bool                           transp
00169     ,int                            n
00170     ,int                            nrhs
00171     ,double                         rhs[]
00172     ,int                            ldrhs
00173     ) const;
00174 
00176 
00177 private:
00178 
00180   void copy_basis_nonzeros(
00181     int                             m_orig
00182     ,int                            n_orig
00183     ,int                            nz_orig
00184     ,const double                   a_orig_val[]
00185     ,const int                      a_orig_row_i[]
00186     ,const int                      a_orig_col_ptr[]
00187     ,const int                      a_orig_perm_r[]
00188     ,const int                      a_orig_perm_c[]
00189     ,const int                      rank
00190     ,double                         b_val[]
00191     ,int                            b_row_i[]
00192     ,int                            b_col_ptr[]
00193     ,int                            *b_nz
00194     ) const;
00195 
00196 }; // end class SuperLUSolver
00197 
00198 //
00199 // SuperLUSolver
00200 //
00201 
00202 Teuchos::RCP<SuperLUSolver>
00203 SuperLUSolver::create_solver()
00204 {
00205   return Teuchos::rcp(new SuperLUSolverImpl());
00206 }
00207 
00208 Teuchos::RCP<SuperLUSolver::FactorizationStructure>
00209 SuperLUSolver::create_fact_struct()
00210 {
00211   return Teuchos::rcp(new SuperLUSolverImpl::FactorizationStructureImpl());
00212 }
00213 
00214 Teuchos::RCP<SuperLUSolver::FactorizationNonzeros>
00215 SuperLUSolver::create_fact_nonzeros()
00216 {
00217   return Teuchos::rcp(new SuperLUSolverImpl::FactorizationNonzerosImpl());
00218 }
00219 
00220 //
00221 // SuperLUSolverImp
00222 //
00223 
00224 // Overridden from SuperLUSolver
00225 
00226 void SuperLUSolverImpl::analyze_and_factor(
00227   int                         m
00228   ,int                        n
00229   ,int                        nz
00230   ,const double               a_val[]
00231   ,const int                  a_row_i[]
00232   ,const int                  a_col_ptr[]
00233   ,FactorizationStructure     *fact_struct
00234   ,FactorizationNonzeros      *fact_nonzeros
00235   ,int                        perm_r[]
00236   ,int                        perm_c[]
00237   ,int                        *rank
00238   )
00239 {
00240   using Teuchos::dyn_cast;
00241   using Teuchos::Workspace;
00242   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00243 
00244   FactorizationStructureImpl
00245     &fs = dyn_cast<FactorizationStructureImpl>(*fact_struct);
00246   FactorizationNonzerosImpl
00247     &fn = dyn_cast<FactorizationNonzerosImpl>(*fact_nonzeros);
00248 
00249   char refact[] = "N";
00250 
00251   // Resize storage.
00252   // Note: if this function was called recursively for m>n on the last call
00253   // then m_orig, n_orig etc. will already be set and should not be
00254   // disturbed.
00255   fs.rank_ = n; // Assume this for now
00256   fs.nz_   = nz;
00257   fs.perm_r_.resize(m);
00258   fs.perm_c_.resize(n);
00259   fs.etree_.resize(n);
00260 
00261     // Create matrix A in the format expected by SuperLU
00262   SuperMatrix A;
00263   dCreate_CompCol_Matrix(
00264     &A, m, n, nz
00265     ,const_cast<double*>(a_val)
00266     ,const_cast<int*>(a_row_i)
00267     ,const_cast<int*>(a_col_ptr)
00268     ,NC, D_, GE
00269     );
00270 
00271   // Get the columm permutations
00272   int permc_spec = 0; // ToDo: Make this an external parameter
00273   get_perm_c(permc_spec, &A, &fs.perm_c_[0]);
00274 
00275   // Permute the columns of the matrix
00276   SuperMatrix AC;
00277   sp_preorder(refact,&A,&fs.perm_c_[0],&fs.etree_[0],&AC);
00278 
00279   int info = -1;
00280   dgstrf(
00281     refact
00282     ,&AC  
00283     ,1.0    /* diag_pivot_thresh */
00284     ,0.0    /* drop_tol */
00285     ,local_relax
00286     ,local_panel_size
00287     ,&fs.etree_[0]
00288     ,NULL   /* work */
00289     ,0      /* lwork */
00290     ,&fs.perm_r_[0]
00291     ,&fs.perm_c_[0]
00292     ,&fn.L_
00293     ,&fn.U_
00294     ,&info
00295     );
00296 
00297   TEUCHOS_TEST_FOR_EXCEPTION(
00298     info != 0, std::runtime_error
00299     ,"SuperLUSolverImpl::analyze_and_factor(...): Error, dgstrf(...) returned info = " << info
00300     );
00301 
00302   std::copy( &fs.perm_r_[0], &fs.perm_r_[0] + m, perm_r );
00303   std::copy( &fs.perm_c_[0], &fs.perm_c_[0] + n, perm_c );
00304   *rank = n; // We must assume this until I can figure out a way to do better!
00305 
00306   if(m > n) {
00307     // Now we must refactor the basis by only passing in the elements for the basis
00308     // determined by SuperLU.  This is wasteful but it is the easiest thing to do
00309     // for now.
00310     fs.rank_        = *rank;
00311     fs.m_orig_      = m;
00312     fs.n_orig_      = n;
00313     fs.nz_orig_     = nz;
00314     fs.perm_r_orig_ = fs.perm_r_;
00315     fs.perm_c_orig_ = fs.perm_c_;
00316     // Copy the nonzeros for the sqare factor into new storage
00317     Workspace<double>       b_val(wss,nz);
00318     Workspace<int>          b_row_i(wss,nz);
00319     Workspace<int>          b_col_ptr(wss,n+1);
00320     copy_basis_nonzeros(
00321       m,n,nz,a_val,a_row_i,a_col_ptr
00322       ,&fs.perm_r_orig_[0],&fs.perm_c_orig_[0],fs.rank_
00323       ,&b_val[0],&b_row_i[0],&b_col_ptr[0]
00324       ,&fs.nz_
00325       );
00326     // Analyze and factor the new matrix
00327     int b_rank = -1;
00328     analyze_and_factor(
00329       fs.rank_, fs.rank_, fs.nz_, &b_val[0], &b_row_i[0], &b_col_ptr[0]
00330       ,fact_struct, fact_nonzeros
00331       ,&fs.perm_r_[0], &fs.perm_c_[0], &b_rank
00332       );
00333     TEUCHOS_TEST_FOR_EXCEPTION(
00334       (b_rank != *rank), std::runtime_error
00335       ,"SuperLUSolverImpl::analyze_and_factor(...): Error, the rank determined by "
00336       "the factorization of the rectangular " << m << " x " << n << " matrix of "
00337       << (*rank) << " is not the same as the refactorization of the basis returned as "
00338       << b_rank << "!"
00339       );
00340   }
00341   else {
00342     fs.m_orig_  = m;
00343     fs.n_orig_  = n;
00344     fs.nz_orig_ = nz;
00345   }
00346 }
00347 
00348 void SuperLUSolverImpl::factor(
00349   int                             m
00350   ,int                            n
00351   ,int                            nz
00352   ,const double                   a_val[]
00353   ,const int                      a_row_i[]
00354   ,const int                      a_col_ptr[]
00355   ,const FactorizationStructure   &fact_struct
00356   ,FactorizationNonzeros          *fact_nonzeros
00357   )
00358 {
00359   using Teuchos::dyn_cast;
00360   using Teuchos::Workspace;
00361   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00362 
00363   const FactorizationStructureImpl
00364     &fs = dyn_cast<const FactorizationStructureImpl>(fact_struct);
00365   FactorizationNonzerosImpl
00366     &fn = dyn_cast<FactorizationNonzerosImpl>(*fact_nonzeros);
00367 
00368   char refact[] = "Y";
00369 
00370   // Copy the nonzeros for the sqare factor into new storage
00371   Workspace<double>       b_val(wss,fs.nz_);
00372   Workspace<int>          b_row_i(wss,fs.nz_);
00373   Workspace<int>          b_col_ptr(wss,fs.rank_+1);
00374   if(fs.m_orig_ > fs.n_orig_) {
00375     int b_nz = -1;
00376     copy_basis_nonzeros(
00377       m,n,nz,a_val,a_row_i,a_col_ptr
00378       ,&cva(fs.perm_r_orig_)[0],&cva(fs.perm_c_orig_)[0],fs.rank_
00379       ,&b_val[0],&b_row_i[0],&b_col_ptr[0]
00380       ,&b_nz
00381       );
00382     TEUCHOS_TEST_FOR_EXCEPTION(
00383       (b_nz != fs.nz_), std::runtime_error
00384       ,"SuperLUSolverImpl::factor(...): Error!"
00385       );
00386   }
00387   else {
00388     std::copy( a_val,     a_val     + nz,  &b_val[0]     );
00389     std::copy( a_row_i,   a_row_i   + nz,  &b_row_i[0]   );
00390     std::copy( a_col_ptr, a_col_ptr + n+1, &b_col_ptr[0] );
00391   }
00392 
00393     // Create matrix A in the format expected by SuperLU
00394   SuperMatrix A;
00395   dCreate_CompCol_Matrix(
00396     &A, fs.rank_, fs.rank_, fs.nz_
00397     ,&b_val[0]
00398     ,&b_row_i[0]
00399     ,&b_col_ptr[0]
00400     ,NC, D_, GE
00401     );
00402 
00403   // Permute the columns
00404   SuperMatrix AC;
00405   sp_preorder(
00406     refact,&A
00407     ,&cva(fs.perm_c_)[0]
00408     ,&cva(fs.etree_)[0]
00409     ,&AC
00410     );
00411 
00412   int info = -1;
00413   dgstrf(
00414     refact
00415     ,&AC  
00416     ,1.0    /* diag_pivot_thresh */
00417     ,0.0    /* drop_tol */
00418     ,local_relax
00419     ,local_panel_size
00420     ,const_cast<int*>(&cva(fs.etree_)[0])
00421     ,NULL   /* work */
00422     ,0      /* lwork */
00423     ,&cva(fs.perm_r_)[0]
00424     ,&cva(fs.perm_c_)[0]
00425     ,&fn.L_
00426     ,&fn.U_
00427     ,&info
00428     );
00429 
00430   TEUCHOS_TEST_FOR_EXCEPTION(
00431     info != 0, std::runtime_error
00432     ,"SuperLUSolverImpl::factor(...): Error, dgstrf(...) returned info = " << info
00433     );
00434 
00435 }
00436 
00437 void SuperLUSolverImpl::solve(
00438   const FactorizationStructure    &fact_struct
00439   ,const FactorizationNonzeros    &fact_nonzeros
00440   ,bool                           transp
00441   ,int                            n
00442   ,int                            nrhs
00443   ,double                         rhs[]
00444   ,int                            ldrhs
00445   ) const
00446 {
00447 
00448   using Teuchos::dyn_cast;
00449 
00450   const FactorizationStructureImpl
00451     &fs = dyn_cast<const FactorizationStructureImpl>(fact_struct);
00452   const FactorizationNonzerosImpl
00453     &fn = dyn_cast<const FactorizationNonzerosImpl>(fact_nonzeros);
00454 
00455   TEUCHOS_TEST_FOR_EXCEPTION(
00456     n != fs.rank_, std::runtime_error
00457     ,"SuperLUSolverImpl::solve(...): Error, the dimmensions n = " << n << " and fs.rank = " << fs.rank_
00458     << " do not match up!"
00459     );
00460 
00461   SuperMatrix B;
00462     dCreate_Dense_Matrix(&B, n, nrhs, rhs, ldrhs, DN, D_, GE);
00463 
00464   char transc[1];
00465   transc[0] = ( transp ? 'T' : 'N' );
00466 
00467   int info = -1;
00468     dgstrs(
00469     transc
00470     ,const_cast<SuperMatrix*>(&fn.L_)
00471     ,const_cast<SuperMatrix*>(&fn.U_)
00472     ,&cva(fs.perm_r_)[0]
00473     ,&cva(fs.perm_c_)[0]
00474     ,&B, &info
00475     );
00476 
00477   TEUCHOS_TEST_FOR_EXCEPTION(
00478     info != 0, std::runtime_error
00479     ,"SuperLUSolverImpl::solve(...): Error, dgssv(...) returned info = " << info
00480     );
00481 
00482 }
00483 
00484 // private
00485 
00486 void SuperLUSolverImpl::copy_basis_nonzeros(
00487   int                             m_orig
00488   ,int                            n_orig
00489   ,int                            nz_orig
00490   ,const double                   a_orig_val[]
00491   ,const int                      a_orig_row_i[]
00492   ,const int                      a_orig_col_ptr[]
00493   ,const int                      a_orig_perm_r[]
00494   ,const int                      a_orig_perm_c[]
00495   ,const int                      rank
00496   ,double                         b_val[]
00497   ,int                            b_row_i[]
00498   ,int                            b_col_ptr[]
00499   ,int                            *b_nz
00500   ) const
00501 {
00502   *b_nz = 0;
00503   b_col_ptr[0] = *b_nz;
00504   for( int j = 0; j < rank; ++j ) {
00505     const int col_start_k = a_orig_col_ptr[j];
00506     const int col_end_k   = a_orig_col_ptr[j+1];
00507     for( int k = col_start_k; k < col_end_k; ++k ) {
00508       const int i_orig = a_orig_row_i[k];
00509       if(i_orig < rank) {
00510         b_val[*b_nz]     = a_orig_val[k];
00511         b_row_i[*b_nz]   = a_orig_row_i[k];
00512         ++(*b_nz);
00513       }
00514     }
00515     b_col_ptr[j+1] = *b_nz;
00516   }
00517 }
00518 
00519 } // end namespace SuperLUPack
00520 
00521 #endif // SPARSE_SOLVER_PACK_USE_SUPERLU
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines