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