AbstractLinAlgPack_DirectSparseSolverImp.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 #include <assert.h>
00030 
00031 #include "AbstractLinAlgPack_DirectSparseSolverImp.hpp"
00032 #include "Teuchos_AbstractFactoryStd.hpp"
00033 #include "Teuchos_TestForException.hpp"
00034 #include "Teuchos_dyn_cast.hpp"
00035 
00036 namespace AbstractLinAlgPack {
00037 
00038 // /////////////////////////////////////////
00039 // DirectSparseSolverImp::BasisMatrixImp
00040 
00041 DirectSparseSolverImp::BasisMatrixImp::BasisMatrixImp()
00042   : dim_(0)
00043 {}
00044 
00045 DirectSparseSolverImp::BasisMatrixImp::BasisMatrixImp(
00046   size_type                      dim
00047   ,const fact_struc_ptr_t        &fact_struc
00048   ,const fact_nonzeros_ptr_t     &fact_nonzeros
00049   )
00050 {
00051   this->initialize(dim,fact_struc,fact_nonzeros);
00052 }
00053 
00054 void DirectSparseSolverImp::BasisMatrixImp::initialize(
00055   size_type                      dim
00056   ,const fact_struc_ptr_t        &fact_struc
00057   ,const fact_nonzeros_ptr_t     &fact_nonzeros
00058   )
00059 {
00060 #ifdef TEUCHOS_DEBUG
00061   const char msg_err[] = "DirectSparseSolverImp::BasisMatrixImp::initialize(...): Error!";
00062   TEST_FOR_EXCEPTION( dim < 0, std::logic_error, msg_err );
00063   TEST_FOR_EXCEPTION( fact_struc.get() == NULL, std::logic_error, msg_err );
00064   TEST_FOR_EXCEPTION( fact_nonzeros.get() == NULL, std::logic_error, msg_err );
00065 #endif
00066   dim_            = dim;
00067   fact_struc_     = fact_struc;
00068   fact_nonzeros_  = fact_nonzeros;
00069   vec_space_.initialize(dim);
00070 }
00071 
00072 void DirectSparseSolverImp::BasisMatrixImp::set_uninitialized()
00073 {
00074   dim_            = 0;
00075   fact_struc_     = Teuchos::null;
00076   fact_nonzeros_  = Teuchos::null;
00077   vec_space_.initialize(0);
00078 }
00079 
00080 const DirectSparseSolverImp::BasisMatrixImp::fact_nonzeros_ptr_t&
00081 DirectSparseSolverImp::BasisMatrixImp::get_fact_nonzeros() const
00082 {
00083   return fact_nonzeros_;
00084 }
00085 
00086 // Overridden from MatrixBase
00087 
00088 const VectorSpace& DirectSparseSolverImp::BasisMatrixImp::space_cols() const
00089 {
00090   return vec_space_;
00091 }
00092 
00093 const VectorSpace& DirectSparseSolverImp::BasisMatrixImp::space_rows() const
00094 {
00095   return vec_space_;
00096 }
00097 
00098 size_type DirectSparseSolverImp::BasisMatrixImp::rows() const
00099 {
00100   return dim_;
00101 }
00102 
00103 size_type DirectSparseSolverImp::BasisMatrixImp::cols() const
00104 {
00105   return dim_;
00106 }
00107 
00108 // Overridden from MatrixNonsinguar
00109 
00110 MatrixNonsing::mat_mns_mut_ptr_t
00111 DirectSparseSolverImp::BasisMatrixImp::clone_mns()
00112 {
00113   namespace rcp = MemMngPack;
00114   Teuchos::RCP<BasisMatrixImp> bm = this->create_matrix();
00115   // A shallow copy is okay if the educated client DirectSparseSolverImp is careful!
00116   bm->initialize(dim_,fact_struc_,fact_nonzeros_);
00117   return bm;
00118 }
00119 
00120 // Overridden from BasisMatrix
00121 
00122 const DirectSparseSolver::BasisMatrix::fact_struc_ptr_t&
00123 DirectSparseSolverImp::BasisMatrixImp::get_fact_struc() const
00124 {
00125   return fact_struc_;
00126 }
00127 
00128 // //////////////////////////
00129 // DirectSparseSolverImp
00130 
00131 // Overridden from DirectSparseSolver
00132 
00133 void DirectSparseSolverImp::analyze_and_factor(
00134   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00135   ,DenseLinAlgPack::IVector                            *row_perm
00136   ,DenseLinAlgPack::IVector                            *col_perm
00137   ,size_type                                      *rank
00138   ,BasisMatrix                                    *basis_matrix
00139   ,std::ostream                                   *out
00140   )
00141 {
00142   namespace mmp = MemMngPack;
00143   using Teuchos::dyn_cast;
00144 #ifdef TEUCHOS_DEBUG
00145   const char msg_err[] = "DirectSparseSolverImp::analyze_and_factor(...): Error!";
00146   TEST_FOR_EXCEPTION( row_perm == NULL, std::logic_error, msg_err );
00147   TEST_FOR_EXCEPTION( col_perm == NULL, std::logic_error, msg_err );
00148   TEST_FOR_EXCEPTION( rank == NULL, std::logic_error, msg_err );
00149   TEST_FOR_EXCEPTION( basis_matrix == NULL, std::logic_error, msg_err );
00150 #endif
00151   BasisMatrixImp
00152     &basis_matrix_imp = dyn_cast<BasisMatrixImp>(*basis_matrix);
00153   // Get reference to factorization structure object or determine that we have to
00154   // allocate a new factorization structure object.
00155   const BasisMatrix::fact_struc_ptr_t        &bm_fact_struc   = basis_matrix_imp.get_fact_struc();
00156   const BasisMatrix::fact_struc_ptr_t        &this_fact_struc = this->get_fact_struc();
00157   BasisMatrix::fact_struc_ptr_t              fact_struc;
00158   bool alloc_new_fact_struc    = false;
00159   if( bm_fact_struc.count() == 1 )
00160     fact_struc = bm_fact_struc;
00161   else if( this_fact_struc.count() == 1 )
00162     fact_struc = this_fact_struc;
00163   else if( bm_fact_struc.get() == this_fact_struc.get() && bm_fact_struc.count() == 2 )
00164     fact_struc = bm_fact_struc;
00165   else
00166     alloc_new_fact_struc = true;
00167   if( alloc_new_fact_struc )
00168     fact_struc = this->create_fact_struc();
00169   // Get references to factorization nonzeros object or allocate a new factorization nonzeros object.
00170   const BasisMatrixImp::fact_nonzeros_ptr_t    &bm_fact_nonzeros  = basis_matrix_imp.get_fact_nonzeros();
00171   BasisMatrixImp::fact_nonzeros_ptr_t          fact_nonzeros;
00172   if( bm_fact_nonzeros.count() == 1 )
00173     fact_nonzeros = bm_fact_nonzeros;
00174   else
00175     fact_nonzeros = this->create_fact_nonzeros();
00176   // Now ask the subclass to do the work
00177   this->imp_analyze_and_factor(
00178     A,fact_struc.get(),fact_nonzeros.get(),row_perm,col_perm,rank,out
00179     );
00180   // Setup the basis matrix
00181   basis_matrix_imp.initialize(*rank,fact_struc,fact_nonzeros);
00182   // Remember rank and factorization structure
00183   rank_       = *rank;
00184   fact_struc_ = fact_struc;
00185 }
00186 
00187 void DirectSparseSolverImp::factor(
00188   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00189   ,BasisMatrix                                    *basis_matrix
00190   ,const BasisMatrix::fact_struc_ptr_t            &fact_struc_in
00191   ,std::ostream                                   *out
00192   )
00193 {
00194   namespace mmp = MemMngPack;
00195   using Teuchos::dyn_cast;
00196 #ifdef TEUCHOS_DEBUG
00197   const char msg_err[] = "DirectSparseSolverImp::analyze_and_factor(...): Error!";
00198   // ToDo: Validate that A is compatible!
00199   TEST_FOR_EXCEPTION( basis_matrix == NULL, std::logic_error, msg_err );
00200 #endif
00201   BasisMatrixImp
00202     &basis_matrix_imp = dyn_cast<BasisMatrixImp>(*basis_matrix);
00203   // Get reference to factorization structure object
00204   const BasisMatrix::fact_struc_ptr_t        &this_fact_struc = this->get_fact_struc();
00205   BasisMatrix::fact_struc_ptr_t              fact_struc;
00206 #ifdef TEUCHOS_DEBUG
00207   TEST_FOR_EXCEPTION(
00208     fact_struc_in.get() == NULL && this_fact_struc.get() == NULL
00209     ,std::logic_error
00210     ,msg_err );
00211 #endif
00212   if( fact_struc_in.get() != NULL )
00213     fact_struc = fact_struc_in;
00214   else
00215     fact_struc = this_fact_struc;
00216   // Get references to factorization nonzeros object or allocate a new factorization nonzeros object.
00217   const BasisMatrixImp::fact_nonzeros_ptr_t    &bm_fact_nonzeros  = basis_matrix_imp.get_fact_nonzeros();
00218   BasisMatrixImp::fact_nonzeros_ptr_t          fact_nonzeros;
00219   if( bm_fact_nonzeros.count() == 1 )
00220     fact_nonzeros = bm_fact_nonzeros;
00221   else
00222     fact_nonzeros = this->create_fact_nonzeros();
00223   // Now ask the subclass to do the work
00224   this->imp_factor(A,*fact_struc,fact_nonzeros.get(),out);
00225   // Setup the basis matrix
00226   basis_matrix_imp.initialize(rank_,fact_struc,fact_nonzeros);
00227 }
00228 
00229 const DirectSparseSolver::BasisMatrix::fact_struc_ptr_t&
00230 DirectSparseSolverImp::get_fact_struc() const
00231 {
00232   return fact_struc_;
00233 }
00234 
00235 void DirectSparseSolverImp::set_uninitialized()
00236 {
00237   fact_struc_ = Teuchos::null;
00238 }
00239 
00240 } // end namespace AbstractLinAlgPack 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:09:15 2011 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.6.3