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