DUNE PDELab (git)

umfpack.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_UMFPACK_HH
6#define DUNE_ISTL_UMFPACK_HH
7
8#if HAVE_SUITESPARSE_UMFPACK || defined DOXYGEN
9
10#include<complex>
11#include<type_traits>
12
13#include<umfpack.h>
14
18#include<dune/istl/bccsmatrixinitializer.hh>
20#include<dune/istl/foreach.hh>
21#include<dune/istl/multitypeblockmatrix.hh>
22#include<dune/istl/multitypeblockvector.hh>
25#include <dune/istl/solverfactory.hh>
26
27
28
29namespace Dune {
41 // FORWARD DECLARATIONS
42 template<class M, class T, class TM, class TD, class TA>
43 class SeqOverlappingSchwarz;
44
45 template<class T, bool tag>
46 struct SeqOverlappingSchwarzAssemblerHelper;
47
48 // wrapper class for C-Function Calls in the backend. Choose the right function namespace
49 // depending on the template parameter used.
50 template<typename T>
51 struct UMFPackMethodChooser
52 {
53 static constexpr bool valid = false ;
54 };
55
56 template<>
57 struct UMFPackMethodChooser<double>
58 {
59 static constexpr bool valid = true ;
60
61 template<typename... A>
62 static void defaults(A... args)
63 {
64 umfpack_dl_defaults(args...);
65 }
66 template<typename... A>
67 static void free_numeric(A... args)
68 {
69 umfpack_dl_free_numeric(args...);
70 }
71 template<typename... A>
72 static void free_symbolic(A... args)
73 {
74 umfpack_dl_free_symbolic(args...);
75 }
76 template<typename... A>
77 static int load_numeric(A... args)
78 {
79 return umfpack_dl_load_numeric(args...);
80 }
81 template<typename... A>
82 static void numeric(A... args)
83 {
84 umfpack_dl_numeric(args...);
85 }
86 template<typename... A>
87 static void report_info(A... args)
88 {
89 umfpack_dl_report_info(args...);
90 }
91 template<typename... A>
92 static void report_status(A... args)
93 {
94 umfpack_dl_report_status(args...);
95 }
96 template<typename... A>
97 static int save_numeric(A... args)
98 {
99 return umfpack_dl_save_numeric(args...);
100 }
101 template<typename... A>
102 static void solve(A... args)
103 {
104 umfpack_dl_solve(args...);
105 }
106 template<typename... A>
107 static void symbolic(A... args)
108 {
109 umfpack_dl_symbolic(args...);
110 }
111 };
112
113 template<>
114 struct UMFPackMethodChooser<std::complex<double> >
115 {
116 static constexpr bool valid = true ;
117 using size_type = SuiteSparse_long;
118
119 template<typename... A>
120 static void defaults(A... args)
121 {
122 umfpack_zl_defaults(args...);
123 }
124 template<typename... A>
125 static void free_numeric(A... args)
126 {
127 umfpack_zl_free_numeric(args...);
128 }
129 template<typename... A>
130 static void free_symbolic(A... args)
131 {
132 umfpack_zl_free_symbolic(args...);
133 }
134 template<typename... A>
135 static int load_numeric(A... args)
136 {
137 return umfpack_zl_load_numeric(args...);
138 }
139 template<typename... A>
140 static void numeric(const size_type* cs, const size_type* ri, const double* val, A... args)
141 {
142 umfpack_zl_numeric(cs,ri,val,NULL,args...);
143 }
144 template<typename... A>
145 static void report_info(A... args)
146 {
147 umfpack_zl_report_info(args...);
148 }
149 template<typename... A>
150 static void report_status(A... args)
151 {
152 umfpack_zl_report_status(args...);
153 }
154 template<typename... A>
155 static int save_numeric(A... args)
156 {
157 return umfpack_zl_save_numeric(args...);
158 }
159 template<typename... A>
160 static void solve(size_type m, const size_type* cs, const size_type* ri, std::complex<double>* val, double* x, const double* b,A... args)
161 {
162 const double* cval = reinterpret_cast<const double*>(val);
163 umfpack_zl_solve(m,cs,ri,cval,NULL,x,NULL,b,NULL,args...);
164 }
165 template<typename... A>
166 static void symbolic(size_type m, size_type n, const size_type* cs, const size_type* ri, const double* val, A... args)
167 {
168 umfpack_zl_symbolic(m,n,cs,ri,val,NULL,args...);
169 }
170 };
171
172 namespace Impl
173 {
174 template<class M, class = void>
175 struct UMFPackVectorChooser;
176
178 template<class M> using UMFPackDomainType = typename UMFPackVectorChooser<M>::domain_type;
179
181 template<class M> using UMFPackRangeType = typename UMFPackVectorChooser<M>::range_type;
182
183 template<class M>
184 struct UMFPackVectorChooser<M,
185 std::enable_if_t<(std::is_same<M,double>::value) || (std::is_same<M,std::complex<double> >::value)>>
186 {
187 using domain_type = M;
188 using range_type = M;
189 };
190
191 template<typename T, int n, int m>
192 struct UMFPackVectorChooser<FieldMatrix<T,n,m>,
193 std::enable_if_t<(std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value)>>
194 {
196 using domain_type = FieldVector<T,m>;
198 using range_type = FieldVector<T,n>;
199 };
200
201 template<typename T, typename A>
202 struct UMFPackVectorChooser<BCRSMatrix<T,A>,
203 std::void_t<UMFPackDomainType<T>, UMFPackRangeType<T>>>
204 {
205 // In case of recursive deduction (e.g., BCRSMatrix<FieldMatrix<...>, Allocator<FieldMatrix<...>>>)
206 // the allocator needs to be converted to the sub-block allocator type too (e.g., Allocator<FieldVector<...>>).
207 // Note that matrix allocator is assumed to be the same as the domain/range type of allocators
209 using domain_type = BlockVector<UMFPackDomainType<T>, typename std::allocator_traits<A>::template rebind_alloc<UMFPackDomainType<T>>>;
211 using range_type = BlockVector<UMFPackRangeType<T>, typename std::allocator_traits<A>::template rebind_alloc<UMFPackRangeType<T>>>;
212 };
213
214 // to make the `UMFPackVectorChooser` work with `MultiTypeBlockMatrix`, we need to add an intermediate step for the rows, which are typically `MultiTypeBlockVector`
215 template<typename FirstBlock, typename... Blocks>
216 struct UMFPackVectorChooser<MultiTypeBlockVector<FirstBlock, Blocks...>,
217 std::void_t<UMFPackDomainType<FirstBlock>, UMFPackRangeType<FirstBlock>, UMFPackDomainType<Blocks>...>>
218 {
220 using domain_type = MultiTypeBlockVector<UMFPackDomainType<FirstBlock>, UMFPackDomainType<Blocks>...>;
222 using range_type = UMFPackRangeType<FirstBlock>;
223 };
224
225 // specialization for `MultiTypeBlockMatrix` with `MultiTypeBlockVector` rows
226 template<typename FirstRow, typename... Rows>
227 struct UMFPackVectorChooser<MultiTypeBlockMatrix<FirstRow, Rows...>,
228 std::void_t<UMFPackDomainType<FirstRow>, UMFPackRangeType<FirstRow>, UMFPackRangeType<Rows>...>>
229 {
231 using domain_type = UMFPackDomainType<FirstRow>;
233 using range_type = MultiTypeBlockVector< UMFPackRangeType<FirstRow>, UMFPackRangeType<Rows>... >;
234 };
235
236 // dummy class to represent no BitVector
237 struct NoBitVector
238 {};
239
240
241 }
242
256 template<typename M>
257 class UMFPack : public InverseOperator<Impl::UMFPackDomainType<M>,Impl::UMFPackRangeType<M>>
258 {
259 using T = typename M::field_type;
260
261 public:
262 using size_type = SuiteSparse_long;
263
265 using Matrix = M;
266 using matrix_type = M;
268 using UMFPackMatrix = ISTL::Impl::BCCSMatrix<typename Matrix::field_type, size_type>;
270 using MatrixInitializer = ISTL::Impl::BCCSMatrixInitializer<M, size_type>;
272 using domain_type = Impl::UMFPackDomainType<M>;
274 using range_type = Impl::UMFPackRangeType<M>;
275
278 {
279 return SolverCategory::Category::sequential;
280 }
281
290 UMFPack(const Matrix& matrix, int verbose=0) : matrixIsLoaded_(false)
291 {
292 //check whether T is a supported type
293 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
294 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
295 Caller::defaults(UMF_Control);
296 setVerbosity(verbose);
297 setMatrix(matrix);
298 }
299
308 UMFPack(const Matrix& matrix, int verbose, bool) : matrixIsLoaded_(false)
309 {
310 //check whether T is a supported type
311 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
312 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
313 Caller::defaults(UMF_Control);
314 setVerbosity(verbose);
315 setMatrix(matrix);
316 }
317
327 UMFPack(const Matrix& mat_, const ParameterTree& config)
328 : UMFPack(mat_, config.get<int>("verbose", 0))
329 {}
330
333 UMFPack() : matrixIsLoaded_(false), verbosity_(0)
334 {
335 //check whether T is a supported type
336 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
337 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
338 Caller::defaults(UMF_Control);
339 }
340
351 UMFPack(const Matrix& mat_, const char* file, int verbose=0)
352 {
353 //check whether T is a supported type
354 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
355 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
356 Caller::defaults(UMF_Control);
357 setVerbosity(verbose);
358 int errcode = Caller::load_numeric(&UMF_Numeric, const_cast<char*>(file));
359 if ((errcode == UMFPACK_ERROR_out_of_memory) || (errcode == UMFPACK_ERROR_file_IO))
360 {
361 matrixIsLoaded_ = false;
362 setMatrix(mat_);
363 saveDecomposition(file);
364 }
365 else
366 {
367 matrixIsLoaded_ = true;
368 std::cout << "UMFPack decomposition successfully loaded from " << file << std::endl;
369 }
370 }
371
378 UMFPack(const char* file, int verbose=0)
379 {
380 //check whether T is a supported type
381 static_assert((std::is_same<T,double>::value) || (std::is_same<T,std::complex<double> >::value),
382 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
383 Caller::defaults(UMF_Control);
384 int errcode = Caller::load_numeric(&UMF_Numeric, const_cast<char*>(file));
385 if (errcode == UMFPACK_ERROR_out_of_memory)
386 DUNE_THROW(Dune::Exception, "ran out of memory while loading UMFPack decomposition");
387 if (errcode == UMFPACK_ERROR_file_IO)
388 DUNE_THROW(Dune::Exception, "IO error while loading UMFPack decomposition");
389 matrixIsLoaded_ = true;
390 std::cout << "UMFPack decomposition successfully loaded from " << file << std::endl;
391 setVerbosity(verbose);
392 }
393
394 virtual ~UMFPack()
395 {
396 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
397 free();
398 }
399
404 {
405 if (umfpackMatrix_.N() != b.dim())
406 DUNE_THROW(Dune::ISTLError, "Size of right-hand-side vector b does not match the number of matrix rows!");
407 if (umfpackMatrix_.M() != x.dim())
408 DUNE_THROW(Dune::ISTLError, "Size of solution vector x does not match the number of matrix columns!");
409 if (b.size() == 0)
410 return;
411
412 // we have to convert x and b into flat structures
413 // however, this is linear in time
414 std::vector<T> xFlat(x.dim()), bFlat(b.dim());
415
416 flatVectorForEach(x, [&](auto&& entry, auto&& offset)
417 {
418 xFlat[ offset ] = entry;
419 });
420
421 flatVectorForEach(b, [&](auto&& entry, auto&& offset)
422 {
423 bFlat[ offset ] = entry;
424 });
425
426 double UMF_Apply_Info[UMFPACK_INFO];
427 Caller::solve(UMFPACK_A,
428 umfpackMatrix_.getColStart(),
429 umfpackMatrix_.getRowIndex(),
430 umfpackMatrix_.getValues(),
431 reinterpret_cast<double*>(&xFlat[0]),
432 reinterpret_cast<double*>(&bFlat[0]),
433 UMF_Numeric,
434 UMF_Control,
435 UMF_Apply_Info);
436
437 // copy back to blocked vector
438 flatVectorForEach(x, [&](auto&& entry, auto offset)
439 {
440 entry = xFlat[offset];
441 });
442
443 //this is a direct solver
444 res.iterations = 1;
445 res.converged = true;
446 res.elapsed = UMF_Apply_Info[UMFPACK_SOLVE_WALLTIME];
447
448 printOnApply(UMF_Apply_Info);
449 }
450
454 void apply (domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
455 {
456 apply(x,b,res);
457 }
458
466 void apply(T* x, T* b)
467 {
468 double UMF_Apply_Info[UMFPACK_INFO];
469 Caller::solve(UMFPACK_A,
470 umfpackMatrix_.getColStart(),
471 umfpackMatrix_.getRowIndex(),
472 umfpackMatrix_.getValues(),
473 x,
474 b,
475 UMF_Numeric,
476 UMF_Control,
477 UMF_Apply_Info);
478 printOnApply(UMF_Apply_Info);
479 }
480
492 void setOption(unsigned int option, double value)
493 {
494 if (option >= UMFPACK_CONTROL)
495 DUNE_THROW(RangeError, "Requested non-existing UMFPack option");
496
497 UMF_Control[option] = value;
498 }
499
503 void saveDecomposition(const char* file)
504 {
505 int errcode = Caller::save_numeric(UMF_Numeric, const_cast<char*>(file));
506 if (errcode != UMFPACK_OK)
507 DUNE_THROW(Dune::Exception,"IO ERROR while trying to save UMFPack decomposition");
508 }
509
519 template<class BitVector = Impl::NoBitVector>
520 void setMatrix(const Matrix& matrix, const BitVector& bitVector = {})
521 {
522 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
523 free();
524 if (matrix.N() == 0 or matrix.M() == 0)
525 return;
526
527 if (umfpackMatrix_.N() + umfpackMatrix_.M() + umfpackMatrix_.nonzeroes() != 0)
528 umfpackMatrix_.free();
529
530 constexpr bool useBitVector = not std::is_same_v<BitVector,Impl::NoBitVector>;
531
532 // use a dynamic flat vector for the bitset
533 std::vector<bool> flatBitVector;
534 // and a mapping from the compressed indices
535 std::vector<size_type> subIndices;
536
537 [[maybe_unused]] int numberOfIgnoredDofs = 0;
538 int nonZeros = 0;
539
540 if constexpr ( useBitVector )
541 {
542 auto flatSize = flatVectorForEach(bitVector, [](auto&&, auto&&){});
543 flatBitVector.resize(flatSize);
544
545 flatVectorForEach(bitVector, [&](auto&& entry, auto&& offset)
546 {
547 flatBitVector[ offset ] = entry;
548 if ( entry )
549 {
550 numberOfIgnoredDofs++;
551 }
552 });
553 }
554
555 // compute the flat dimension and the number of nonzeros of the matrix
556 auto [flatRows,flatCols] = flatMatrixForEach( matrix, [&](auto&& /*entry*/, auto&& row, auto&& col){
557 // do not count ignored entries
558 if constexpr ( useBitVector )
559 if ( flatBitVector[row] or flatBitVector[col] )
560 return;
561
562 nonZeros++;
563 });
564
565 if constexpr ( useBitVector )
566 {
567 // use the original flatRows!
568 subIndices.resize(flatRows,std::numeric_limits<std::size_t>::max());
569
570 size_type subIndexCounter = 0;
571 for ( size_type i=0; i<size_type(flatRows); i++ )
572 if ( not flatBitVector[ i ] )
573 subIndices[ i ] = subIndexCounter++;
574
575 // update the original matrix size
576 flatRows -= numberOfIgnoredDofs;
577 flatCols -= numberOfIgnoredDofs;
578 }
579
580
581 umfpackMatrix_.setSize(flatRows,flatCols);
582 umfpackMatrix_.Nnz_ = nonZeros;
583
584 // prepare the arrays
585 umfpackMatrix_.colstart = new size_type[flatCols+1];
586 umfpackMatrix_.rowindex = new size_type[nonZeros];
587 umfpackMatrix_.values = new T[nonZeros];
588
589 for ( size_type i=0; i<size_type(flatCols+1); i++ )
590 {
591 umfpackMatrix_.colstart[i] = 0;
592 }
593
594 // at first, we need to compute the column start indices
595 // therefore, we count all entries in each column and in the end we accumulate everything
596 flatMatrixForEach(matrix, [&](auto&& /*entry*/, auto&& flatRowIndex, auto&& flatColIndex)
597 {
598 // do nothing if entry is excluded
599 if constexpr ( useBitVector )
600 if ( flatBitVector[flatRowIndex] or flatBitVector[flatColIndex] )
601 return;
602
603 // pick compressed or uncompressed index
604 // compiler will hopefully do some constexpr optimization here
605 auto colIdx = useBitVector ? subIndices[flatColIndex] : flatColIndex;
606
607 umfpackMatrix_.colstart[colIdx+1]++;
608 });
609
610 // now accumulate
611 for ( size_type i=0; i<(size_type)flatCols; i++ )
612 {
613 umfpackMatrix_.colstart[i+1] += umfpackMatrix_.colstart[i];
614 }
615
616 // we need a compressed position counter in each column
617 std::vector<size_type> colPosition(flatCols,0);
618
619 // now we can set the entries: the procedure below works with both row- or column major index ordering
620 flatMatrixForEach(matrix, [&](auto&& entry, auto&& flatRowIndex, auto&& flatColIndex)
621 {
622 // do nothing if entry is excluded
623 if constexpr ( useBitVector )
624 if ( flatBitVector[flatRowIndex] or flatBitVector[flatColIndex] )
625 return;
626
627 // pick compressed or uncompressed index
628 // compiler will hopefully do some constexpr optimization here
629 auto rowIdx = useBitVector ? subIndices[flatRowIndex] : flatRowIndex;
630 auto colIdx = useBitVector ? subIndices[flatColIndex] : flatColIndex;
631
632 // the start index of each column is already fixed
633 auto colStart = umfpackMatrix_.colstart[colIdx];
634 // get the current number of picked elements in this column
635 auto colPos = colPosition[colIdx];
636 // assign the corresponding row index and the value of this element
637 umfpackMatrix_.rowindex[ colStart + colPos ] = rowIdx;
638 umfpackMatrix_.values[ colStart + colPos ] = entry;
639 // increase the number of picked elements in this column
640 colPosition[colIdx]++;
641 });
642
643 decompose();
644 }
645
646 // Keep legacy version using a set of scalar indices
647 // The new version using a bitVector type for marking the active matrix indices is
648 // directly given in `setMatrix` with an additional BitVector argument.
649 // The new version is more flexible and allows, e.g., marking single components of a matrix block.
650 template<typename S>
651 void setSubMatrix(const Matrix& _mat, const S& rowIndexSet)
652 {
653 if ((umfpackMatrix_.N() + umfpackMatrix_.M() > 0) || matrixIsLoaded_)
654 free();
655
656 if (umfpackMatrix_.N() + umfpackMatrix_.M() + umfpackMatrix_.nonzeroes() != 0)
657 umfpackMatrix_.free();
658
659 umfpackMatrix_.setSize(rowIndexSet.size()*MatrixDimension<Matrix>::rowdim(_mat) / _mat.N(),
660 rowIndexSet.size()*MatrixDimension<Matrix>::coldim(_mat) / _mat.M());
661 ISTL::Impl::BCCSMatrixInitializer<Matrix, SuiteSparse_long> initializer(umfpackMatrix_);
662
663 copyToBCCSMatrix(initializer, ISTL::Impl::MatrixRowSubset<Matrix,std::set<std::size_t> >(_mat,rowIndexSet));
664
665 decompose();
666 }
667
675 void setVerbosity(int v)
676 {
677 verbosity_ = v;
678 // set the verbosity level in UMFPack
679 if (verbosity_ == 0)
680 UMF_Control[UMFPACK_PRL] = 1;
681 if (verbosity_ == 1)
682 UMF_Control[UMFPACK_PRL] = 2;
683 if (verbosity_ == 2)
684 UMF_Control[UMFPACK_PRL] = 4;
685 }
686
692 {
693 return UMF_Numeric;
694 }
695
701 {
702 return umfpackMatrix_;
703 }
704
709 void free()
710 {
711 if (!matrixIsLoaded_)
712 {
713 Caller::free_symbolic(&UMF_Symbolic);
714 umfpackMatrix_.free();
715 }
716 Caller::free_numeric(&UMF_Numeric);
717 matrixIsLoaded_ = false;
718 }
719
720 const char* name() { return "UMFPACK"; }
721
722 private:
723 typedef typename Dune::UMFPackMethodChooser<T> Caller;
724
725 template<class Mat,class X, class TM, class TD, class T1>
726 friend class SeqOverlappingSchwarz;
727 friend struct SeqOverlappingSchwarzAssemblerHelper<UMFPack<Matrix>,true>;
728
730 void decompose()
731 {
732 double UMF_Decomposition_Info[UMFPACK_INFO];
733 Caller::symbolic(static_cast<SuiteSparse_long>(umfpackMatrix_.N()),
734 static_cast<SuiteSparse_long>(umfpackMatrix_.N()),
735 umfpackMatrix_.getColStart(),
736 umfpackMatrix_.getRowIndex(),
737 reinterpret_cast<double*>(umfpackMatrix_.getValues()),
738 &UMF_Symbolic,
739 UMF_Control,
740 UMF_Decomposition_Info);
741 Caller::numeric(umfpackMatrix_.getColStart(),
742 umfpackMatrix_.getRowIndex(),
743 reinterpret_cast<double*>(umfpackMatrix_.getValues()),
744 UMF_Symbolic,
745 &UMF_Numeric,
746 UMF_Control,
747 UMF_Decomposition_Info);
748 Caller::report_status(UMF_Control,UMF_Decomposition_Info[UMFPACK_STATUS]);
749 if (verbosity_ == 1)
750 {
751 std::cout << "[UMFPack Decomposition]" << std::endl;
752 std::cout << "Wallclock Time taken: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_WALLTIME] << " (CPU Time: " << UMF_Decomposition_Info[UMFPACK_NUMERIC_TIME] << ")" << std::endl;
753 std::cout << "Flops taken: " << UMF_Decomposition_Info[UMFPACK_FLOPS] << std::endl;
754 std::cout << "Peak Memory Usage: " << UMF_Decomposition_Info[UMFPACK_PEAK_MEMORY]*UMF_Decomposition_Info[UMFPACK_SIZE_OF_UNIT] << " bytes" << std::endl;
755 std::cout << "Condition number estimate: " << 1./UMF_Decomposition_Info[UMFPACK_RCOND] << std::endl;
756 std::cout << "Numbers of non-zeroes in decomposition: L: " << UMF_Decomposition_Info[UMFPACK_LNZ] << " U: " << UMF_Decomposition_Info[UMFPACK_UNZ] << std::endl;
757 }
758 if (verbosity_ == 2)
759 {
760 Caller::report_info(UMF_Control,UMF_Decomposition_Info);
761 }
762 }
763
764 void printOnApply(double* UMF_Info)
765 {
766 Caller::report_status(UMF_Control,UMF_Info[UMFPACK_STATUS]);
767 if (verbosity_ > 0)
768 {
769 std::cout << "[UMFPack Solve]" << std::endl;
770 std::cout << "Wallclock Time: " << UMF_Info[UMFPACK_SOLVE_WALLTIME] << " (CPU Time: " << UMF_Info[UMFPACK_SOLVE_TIME] << ")" << std::endl;
771 std::cout << "Flops Taken: " << UMF_Info[UMFPACK_SOLVE_FLOPS] << std::endl;
772 std::cout << "Iterative Refinement steps taken: " << UMF_Info[UMFPACK_IR_TAKEN] << std::endl;
773 std::cout << "Error Estimate: " << UMF_Info[UMFPACK_OMEGA1] << " resp. " << UMF_Info[UMFPACK_OMEGA2] << std::endl;
774 }
775 }
776
777 UMFPackMatrix umfpackMatrix_;
778 bool matrixIsLoaded_;
779 int verbosity_;
780 void *UMF_Symbolic;
781 void *UMF_Numeric;
782 double UMF_Control[UMFPACK_CONTROL];
783 };
784
785 template<typename T, typename A, int n, int m>
786 struct IsDirectSolver<UMFPack<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
787 {
788 enum { value=true};
789 };
790
791 template<typename T, typename A>
792 struct StoresColumnCompressed<UMFPack<BCRSMatrix<T,A> > >
793 {
794 enum { value = true };
795 };
796
797 namespace UMFPackImpl {
798 template<class OpTraits, class=void> struct isValidBlock : std::false_type{};
799 template<class OpTraits> struct isValidBlock<OpTraits,
800 std::enable_if_t<
801 std::is_same_v<Impl::UMFPackDomainType<typename OpTraits::matrix_type>, typename OpTraits::domain_type>
802 && std::is_same_v<Impl::UMFPackRangeType<typename OpTraits::matrix_type>, typename OpTraits::range_type>
803 && std::is_same_v<typename FieldTraits<typename OpTraits::domain_type::field_type>::real_type, double>
804 && std::is_same_v<typename FieldTraits<typename OpTraits::range_type::field_type>::real_type, double>
805 >> : std::true_type {};
806 }
807 DUNE_REGISTER_SOLVER("umfpack",
808 [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
809 -> std::shared_ptr<typename decltype(opTraits)::solver_type>
810 {
811 using OpTraits = decltype(opTraits);
812 // works only for sequential operators
813 if constexpr (OpTraits::isParallel){
814 if(opTraits.getCommOrThrow(op).communicator().size() > 1)
815 DUNE_THROW(Dune::InvalidStateException, "UMFPack works only for sequential operators.");
816 }
817 if constexpr (OpTraits::isAssembled){
818 using M = typename OpTraits::matrix_type;
819 // check if UMFPack<M>* is convertible to
820 // InverseOperator*. This checks compatibility of the
821 // domain and range types
822 if constexpr (UMFPackImpl::isValidBlock<OpTraits>::value) {
823 const auto& A = opTraits.getAssembledOpOrThrow(op);
824 const M& mat = A->getmat();
825 int verbose = config.get("verbose", 0);
826 return std::make_shared<Dune::UMFPack<M>>(mat,verbose);
827 }
828 }
829 DUNE_THROW(UnsupportedType,
830 "Unsupported Type in UMFPack (only double and std::complex<double> supported)");
831 return nullptr;
832 });
833} // end namespace Dune
834
835#endif // HAVE_SUITESPARSE_UMFPACK
836
837#endif //DUNE_ISTL_UMFPACK_HH
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
derive error class from the base class in common
Definition: istlexception.hh:19
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:281
Abstract base class for all solvers.
Definition: solver.hh:101
Hierarchical structure of string parameters.
Definition: parametertree.hh:37
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:188
Default exception class for range errors.
Definition: exceptions.hh:254
The UMFPack direct sparse solver.
Definition: umfpack.hh:258
A few common exception classes.
Implementation of the BCRSMatrix class.
Implementations of the inverse operator interface.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
typename Impl::voider< Types... >::type void_t
Is void for all valid input types. The workhorse for C++11 SFINAE-techniques.
Definition: typetraits.hh:40
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
void free()
free allocated space.
Definition: umfpack.hh:709
UMFPack(const Matrix &mat_, const ParameterTree &config)
Construct a solver object from a matrix.
Definition: umfpack.hh:327
UMFPack(const Matrix &mat_, const char *file, int verbose=0)
Try loading a decomposition from file and do a decomposition if unsuccessful.
Definition: umfpack.hh:351
Impl::UMFPackRangeType< M > range_type
The type of the range of the solver.
Definition: umfpack.hh:274
UMFPack()
default constructor
Definition: umfpack.hh:333
Impl::UMFPackDomainType< M > domain_type
The type of the domain of the solver.
Definition: umfpack.hh:272
void apply(T *x, T *b)
additional apply method with c-arrays in analogy to superlu
Definition: umfpack.hh:466
SolverCategory::Category category() const override
Category of the solver (see SolverCategory::Category)
Definition: umfpack.hh:277
void setVerbosity(int v)
sets the verbosity level for the UMFPack solver
Definition: umfpack.hh:675
UMFPack(const char *file, int verbose=0)
try loading a decomposition from file
Definition: umfpack.hh:378
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res) override
apply inverse operator, with given convergence criteria.
Definition: umfpack.hh:454
void setMatrix(const Matrix &matrix, const BitVector &bitVector={})
Initialize data from given matrix.
Definition: umfpack.hh:520
void saveDecomposition(const char *file)
saves a decomposition to a file
Definition: umfpack.hh:503
UMFPackMatrix & getInternalMatrix()
Return the column compress matrix from UMFPack.
Definition: umfpack.hh:700
UMFPack(const Matrix &matrix, int verbose=0)
Construct a solver object from a matrix.
Definition: umfpack.hh:290
ISTL::Impl::BCCSMatrixInitializer< M, size_type > MatrixInitializer
Type of an associated initializer class.
Definition: umfpack.hh:270
ISTL::Impl::BCCSMatrix< typename Matrix::field_type, size_type > UMFPackMatrix
The corresponding (scalar) UMFPack matrix type.
Definition: umfpack.hh:268
void setOption(unsigned int option, double value)
Set UMFPack-specific options.
Definition: umfpack.hh:492
M Matrix
The matrix type.
Definition: umfpack.hh:265
void apply(domain_type &x, range_type &b, InverseOperatorResult &res) override
Apply inverse operator,.
Definition: umfpack.hh:403
UMFPack(const Matrix &matrix, int verbose, bool)
Constructor for compatibility with SuperLU standard constructor.
Definition: umfpack.hh:308
void * getFactorization()
Return the matrix factorization.
Definition: umfpack.hh:691
Dune namespace.
Definition: alignedallocator.hh:13
std::pair< std::size_t, std::size_t > flatMatrixForEach(Matrix &&matrix, F &&f, std::size_t rowOffset=0, std::size_t colOffset=0)
Traverse a blocked matrix and call a functor at each scalar entry.
Definition: foreach.hh:132
std::size_t flatVectorForEach(Vector &&vector, F &&f, std::size_t offset=0)
Traverse a blocked vector and call a functor at each scalar entry.
Definition: foreach.hh:95
constexpr auto get(std::integer_sequence< T, II... >, std::integral_constant< std::size_t, pos >={})
Return the entry at position pos of the given sequence.
Definition: integersequence.hh:22
STL namespace.
Templates characterizing the type of a solver.
Statistics about the application of an inverse operator.
Definition: solver.hh:50
double elapsed
Elapsed time in seconds.
Definition: solver.hh:84
int iterations
Number of iterations.
Definition: solver.hh:69
bool converged
True if convergence criterion has been met.
Definition: solver.hh:75
Category
Definition: solvercategory.hh:23
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)