AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_DirectSparseSolverDense.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 // 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 #include <assert.h>
00043 
00044 #include <fstream>
00045 #include <algorithm>
00046 
00047 #include "AbstractLinAlgPack_DirectSparseSolverDense.hpp"
00048 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00049 #include "DenseLinAlgLAPack.hpp"
00050 #include "DenseLinAlgPack_PermVecMat.hpp"
00051 #include "Teuchos_AbstractFactoryStd.hpp"
00052 #include "Teuchos_Assert.hpp"
00053 #include "Teuchos_Workspace.hpp"
00054 #include "Teuchos_dyn_cast.hpp"
00055 
00056 namespace {
00057 
00058 // My swap function
00059 template<class T>
00060 inline
00061 void my_swap( T* v1, T* v2 )
00062 {
00063   T tmp = *v1;
00064   *v1 = *v2;
00065   *v2 = tmp;
00066 }
00067 
00068 // A cast to const is needed because the standard does not return a reference from
00069 // valarray<>::operator[]() const.
00070 template <class T>
00071 std::valarray<T>& cva(const std::valarray<T>& va )
00072 {
00073   return const_cast<std::valarray<T>&>(va);
00074 }
00075 
00076 } // end namespace
00077 
00078 namespace AbstractLinAlgPack {
00079 
00080 //
00081 // Implementation of DirectSparseSolver(Imp) interface using dense LAPACK routines.
00082 //
00083 
00084 // //////////////////////////////////////////////////
00085 // DirectSparseSolverDense::BasisMatrixDense
00086 
00087 // Overridden from BasisMatrixImp
00088 
00089 Teuchos::RCP<DirectSparseSolverImp::BasisMatrixImp>
00090 DirectSparseSolverDense::BasisMatrixDense::create_matrix() const
00091 {
00092   return Teuchos::rcp(new BasisMatrixDense);
00093 }
00094 
00095 void DirectSparseSolverDense::BasisMatrixDense::V_InvMtV(
00096   VectorMutable* y, BLAS_Cpp::Transp M_trans, const Vector& x
00097   ) const 
00098 {
00099   using Teuchos::dyn_cast;
00100   using Teuchos::Workspace;
00101   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00102   size_type k;
00103 
00104   // Get concrete objects
00105   const FactorizationStructureDense
00106     &fs = dyn_cast<const FactorizationStructureDense>(*this->get_fact_struc());
00107   const FactorizationNonzerosDense
00108     &fn = dyn_cast<const FactorizationNonzerosDense>(*this->get_fact_nonzeros());
00109 
00110   VectorDenseMutableEncap  yd(*y);
00111   VectorDenseEncap         xd(x);
00112 
00113   TEUCHOS_TEST_FOR_EXCEPTION(
00114     yd().dim() != xd().dim(), std::invalid_argument
00115     ,"DirectSparseSolverDense::BasisMatrixDense::V_InvMtV(...) : Error, "
00116     " y.dim() = " << yd().dim() << " != x.dim() = " << xd().dim() << "!"
00117     );
00118 
00119   // Get temp storage for rhs and solution to communicate with xGESTRS
00120   Workspace<value_type>   B_store(wss,xd().dim());
00121   DMatrixSlice  B(&B_store[0],B_store.size(),B_store.size(),B_store.size(),1);
00122 
00123   //
00124   // Now we must permute the rhs or solution vectors based on our own
00125   // permutation fn.basis_perm_.
00126   //
00127   // xGETRF(...) factored the transpose of the basis matrix C' = Ct = P*L*U
00128   // where the permtuation P is stored in the array fn.basis_perm_ where
00129   //
00130   //    q = P * v
00131   //
00132   // is given by
00133   //
00134   //    q(i) = v(fn.basis_perm_(i)), for i = 1...n
00135   //
00136   // and q = P' * v is given by
00137   //
00138   //    q(fn.basis_perm_(i)) = v(i), for i = 1...n
00139   //
00140   // The system we are solving is therefore:
00141   //
00142   //   C * y = x   =>  U'*L'*P'*y = x
00143   //   
00144   //        for M_trans == no_trans
00145   //
00146   //   C'* y = x   =>  P*L*U*y = x   =>  L*U*y = P'*x 
00147   //
00148   //        for M_trans == trans
00149   // 
00150 
00151   // Copy rsh
00152   if( M_trans == BLAS_Cpp::trans && fn.rect_analyze_and_factor_ ) {
00153     // b = P'*x =
00154     DVectorSlice b = B.col(1);
00155 //    DenseLinAlgPack::inv_perm_ele(xd(),fn.basis_perm_,&b);
00156     DenseLinAlgPack::perm_ele(xd(),fn.basis_perm_,&b);
00157   }
00158   else {
00159     B.col(1) = xd();
00160   }
00161 
00162   // Solve
00163   DenseLinAlgLAPack::getrs(
00164     fn.LU_(1,fs.rank_,1,fs.rank_), &cva(fn.ipiv_)[0], BLAS_Cpp::trans_not(M_trans)
00165     ,&B
00166     );
00167 
00168   // Copy solution
00169   if( M_trans == BLAS_Cpp::no_trans  && fn.rect_analyze_and_factor_ ) {
00170     // y = P*b = P*(P'*y)
00171     const DVectorSlice b = B.col(1);
00172 //    DenseLinAlgPack::perm_ele(b,fn.basis_perm_,&yd());
00173     DenseLinAlgPack::inv_perm_ele(b,fn.basis_perm_,&yd());
00174   }
00175   else {
00176     yd() = B.col(1);
00177   }
00178 
00179 }
00180 
00181 // //////////////////////////////////////////////////
00182 // DirectSparseSolverDense::FactorizationStructureDense
00183 
00184 DirectSparseSolverDense::FactorizationStructureDense::FactorizationStructureDense()
00185 {}
00186 
00187 // //////////////////////////////////////////////////
00188 // DirectSparseSolverDense
00189 
00190 // Constructors/initializers
00191 
00192 DirectSparseSolverDense::DirectSparseSolverDense()
00193 {}
00194 
00195 // Overridden from DirectSparseSolver
00196 
00197 const DirectSparseSolver::basis_matrix_factory_ptr_t
00198 DirectSparseSolverDense::basis_matrix_factory() const
00199 {
00200   namespace mmp = MemMngPack;
00201   return Teuchos::rcp(new Teuchos::AbstractFactoryStd<BasisMatrix,BasisMatrixDense>());
00202 }
00203 
00204 void DirectSparseSolverDense::estimated_fillin_ratio(
00205   value_type estimated_fillin_ratio
00206   )
00207 {
00208   // We ignore this!
00209 }
00210 
00211 // Overridden from DirectSparseSolverImp
00212 
00213 const Teuchos::RCP<DirectSparseSolver::FactorizationStructure>
00214 DirectSparseSolverDense::create_fact_struc() const
00215 {
00216   return Teuchos::rcp(new FactorizationStructureDense);
00217 }
00218 
00219 const Teuchos::RCP<DirectSparseSolverImp::FactorizationNonzeros>
00220 DirectSparseSolverDense::create_fact_nonzeros() const
00221 {
00222   return Teuchos::rcp(new FactorizationNonzerosDense);
00223 }
00224 
00225 void DirectSparseSolverDense::imp_analyze_and_factor(
00226   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00227   ,FactorizationStructure                         *fact_struc
00228   ,FactorizationNonzeros                          *fact_nonzeros
00229   ,DenseLinAlgPack::IVector                            *row_perm
00230   ,DenseLinAlgPack::IVector                            *col_perm
00231   ,size_type                                      *rank
00232   ,std::ostream                                   *out
00233   )
00234 {
00235   namespace mmp = MemMngPack;
00236   using Teuchos::dyn_cast;
00237   typedef MatrixConvertToSparse MCTS;
00238   using Teuchos::Workspace;
00239   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00240 
00241   if(out)
00242     *out << "\nUsing LAPACK xGETRF to analyze and factor a new matrix ...\n";
00243 
00244   // Get the concrete factorization and nonzeros objects
00245   FactorizationStructureDense
00246     &fs = dyn_cast<FactorizationStructureDense>(*fact_struc);
00247   FactorizationNonzerosDense
00248     &fn = dyn_cast<FactorizationNonzerosDense>(*fact_nonzeros);
00249 
00250   // Get the dimensions of things.
00251   const index_type
00252     m = A.rows(),
00253     n = A.cols(),
00254     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00255 
00256   // Validate input
00257   TEUCHOS_TEST_FOR_EXCEPTION(
00258     n <= 0 || m <= 0 || m > n, std::invalid_argument
00259     ,"DirectSparseSolverDense::imp_analyze_and_factor(...) : Error!" );
00260 
00261   // Extract the matrix in coordinate format
00262   Workspace<value_type>   a_val(wss,nz);
00263   Workspace<index_type>   a_row_i(wss,nz);
00264   Workspace<index_type>   a_col_j(wss,nz);
00265   A.coor_extract_nonzeros(
00266     MCTS::EXTRACT_FULL_MATRIX
00267     ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00268     ,nz
00269     ,&a_val[0]
00270     ,nz
00271     ,&a_row_i[0]
00272     ,&a_col_j[0]
00273     );
00274 
00275   //
00276   // Fill the matrix LU = A' so that xGETRF will pivot by row to find
00277   // the basis.
00278   //
00279   // Here we will form the factor of A' = P*L*U where L will be
00280   // a n x m upper trapizodial matrix containing the factor lower
00281   // triangular factor in LU(1:rank,1:rank) and junk below this.
00282   //
00283   // Note that xGETRF() pivots by row so if any dependent columns
00284   // are found they will always be the last few columns.
00285   //
00286 
00287   // Resize the storage
00288   fn.LU_.resize(n,m);
00289   fn.ipiv_.resize(n);
00290 
00291   // Add in the nonzero entires transposed (allows for multiple entries with same
00292   // row and column indexes).
00293   fn.LU_ = 0.0;
00294   for( size_type k = 0; k < nz; ++k )
00295     fn.LU_(a_col_j[k],a_row_i[k]) += a_val[k];
00296 
00297   //
00298   // Have xGETRF factor this matrix.
00299   //
00300 
00301   DenseLinAlgLAPack::getrf( &fn.LU_(), &fn.ipiv_[0], &fs.rank_ );
00302 
00303   // Remember the dimensions
00304   fs.m_  = m;
00305   fs.n_  = n;
00306   fs.nz_ = nz;
00307 
00308   //
00309   // At this point it is important to understand exactly what
00310   // ipiv() represents.  Each entry in ipiv(k) represents a row
00311   // interchange A(k) <=> A(ipiv(k)).  Therefore, we have to
00312   // do the same row interchanges to the identity permutation
00313   // of col_perm to form the permutation array expected by the
00314   // DSS interface.
00315   //
00316 
00317   // Form fs.col_perm_
00318   fs.col_perm_.resize(n);
00319   DenseLinAlgPack::identity_perm(&fs.col_perm_);
00320   Workspace<index_type> col_perm_unsorted(wss,fs.rank_);
00321   if( m == n && n == fs.rank_ ) {
00322     // Leave fs.col_perm_ as identity
00323     fn.rect_analyze_and_factor_ = false;
00324   }
00325   else {
00326     fn.rect_analyze_and_factor_ = true;
00327     // Permute fs.col_perm_ and save them in col_perm_unsorted
00328     for( index_type k = 1; k <= fs.rank_; ++k ) {
00329       my_swap( &fs.col_perm_(k), &fs.col_perm_(fn.ipiv_[k-1]) );
00330       col_perm_unsorted[k-1] = fs.col_perm_(k);
00331     }
00332     // Sort the basis selection
00333     std::sort(&(fs.col_perm_)[0]           , &(fs.col_perm_)[0] + fs.rank_ );
00334     std::sort(&(fs.col_perm_)[0] + fs.rank_, &(fs.col_perm_)[0] + n        );
00335   }
00336 
00337   // Form the inverse permutation
00338   fs.inv_col_perm_.resize(n);
00339   DenseLinAlgPack::inv_perm( fs.col_perm_, &fs.inv_col_perm_ );
00340 
00341   if( !(m == n && n == fs.rank_) ) {
00342     // Form fn.basis_perm_ and set fs.ipiv_ to identity
00343     fn.basis_perm_.resize(fs.rank_);
00344     for( size_type k = 1; k <= fs.rank_; ++k ) {
00345       fn.basis_perm_(k) = fs.inv_col_perm_(col_perm_unsorted[k-1]);
00346       fn.ipiv_[k-1] = k;
00347     }
00348   }
00349 
00350   // Copy the data to the output
00351   row_perm->resize(m);
00352   col_perm->resize(n);
00353   *rank = fs.rank_;
00354   DenseLinAlgPack::identity_perm(row_perm);
00355   *col_perm = fs.col_perm_;
00356 
00357 }
00358 
00359 void DirectSparseSolverDense::imp_factor(
00360   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00361   ,const FactorizationStructure                   &fact_struc
00362   ,FactorizationNonzeros                          *fact_nonzeros
00363   ,std::ostream                                   *out
00364   )
00365 {
00366   namespace mmp = MemMngPack;
00367   using Teuchos::dyn_cast;
00368   typedef MatrixConvertToSparse MCTS;
00369   using Teuchos::Workspace;
00370   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00371 
00372   if(out)
00373     *out << "\nUsing LAPACK xGETRF to refactor the basis matrix ...\n";
00374 
00375   // Get the concrete factorization and nonzeros objects
00376   const FactorizationStructureDense
00377     &fs = dyn_cast<const FactorizationStructureDense>(fact_struc);
00378   FactorizationNonzerosDense
00379     &fn = dyn_cast<FactorizationNonzerosDense>(*fact_nonzeros);
00380 
00381   // Get the dimensions of things.
00382   const index_type
00383     m = A.rows(),
00384     n = A.cols(),
00385     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00386 
00387   // Validate input
00388   TEUCHOS_TEST_FOR_EXCEPTION(
00389     (m != fs.m_ || n != fs.n_ || nz != fs.nz_), std::invalid_argument
00390     ,"DirectSparseSolverDense::imp_factor(...): Error!"
00391     );
00392 
00393   // Extract the matrix in coordinate format
00394   Workspace<value_type>   a_val(wss,nz);
00395   Workspace<index_type>   a_row_i(wss,nz);
00396   Workspace<index_type>   a_col_j(wss,nz);
00397   A.coor_extract_nonzeros(
00398     MCTS::EXTRACT_FULL_MATRIX
00399     ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00400     ,nz
00401     ,&a_val[0]
00402     ,nz
00403     ,&a_row_i[0]
00404     ,&a_col_j[0]
00405     );
00406   
00407   //
00408   // Fill the matrix LU = B so that xGETRF will pivot by row to find
00409   // the basis.  Here B is the basis matrix from A'.
00410   //
00411   // Here we will form the factor of B = P*L*U where L will be
00412   // a rank x rank lower triangular.
00413   //
00414 
00415   // Resize the storage
00416   fn.rect_analyze_and_factor_ = false;
00417   fn.LU_.resize(fs.rank_,fs.rank_);
00418   fn.ipiv_.resize(fs.rank_);
00419 
00420   // Copy only the basis entries (transposed)
00421   fn.LU_ = 0.0;
00422   for( size_type k = 0; k < nz; ++k ) {
00423     const index_type B_i = fs.inv_col_perm_(a_col_j[k]);
00424     const index_type B_j = a_row_i[k];
00425     if( B_i <= fs.rank_ && B_j <= fs.rank_ )
00426       fn.LU_(B_i,B_j) = a_val[k];
00427   }
00428 
00429   //
00430   // Have xGETRF factor this matrix.
00431   //
00432 
00433   FortranTypes::f_int B_rank = 0;
00434   DenseLinAlgLAPack::getrf( &fn.LU_(), &fn.ipiv_[0], &B_rank );
00435 
00436   TEUCHOS_TEST_FOR_EXCEPTION(
00437     B_rank != fs.rank_, FactorizationFailure
00438     ,"DirectSparseSolverDense::imp_factor(...): Error, the basis matrix is no "
00439     "longer full rank with B_rank = " << B_rank << " != fs.rank = " << fs.rank_ << "!"
00440     );
00441 
00442 }
00443 
00444 // private
00445 
00446 } // end namespace AbstractLinAlgPack 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends