61 void my_swap( T* v1, T* v2 )
71 std::valarray<T>& cva(
const std::valarray<T>& va )
73 return const_cast<std::valarray<T>&
>(va);
78 namespace AbstractLinAlgPack {
114 yd().dim() != xd().dim(), std::invalid_argument
115 ,
"DirectSparseSolverDense::BasisMatrixDense::V_InvMtV(...) : Error, "
116 " y.dim() = " << yd().dim() <<
" != x.dim() = " << xd().dim() <<
"!"
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);
200 namespace mmp = MemMngPack;
235 namespace mmp = MemMngPack;
242 *out <<
"\nUsing LAPACK xGETRF to analyze and factor a new matrix ...\n";
254 nz = A.
num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
258 n <= 0 || m <= 0 || m >
n, std::invalid_argument
259 ,
"DirectSparseSolverDense::imp_analyze_and_factor(...) : Error!" );
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
295 fn.
LU_(a_col_j[k],a_row_i[k]) += a_val[k];
320 Workspace<index_type> col_perm_unsorted(wss,fs.
rank_);
321 if( m == n && n == fs.
rank_ ) {
330 col_perm_unsorted[k-1] = fs.
col_perm_(k);
341 if( !(m == n && n == fs.
rank_) ) {
366 namespace mmp = MemMngPack;
373 *out <<
"\nUsing LAPACK xGETRF to refactor the basis matrix ...\n";
385 nz = A.
num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
389 (m != fs.
m_ ||
n != fs.
n_ ||
nz != fs.
nz_), std::invalid_argument
390 ,
"DirectSparseSolverDense::imp_factor(...): Error!"
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
416 fn.rect_analyze_and_factor_ =
false;
418 fn.ipiv_.resize(fs.
rank_);
426 fn.LU_(B_i,B_j) = a_val[k];
438 ,
"DirectSparseSolverDense::imp_factor(...): Error, the basis matrix is no "
439 "longer full rank with B_rank = " << B_rank <<
" != fs.rank = " << fs.
rank_ <<
"!"
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.
FactorizationStructureDense()
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
RTOp_index_type index_type
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.
RTOp_value_type value_type
Stores the factorization structure for Dense.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void estimated_fillin_ratio(value_type estimated_fillin_ratio)
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.
Stores the factorization nonzeros for Dense.
std::valarray< f_int > ipiv_
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)
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)
FortranTypes::f_int rank_
RTOp_index_type size_type
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
const basis_matrix_factory_ptr_t basis_matrix_factory() const
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.
bool rect_analyze_and_factor_
FactorizationStructure * fact_struc
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.
DirectSparseSolverDense()
Default constructor.
Abstract interface for mutable coordinate vectors {abstract}.
Implements the BasisMatrix object for Dense.
Teuchos::RCP< BasisMatrixImp > create_matrix() const
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.
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
void perm_ele(const IVector &perm, DVectorSlice *vs)
Permute a DVectorSlice in place.