MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbstractLinAlgPack_DirectSparseSolverDense.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 
44 #include <fstream>
45 #include <algorithm>
46 
49 #include "DenseLinAlgLAPack.hpp"
52 #include "Teuchos_Assert.hpp"
53 #include "Teuchos_Workspace.hpp"
54 #include "Teuchos_dyn_cast.hpp"
55 
56 namespace {
57 
58 // My swap function
59 template<class T>
60 inline
61 void my_swap( T* v1, T* v2 )
62 {
63  T tmp = *v1;
64  *v1 = *v2;
65  *v2 = tmp;
66 }
67 
68 // A cast to const is needed because the standard does not return a reference from
69 // valarray<>::operator[]() const.
70 template <class T>
71 std::valarray<T>& cva(const std::valarray<T>& va )
72 {
73  return const_cast<std::valarray<T>&>(va);
74 }
75 
76 } // end namespace
77 
78 namespace AbstractLinAlgPack {
79 
80 //
81 // Implementation of DirectSparseSolver(Imp) interface using dense LAPACK routines.
82 //
83 
84 // //////////////////////////////////////////////////
85 // DirectSparseSolverDense::BasisMatrixDense
86 
87 // Overridden from BasisMatrixImp
88 
91 {
92  return Teuchos::rcp(new BasisMatrixDense);
93 }
94 
96  VectorMutable* y, BLAS_Cpp::Transp M_trans, const Vector& x
97  ) const
98 {
99  using Teuchos::dyn_cast;
100  using Teuchos::Workspace;
102  size_type k;
103 
104  // Get concrete objects
106  &fs = dyn_cast<const FactorizationStructureDense>(*this->get_fact_struc());
108  &fn = dyn_cast<const FactorizationNonzerosDense>(*this->get_fact_nonzeros());
109 
111  VectorDenseEncap xd(x);
112 
114  yd().dim() != xd().dim(), std::invalid_argument
115  ,"DirectSparseSolverDense::BasisMatrixDense::V_InvMtV(...) : Error, "
116  " y.dim() = " << yd().dim() << " != x.dim() = " << xd().dim() << "!"
117  );
118 
119  // Get temp storage for rhs and solution to communicate with xGESTRS
120  Workspace<value_type> B_store(wss,xd().dim());
121  DMatrixSlice B(&B_store[0],B_store.size(),B_store.size(),B_store.size(),1);
122 
123  //
124  // Now we must permute the rhs or solution vectors based on our own
125  // permutation fn.basis_perm_.
126  //
127  // xGETRF(...) factored the transpose of the basis matrix C' = Ct = P*L*U
128  // where the permtuation P is stored in the array fn.basis_perm_ where
129  //
130  // q = P * v
131  //
132  // is given by
133  //
134  // q(i) = v(fn.basis_perm_(i)), for i = 1...n
135  //
136  // and q = P' * v is given by
137  //
138  // q(fn.basis_perm_(i)) = v(i), for i = 1...n
139  //
140  // The system we are solving is therefore:
141  //
142  // C * y = x => U'*L'*P'*y = x
143  //
144  // for M_trans == no_trans
145  //
146  // C'* y = x => P*L*U*y = x => L*U*y = P'*x
147  //
148  // for M_trans == trans
149  //
150 
151  // Copy rsh
152  if( M_trans == BLAS_Cpp::trans && fn.rect_analyze_and_factor_ ) {
153  // b = P'*x =
154  DVectorSlice b = B.col(1);
155 // DenseLinAlgPack::inv_perm_ele(xd(),fn.basis_perm_,&b);
156  DenseLinAlgPack::perm_ele(xd(),fn.basis_perm_,&b);
157  }
158  else {
159  B.col(1) = xd();
160  }
161 
162  // Solve
164  fn.LU_(1,fs.rank_,1,fs.rank_), &cva(fn.ipiv_)[0], BLAS_Cpp::trans_not(M_trans)
165  ,&B
166  );
167 
168  // Copy solution
169  if( M_trans == BLAS_Cpp::no_trans && fn.rect_analyze_and_factor_ ) {
170  // y = P*b = P*(P'*y)
171  const DVectorSlice b = B.col(1);
172 // DenseLinAlgPack::perm_ele(b,fn.basis_perm_,&yd());
173  DenseLinAlgPack::inv_perm_ele(b,fn.basis_perm_,&yd());
174  }
175  else {
176  yd() = B.col(1);
177  }
178 
179 }
180 
181 // //////////////////////////////////////////////////
182 // DirectSparseSolverDense::FactorizationStructureDense
183 
185 {}
186 
187 // //////////////////////////////////////////////////
188 // DirectSparseSolverDense
189 
190 // Constructors/initializers
191 
193 {}
194 
195 // Overridden from DirectSparseSolver
196 
199 {
200  namespace mmp = MemMngPack;
202 }
203 
206  )
207 {
208  // We ignore this!
209 }
210 
211 // Overridden from DirectSparseSolverImp
212 
215 {
217 }
218 
221 {
223 }
224 
228  ,FactorizationNonzeros *fact_nonzeros
229  ,DenseLinAlgPack::IVector *row_perm
230  ,DenseLinAlgPack::IVector *col_perm
231  ,size_type *rank
232  ,std::ostream *out
233  )
234 {
235  namespace mmp = MemMngPack;
236  using Teuchos::dyn_cast;
237  typedef MatrixConvertToSparse MCTS;
238  using Teuchos::Workspace;
240 
241  if(out)
242  *out << "\nUsing LAPACK xGETRF to analyze and factor a new matrix ...\n";
243 
244  // Get the concrete factorization and nonzeros objects
246  &fs = dyn_cast<FactorizationStructureDense>(*fact_struc);
248  &fn = dyn_cast<FactorizationNonzerosDense>(*fact_nonzeros);
249 
250  // Get the dimensions of things.
251  const index_type
252  m = A.rows(),
253  n = A.cols(),
254  nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
255 
256  // Validate input
258  n <= 0 || m <= 0 || m > n, std::invalid_argument
259  ,"DirectSparseSolverDense::imp_analyze_and_factor(...) : Error!" );
260 
261  // Extract the matrix in coordinate format
262  Workspace<value_type> a_val(wss,nz);
263  Workspace<index_type> a_row_i(wss,nz);
264  Workspace<index_type> a_col_j(wss,nz);
266  MCTS::EXTRACT_FULL_MATRIX
267  ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
268  ,nz
269  ,&a_val[0]
270  ,nz
271  ,&a_row_i[0]
272  ,&a_col_j[0]
273  );
274 
275  //
276  // Fill the matrix LU = A' so that xGETRF will pivot by row to find
277  // the basis.
278  //
279  // Here we will form the factor of A' = P*L*U where L will be
280  // a n x m upper trapizodial matrix containing the factor lower
281  // triangular factor in LU(1:rank,1:rank) and junk below this.
282  //
283  // Note that xGETRF() pivots by row so if any dependent columns
284  // are found they will always be the last few columns.
285  //
286 
287  // Resize the storage
288  fn.LU_.resize(n,m);
289  fn.ipiv_.resize(n);
290 
291  // Add in the nonzero entires transposed (allows for multiple entries with same
292  // row and column indexes).
293  fn.LU_ = 0.0;
294  for( size_type k = 0; k < nz; ++k )
295  fn.LU_(a_col_j[k],a_row_i[k]) += a_val[k];
296 
297  //
298  // Have xGETRF factor this matrix.
299  //
300 
301  DenseLinAlgLAPack::getrf( &fn.LU_(), &fn.ipiv_[0], &fs.rank_ );
302 
303  // Remember the dimensions
304  fs.m_ = m;
305  fs.n_ = n;
306  fs.nz_ = nz;
307 
308  //
309  // At this point it is important to understand exactly what
310  // ipiv() represents. Each entry in ipiv(k) represents a row
311  // interchange A(k) <=> A(ipiv(k)). Therefore, we have to
312  // do the same row interchanges to the identity permutation
313  // of col_perm to form the permutation array expected by the
314  // DSS interface.
315  //
316 
317  // Form fs.col_perm_
318  fs.col_perm_.resize(n);
320  Workspace<index_type> col_perm_unsorted(wss,fs.rank_);
321  if( m == n && n == fs.rank_ ) {
322  // Leave fs.col_perm_ as identity
323  fn.rect_analyze_and_factor_ = false;
324  }
325  else {
326  fn.rect_analyze_and_factor_ = true;
327  // Permute fs.col_perm_ and save them in col_perm_unsorted
328  for( index_type k = 1; k <= fs.rank_; ++k ) {
329  my_swap( &fs.col_perm_(k), &fs.col_perm_(fn.ipiv_[k-1]) );
330  col_perm_unsorted[k-1] = fs.col_perm_(k);
331  }
332  // Sort the basis selection
333  std::sort(&(fs.col_perm_)[0] , &(fs.col_perm_)[0] + fs.rank_ );
334  std::sort(&(fs.col_perm_)[0] + fs.rank_, &(fs.col_perm_)[0] + n );
335  }
336 
337  // Form the inverse permutation
338  fs.inv_col_perm_.resize(n);
340 
341  if( !(m == n && n == fs.rank_) ) {
342  // Form fn.basis_perm_ and set fs.ipiv_ to identity
343  fn.basis_perm_.resize(fs.rank_);
344  for( size_type k = 1; k <= fs.rank_; ++k ) {
345  fn.basis_perm_(k) = fs.inv_col_perm_(col_perm_unsorted[k-1]);
346  fn.ipiv_[k-1] = k;
347  }
348  }
349 
350  // Copy the data to the output
351  row_perm->resize(m);
352  col_perm->resize(n);
353  *rank = fs.rank_;
355  *col_perm = fs.col_perm_;
356 
357 }
358 
362  ,FactorizationNonzeros *fact_nonzeros
363  ,std::ostream *out
364  )
365 {
366  namespace mmp = MemMngPack;
367  using Teuchos::dyn_cast;
368  typedef MatrixConvertToSparse MCTS;
369  using Teuchos::Workspace;
371 
372  if(out)
373  *out << "\nUsing LAPACK xGETRF to refactor the basis matrix ...\n";
374 
375  // Get the concrete factorization and nonzeros objects
379  &fn = dyn_cast<FactorizationNonzerosDense>(*fact_nonzeros);
380 
381  // Get the dimensions of things.
382  const index_type
383  m = A.rows(),
384  n = A.cols(),
385  nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
386 
387  // Validate input
389  (m != fs.m_ || n != fs.n_ || nz != fs.nz_), std::invalid_argument
390  ,"DirectSparseSolverDense::imp_factor(...): Error!"
391  );
392 
393  // Extract the matrix in coordinate format
394  Workspace<value_type> a_val(wss,nz);
395  Workspace<index_type> a_row_i(wss,nz);
396  Workspace<index_type> a_col_j(wss,nz);
398  MCTS::EXTRACT_FULL_MATRIX
399  ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
400  ,nz
401  ,&a_val[0]
402  ,nz
403  ,&a_row_i[0]
404  ,&a_col_j[0]
405  );
406 
407  //
408  // Fill the matrix LU = B so that xGETRF will pivot by row to find
409  // the basis. Here B is the basis matrix from A'.
410  //
411  // Here we will form the factor of B = P*L*U where L will be
412  // a rank x rank lower triangular.
413  //
414 
415  // Resize the storage
416  fn.rect_analyze_and_factor_ = false;
417  fn.LU_.resize(fs.rank_,fs.rank_);
418  fn.ipiv_.resize(fs.rank_);
419 
420  // Copy only the basis entries (transposed)
421  fn.LU_ = 0.0;
422  for( size_type k = 0; k < nz; ++k ) {
423  const index_type B_i = fs.inv_col_perm_(a_col_j[k]);
424  const index_type B_j = a_row_i[k];
425  if( B_i <= fs.rank_ && B_j <= fs.rank_ )
426  fn.LU_(B_i,B_j) = a_val[k];
427  }
428 
429  //
430  // Have xGETRF factor this matrix.
431  //
432 
433  FortranTypes::f_int B_rank = 0;
434  DenseLinAlgLAPack::getrf( &fn.LU_(), &fn.ipiv_[0], &B_rank );
435 
437  B_rank != fs.rank_, FactorizationFailure
438  ,"DirectSparseSolverDense::imp_factor(...): Error, the basis matrix is no "
439  "longer full rank with B_rank = " << B_rank << " != fs.rank = " << fs.rank_ << "!"
440  );
441 
442 }
443 
444 // private
445 
446 } // end namespace AbstractLinAlgPack
void imp_factor(const AbstractLinAlgPack::MatrixConvertToSparse &A, const FactorizationStructure &fact_struc, FactorizationNonzeros *fact_nonzeros, std::ostream *out)
virtual void coor_extract_nonzeros(EExtractRegion extract_region, EElementUniqueness element_uniqueness, const index_type len_Aval, value_type Aval[], const index_type len_Aij, index_type Arow[], index_type Acol[], const index_type row_offset=0, const index_type col_offset=0) const =0
Extract sparse elements in a coordinate data structure.
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
FortranTypes::f_int f_int
void getrs(const DMatrixSlice &LU, const FortranTypes::f_int ipiv[], BLAS_Cpp::Transp transp, DMatrixSlice *B)
Calls xGETRS to solve linear systems with the factors of P'*A = L*U generated by xGETRF.
void inv_perm(const IVector &perm, IVector *inv_perm)
const BasisMatrix::fact_struc_ptr_t & get_fact_struc() const
virtual index_type num_nonzeros(EExtractRegion extract_region, EElementUniqueness element_uniqueness) const =0
Returns the number of nonzeros in the matrix.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void getrf(DMatrixSlice *A, FortranTypes::f_int ipiv[], FortranTypes::f_int *rank)
Calls xGETRF to compute the P'*A = L*U factorization of a general rectuangular matrix.
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)
Transposed.
Not transposed.
virtual size_type cols() const
Return the number of columns in the matrix.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T_To & dyn_cast(T_From &from)
Mix-in interface for extracing explicit elements from a sparse matrix in one of several Fortran compa...
const Teuchos::RCP< FactorizationStructure > create_fact_struc() const
Extract a constant DenseLinAlgPack::DVectorSlice view of a Vector object.
const Teuchos::RCP< FactorizationNonzeros > create_fact_nonzeros() const
std::ostream * out
const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_dbl_prec const LAPACK_C_Decl::f_int const LAPACK_C_Decl::f_int LAPACK_C_Decl::f_dbl_prec B[]
void resize(size_type rows, size_type cols, value_type val=value_type())
Resize matrix to a (rows x cols) matrix and initializes any added elements by val.
Abstract class for objects that represent the factorization structure of a particular matrix...
void identity_perm(IVector *perm)
Transp trans_not(Transp _trans)
Return the opposite of the transpose argument.
Abstract interface for mutable coordinate vectors {abstract}.
Extract a non-const DenseLinAlgPack::DVectorSlice view of a VectorMutable object. ...
Abstract class for objects that represent the factorization nonzeros of a particular matrix...
virtual size_type rows() const
Return the number of rows in the matrix.
void inv_perm_ele(const DVectorSlice &y, const IVector &perm, DVectorSlice *x)
Perform x = P'*y.
Transp
TRANS.
int n
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
void perm_ele(const IVector &perm, DVectorSlice *vs)
Permute a DVectorSlice in place.