5#ifndef DUNE_ISTL_SPQR_HH
6#define DUNE_ISTL_SPQR_HH
8#if HAVE_SUITESPARSE_SPQR || defined DOXYGEN
13#include <SuiteSparseQR.hpp>
15#include <dune/common/exceptions.hh>
35 template<
class M,
class T,
class TM,
class TD,
class TA>
36 class SeqOverlappingSchwarz;
38 template<
class T,
bool tag>
39 struct SeqOverlappingSchwarzAssemblerHelper;
46 template<
class Matrix>
63 template<
typename T,
typename A,
int n,
int m>
65 :
public InverseOperator<BlockVector<FieldVector<T,m>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >,
66 BlockVector<FieldVector<T,n>, typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > > >
75 typedef ISTL::Impl::BCCSMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A>,
int>
MatrixInitializer;
95 SPQR(
const Matrix& matrix,
int verbose=0) : matrixIsLoaded_(false), verbose_(verbose)
98 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
99 "Unsupported Type in SPQR (only double and std::complex<double> supported)");
100 cc_ =
new cholmod_common();
101 cholmod_l_start(cc_);
113 SPQR(
const Matrix& matrix,
int verbose,
bool) : matrixIsLoaded_(false), verbose_(verbose)
116 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
117 "Unsupported Type in SPQR (only double and std::complex<double> supported)");
118 cc_ =
new cholmod_common();
119 cholmod_l_start(cc_);
133 :
SPQR(matrix, config.
get<int>(
"verbose", 0))
137 SPQR() : matrixIsLoaded_(false), verbose_(0)
140 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
141 "Unsupported Type in SPQR (only double and std::complex<double> supported)");
142 cc_ =
new cholmod_common();
143 cholmod_l_start(cc_);
149 if ((spqrMatrix_.N() + spqrMatrix_.M() > 0) || matrixIsLoaded_)
151 cholmod_l_finish(cc_);
157 const std::size_t numRows(spqrMatrix_.N());
159 for(std::size_t k = 0; k != numRows/n; ++k)
160 for (
int l = 0; l < n; ++l)
161 (
static_cast<T*
>(B_->x))[n*k+l] = b[k][l];
163 cholmod_dense* BTemp = B_;
164 B_ = SuiteSparseQR_qmult<T>(0, spqrfactorization_, B_, cc_);
165 cholmod_dense* X = SuiteSparseQR_solve<T>(1, spqrfactorization_, B_, cc_);
166 cholmod_l_free_dense(&BTemp, cc_);
168 const std::size_t numCols(spqrMatrix_.M());
170 for(std::size_t k = 0; k != numCols/m; ++k)
171 for (
int l = 0; l < m; ++l)
172 x[k][l] = (
static_cast<T*
>(X->x))[m*k+l];
174 cholmod_l_free_dense(&X, cc_);
180 std::cout<<std::endl<<
"Solving with SuiteSparseQR"<<std::endl;
181 std::cout<<
"Flops Taken: "<<cc_->SPQR_flopcount<<std::endl;
182 std::cout<<
"Analysis Time: "<<cc_->SPQR_analyze_time<<
" s"<<std::endl;
183 std::cout<<
"Factorize Time: "<<cc_->SPQR_factorize_time<<
" s"<<std::endl;
184 std::cout<<
"Backsolve Time: "<<cc_->SPQR_solve_time<<
" s"<<std::endl;
185 std::cout<<
"Peak Memory Usage: "<<cc_->memory_usage<<
" bytes"<<std::endl;
186 std::cout<<
"Rank Estimate: "<<cc_->SPQR_istat[4]<<std::endl<<std::endl;
196 void setOption([[maybe_unused]]
unsigned int option, [[maybe_unused]]
double value)
202 if ((spqrMatrix_.N() + spqrMatrix_.M() > 0) || matrixIsLoaded_)
205 if (spqrMatrix_.N() + spqrMatrix_.M() + spqrMatrix_.nonzeroes() != 0)
209 ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(spqrMatrix_);
211 copyToBCCSMatrix(initializer, matrix);
219 if ((spqrMatrix_.N() + spqrMatrix_.M() > 0) || matrixIsLoaded_)
222 if (spqrMatrix_.N() + spqrMatrix_.M() + spqrMatrix_.nonzeroes() != 0)
227 ISTL::Impl::BCCSMatrixInitializer<Matrix, int> initializer(spqrMatrix_);
229 copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<
Matrix,std::set<std::size_t> >(matrix,rowIndexSet));
249 return spqrfactorization_;
267 cholmod_l_free_sparse(&A_, cc_);
268 cholmod_l_free_dense(&B_, cc_);
269 SuiteSparseQR_free<T>(&spqrfactorization_, cc_);
271 matrixIsLoaded_ =
false;
281 template<
class M,
class X,
class TM,
class TD,
class T1>
289 const std::size_t nrows(spqrMatrix_.N());
290 const std::size_t ncols(spqrMatrix_.M());
291 const std::size_t nnz(spqrMatrix_.getColStart()[ncols]);
294 A_ = cholmod_l_allocate_sparse(nrows, ncols, nnz, 1, 1, 0, 1, cc_);
297 for(std::size_t k = 0; k != (ncols+1); ++k)
298 (
static_cast<long int *
>(A_->p))[k] = spqrMatrix_.getColStart()[k];
300 for(std::size_t k = 0; k != nnz; ++k)
302 (
static_cast<long int*
>(A_->i))[k] = spqrMatrix_.getRowIndex()[k];
303 (
static_cast<T*
>(A_->x))[k] = spqrMatrix_.getValues()[k];
307 B_ = cholmod_l_allocate_dense(nrows, 1, nrows, A_->xtype, cc_);
309 spqrfactorization_=SuiteSparseQR_factorize<T>(SPQR_ORDERING_DEFAULT,SPQR_DEFAULT_TOL,A_,cc_);
312 SPQRMatrix spqrMatrix_;
313 bool matrixIsLoaded_;
318 SuiteSparseQR_factorization<T>* spqrfactorization_;
321 template<
typename T,
typename A>
327 template<
typename T,
typename A>
336 template<
typename TL,
typename M>
337 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
338 typename Dune::TypeListElement<2, TL>::type>>
341 isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
343 int verbose = config.get(
"verbose", 0);
344 return std::make_shared<Dune::SPQR<M>>(
mat,verbose);
348 template<
typename TL,
typename M>
349 std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
350 typename Dune::TypeListElement<2, TL>::type>>
352 std::enable_if_t<!
isValidBlock<
typename Dune::TypeListElement<1, TL>::type::block_type>::value,
int> = 0)
const
355 "Unsupported Type in SPQR (only double and std::complex<double> supported)");
Templates characterizing the type of a solver.
#define DUNE_REGISTER_DIRECT_SOLVER(name,...)
Definition solverregistry.hh:13
Implementations of the inverse operator interface.
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A > matrix_type
Definition spqr.hh:71
virtual ~SPQR()
Destructor.
Definition spqr.hh:147
SPQR(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition spqr.hh:113
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition spqr.hh:82
SPQRMatrix & getInternalMatrix()
Return the column coppressed matrix.
Definition spqr.hh:256
void setMatrix(const Matrix &matrix)
Initialize data from given matrix.
Definition spqr.hh:200
ISTL::Impl::BCCSMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A >, int > MatrixInitializer
Type of an associated initializer class.
Definition spqr.hh:75
Dune::BCRSMatrix< FieldMatrix< T, n, m >, A > Matrix
The matrix type.
Definition spqr.hh:70
SPQR()
Default constructor.
Definition spqr.hh:137
std::shared_ptr< Dune::InverseOperator< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL, const M &mat, const Dune::ParameterTree &config, std::enable_if_t< isValidBlock< typename Dune::TypeListElement< 1, TL >::type::block_type >::value, int >=0) const
Definition spqr.hh:339
const char * name()
Get method name.
Definition spqr.hh:275
SuiteSparseQR_factorization< T > * getFactorization()
Return the matrix factorization.
Definition spqr.hh:247
void setVerbosity(int v)
Sets the verbosity level for the solver.
Definition spqr.hh:238
ISTL::Impl::BCCSMatrix< T, int > SPQRMatrix
The corresponding SuperLU Matrix type.
Definition spqr.hh:73
virtual void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition spqr.hh:191
void free()
Free allocated space.
Definition spqr.hh:265
SPQR(const Matrix &matrix, const ParameterTree &config)
Constructs the SPQR solver.
Definition spqr.hh:132
Dune::BlockVector< FieldVector< T, n >, typename std::allocator_traits< A >::template rebind_alloc< FieldVector< T, n > > > range_type
The type of the range of the solver.
Definition spqr.hh:79
void setSubMatrix(const Matrix &matrix, const S &rowIndexSet)
Definition spqr.hh:217
Dune::BlockVector< FieldVector< T, m >, typename std::allocator_traits< A >::template rebind_alloc< FieldVector< T, m > > > domain_type
The type of the domain of the solver.
Definition spqr.hh:77
virtual void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition spqr.hh:155
void setOption(unsigned int option, double value)
Definition spqr.hh:196
SPQR(const Matrix &matrix, int verbose=0)
Construct a solver object from a BCRSMatrix.
Definition spqr.hh:95
Matrix & mat
Definition matrixmatrix.hh:347
Definition allocator.hh:11
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get(const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition dependency.hh:293
Definition matrixutils.hh:211
A sparse block matrix with compressed row storage.
Definition bcrsmatrix.hh:466
size_type M() const
number of columns (counted in blocks)
Definition bcrsmatrix.hh:1978
size_type N() const
number of rows (counted in blocks)
Definition bcrsmatrix.hh:1972
A vector of blocks with memory management.
Definition bvector.hh:395
Sequential overlapping Schwarz preconditioner.
Definition overlappingschwarz.hh:755
Definition overlappingschwarz.hh:694
Definition matrixutils.hh:27
Statistics about the application of an inverse operator.
Definition solver.hh:48
int iterations
Number of iterations.
Definition solver.hh:67
bool converged
True if convergence criterion has been met.
Definition solver.hh:73
Abstract base class for all solvers.
Definition solver.hh:99
Category
Definition solvercategory.hh:23
@ sequential
Category for sequential solvers.
Definition solvercategory.hh:25
Definition solverregistry.hh:77
Definition solvertype.hh:16
@ value
Whether this is a direct solver.
Definition solvertype.hh:24
Definition solvertype.hh:30
@ value
whether the solver internally uses column compressed storage
Definition solvertype.hh:36
Use the SPQR package to directly solve linear systems – empty default class.
Definition spqr.hh:48