MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_DirectSparseSolverImp.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
46 #include "Teuchos_Assert.hpp"
47 #include "Teuchos_dyn_cast.hpp"
48 
49 namespace AbstractLinAlgPack {
50 
51 // /////////////////////////////////////////
52 // DirectSparseSolverImp::BasisMatrixImp
53 
55  : dim_(0)
56 {}
57 
59  size_type dim
60  ,const fact_struc_ptr_t &fact_struc
61  ,const fact_nonzeros_ptr_t &fact_nonzeros
62  )
63 {
64  this->initialize(dim,fact_struc,fact_nonzeros);
65 }
66 
68  size_type dim
69  ,const fact_struc_ptr_t &fact_struc
70  ,const fact_nonzeros_ptr_t &fact_nonzeros
71  )
72 {
73 #ifdef TEUCHOS_DEBUG
74  const char msg_err[] = "DirectSparseSolverImp::BasisMatrixImp::initialize(...): Error!";
75  TEUCHOS_TEST_FOR_EXCEPTION( dim < 0, std::logic_error, msg_err );
76  TEUCHOS_TEST_FOR_EXCEPTION( fact_struc.get() == NULL, std::logic_error, msg_err );
77  TEUCHOS_TEST_FOR_EXCEPTION( fact_nonzeros.get() == NULL, std::logic_error, msg_err );
78 #endif
79  dim_ = dim;
80  fact_struc_ = fact_struc;
81  fact_nonzeros_ = fact_nonzeros;
82  vec_space_.initialize(dim);
83 }
84 
86 {
87  dim_ = 0;
88  fact_struc_ = Teuchos::null;
89  fact_nonzeros_ = Teuchos::null;
90  vec_space_.initialize(0);
91 }
92 
95 {
96  return fact_nonzeros_;
97 }
98 
99 // Overridden from MatrixBase
100 
102 {
103  return vec_space_;
104 }
105 
107 {
108  return vec_space_;
109 }
110 
112 {
113  return dim_;
114 }
115 
117 {
118  return dim_;
119 }
120 
121 // Overridden from MatrixNonsinguar
122 
123 MatrixNonsing::mat_mns_mut_ptr_t
125 {
126  namespace rcp = MemMngPack;
127  Teuchos::RCP<BasisMatrixImp> bm = this->create_matrix();
128  // A shallow copy is okay if the educated client DirectSparseSolverImp is careful!
129  bm->initialize(dim_,fact_struc_,fact_nonzeros_);
130  return bm;
131 }
132 
133 // Overridden from BasisMatrix
134 
137 {
138  return fact_struc_;
139 }
140 
141 // //////////////////////////
142 // DirectSparseSolverImp
143 
144 // Overridden from DirectSparseSolver
145 
148  ,DenseLinAlgPack::IVector *row_perm
149  ,DenseLinAlgPack::IVector *col_perm
150  ,size_type *rank
151  ,BasisMatrix *basis_matrix
152  ,std::ostream *out
153  )
154 {
155  namespace mmp = MemMngPack;
156  using Teuchos::dyn_cast;
157 #ifdef TEUCHOS_DEBUG
158  const char msg_err[] = "DirectSparseSolverImp::analyze_and_factor(...): Error!";
159  TEUCHOS_TEST_FOR_EXCEPTION( row_perm == NULL, std::logic_error, msg_err );
160  TEUCHOS_TEST_FOR_EXCEPTION( col_perm == NULL, std::logic_error, msg_err );
161  TEUCHOS_TEST_FOR_EXCEPTION( rank == NULL, std::logic_error, msg_err );
162  TEUCHOS_TEST_FOR_EXCEPTION( basis_matrix == NULL, std::logic_error, msg_err );
163 #endif
165  &basis_matrix_imp = dyn_cast<BasisMatrixImp>(*basis_matrix);
166  // Get reference to factorization structure object or determine that we have to
167  // allocate a new factorization structure object.
168  const BasisMatrix::fact_struc_ptr_t &bm_fact_struc = basis_matrix_imp.get_fact_struc();
169  const BasisMatrix::fact_struc_ptr_t &this_fact_struc = this->get_fact_struc();
171  bool alloc_new_fact_struc = false;
172  if( bm_fact_struc.total_count() == 1 )
173  fact_struc = bm_fact_struc;
174  else if( this_fact_struc.total_count() == 1 )
175  fact_struc = this_fact_struc;
176  else if( bm_fact_struc.get() == this_fact_struc.get() && bm_fact_struc.total_count() == 2 )
177  fact_struc = bm_fact_struc;
178  else
179  alloc_new_fact_struc = true;
180  if( alloc_new_fact_struc )
181  fact_struc = this->create_fact_struc();
182  // Get references to factorization nonzeros object or allocate a new factorization nonzeros object.
183  const BasisMatrixImp::fact_nonzeros_ptr_t &bm_fact_nonzeros = basis_matrix_imp.get_fact_nonzeros();
185  if( bm_fact_nonzeros.total_count() == 1 )
186  fact_nonzeros = bm_fact_nonzeros;
187  else
188  fact_nonzeros = this->create_fact_nonzeros();
189  // Now ask the subclass to do the work
191  A,fact_struc.get(),fact_nonzeros.get(),row_perm,col_perm,rank,out
192  );
193  // Setup the basis matrix
194  basis_matrix_imp.initialize(*rank,fact_struc,fact_nonzeros);
195  // Remember rank and factorization structure
196  rank_ = *rank;
197  fact_struc_ = fact_struc;
198 }
199 
202  ,BasisMatrix *basis_matrix
203  ,const BasisMatrix::fact_struc_ptr_t &fact_struc_in
204  ,std::ostream *out
205  )
206 {
207  namespace mmp = MemMngPack;
208  using Teuchos::dyn_cast;
209 #ifdef TEUCHOS_DEBUG
210  const char msg_err[] = "DirectSparseSolverImp::analyze_and_factor(...): Error!";
211  // ToDo: Validate that A is compatible!
212  TEUCHOS_TEST_FOR_EXCEPTION( basis_matrix == NULL, std::logic_error, msg_err );
213 #endif
215  &basis_matrix_imp = dyn_cast<BasisMatrixImp>(*basis_matrix);
216  // Get reference to factorization structure object
217  const BasisMatrix::fact_struc_ptr_t &this_fact_struc = this->get_fact_struc();
219 #ifdef TEUCHOS_DEBUG
221  fact_struc_in.get() == NULL && this_fact_struc.get() == NULL
222  ,std::logic_error
223  ,msg_err );
224 #endif
225  if( fact_struc_in.get() != NULL )
226  fact_struc = fact_struc_in;
227  else
228  fact_struc = this_fact_struc;
229  // Get references to factorization nonzeros object or allocate a new factorization nonzeros object.
230  const BasisMatrixImp::fact_nonzeros_ptr_t &bm_fact_nonzeros = basis_matrix_imp.get_fact_nonzeros();
232  if( bm_fact_nonzeros.total_count() == 1 )
233  fact_nonzeros = bm_fact_nonzeros;
234  else
235  fact_nonzeros = this->create_fact_nonzeros();
236  // Now ask the subclass to do the work
237  this->imp_factor(A,*fact_struc,fact_nonzeros.get(),out);
238  // Setup the basis matrix
239  basis_matrix_imp.initialize(rank_,fact_struc,fact_nonzeros);
240 }
241 
244 {
245  return fact_struc_;
246 }
247 
249 {
250  fact_struc_ = Teuchos::null;
251 }
252 
253 } // end namespace AbstractLinAlgPack
const BasisMatrix::fact_struc_ptr_t & get_fact_struc() const
Abstract class for objects that represent the factorized matrix and can be used to solve for differen...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T * get() const
virtual const fact_nonzeros_ptr_t & get_fact_nonzeros() const
Return a reference to a smart pointer to the object that represents the factorization nonzeros...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual const Teuchos::RCP< FactorizationStructure > create_fact_struc() const =0
Create a new, uninitialized FactorizationStructure object.
virtual const Teuchos::RCP< FactorizationNonzeros > create_fact_nonzeros() const =0
Create a new, uninitialized FactorizationNonzeros object.
void factor(const AbstractLinAlgPack::MatrixConvertToSparse &A, BasisMatrix *basis_matrix, const BasisMatrix::fact_struc_ptr_t &fact_struc, std::ostream *out)
Abstract interface for objects that represent a space for mutable coordinate vectors.
virtual void imp_analyze_and_factor(const AbstractLinAlgPack::MatrixConvertToSparse &A, FactorizationStructure *fact_struc, FactorizationNonzeros *fact_nonzeros, DenseLinAlgPack::IVector *row_perm, DenseLinAlgPack::IVector *col_perm, size_type *rank, std::ostream *out=NULL)=0
Called to implement the analyze_and_factor() without having to worry about memory mangagment details...
T_To & dyn_cast(T_From &from)
Mix-in interface for extracing explicit elements from a sparse matrix in one of several Fortran compa...
virtual void initialize(size_type dim, const fact_struc_ptr_t &fact_struc, const fact_nonzeros_ptr_t &fact_nonzeros)
Initialize given initialized factorization structure and factorization nonzeros objects.
std::ostream * out
void analyze_and_factor(const AbstractLinAlgPack::MatrixConvertToSparse &A, DenseLinAlgPack::IVector *row_perm, DenseLinAlgPack::IVector *col_perm, size_type *rank, BasisMatrix *basis_matrix, std::ostream *out)
Implementation node subclass that combines factorization structure and factorization nonzeros into a ...
virtual void imp_factor(const AbstractLinAlgPack::MatrixConvertToSparse &A, const FactorizationStructure &fact_struc, FactorizationNonzeros *fact_nonzeros, std::ostream *out=NULL)=0
Called to implement the analyze_and_factor() without having to worry about memory mangagment details...
int total_count() const