27 using namespace boost::numeric::ublas;
31 : L_(A.size1(),A.size2()),
32 U_(A.size1(),A.size2()) {
35 "sparse ILU preconditioner works only with square matrices");
37 for (SparseMatrix::size_type i=0; i <
L_.size1(); ++i)
41 std::set<Integer> lBandSet, uBandSet;
43 compressed_matrix<Integer> levs(
n,
n);
52 std::vector<Integer> levii(
n, 0);
70 std::vector<Integer> nonZeros;
71 std::vector<Real> nonZeroEntries;
72 nonZeros.reserve(uBandSet.size()+1);
73 nonZeroEntries.reserve(uBandSet.size()+1);
74 const Real entry =
U_(jj,jj);
76 nonZeros.push_back(jj);
77 nonZeroEntries.push_back(entry);
79 auto iter = uBandSet.begin();
80 auto end = uBandSet.end();
81 for (; iter != end; ++iter) {
82 const Real entry =
U_(jj,jj+*iter);
84 nonZeros.push_back(jj+*iter);
85 nonZeroEntries.push_back(entry);
89 if(!nonZeroEntries.empty()) {
90 fact /= nonZeroEntries[0];
92 for (
Size k=0; k<nonZeros.size(); ++k) {
94 const Integer temp = levs(jj,j) + jlev ;
97 w[j] = - fact*nonZeroEntries[k];
102 w[j] -= fact*nonZeroEntries[k];
103 levii[j] = std::min(levii[j],temp);
109 std::vector<Integer> wNonZeros;
110 std::vector<Real> wNonZeroEntries;
111 wNonZeros.reserve(w.
size());
112 wNonZeroEntries.reserve(w.
size());
114 const Real entry = w[i];
116 wNonZeros.push_back(i);
117 wNonZeroEntries.push_back(entry);
120 std::vector<Integer> leviiNonZeroEntries;
121 leviiNonZeroEntries.reserve(levii.size());
122 for (
int entry : levii) {
124 leviiNonZeroEntries.push_back(entry);
127 for (
Size k=0; k<wNonZeros.size(); ++k) {
130 L_(ii,j) = wNonZeroEntries[k];
131 lBandSet.insert(ii-j);
134 U_(ii,j) = wNonZeroEntries[k];
135 levs(ii,j) = leviiNonZeroEntries[k];
137 uBandSet.insert(j-ii);
142 lBands_.resize(lBandSet.size());
143 uBands_.resize(uBandSet.size());
144 std::copy(lBandSet.begin(), lBandSet.end(),
lBands_.begin());
145 std::copy(uBandSet.begin(), uBandSet.end(),
uBands_.begin());
170 y[i]-=
L_(i,k)*
y[k]/
L_(i,i);
1-D array used in linear algebra.
Size size() const
dimension of the array
SparseILUPreconditioner(const SparseMatrix &A, Integer lfil=1)
std::vector< Size > lBands_
Array backwardSolve(const Array &y) const
Array apply(const Array &b) const
std::vector< Size > uBands_
Array forwardSolve(const Array &b) const
const SparseMatrix & U() const
const SparseMatrix & L() const
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
ext::function< Real(Real)> b
QL_INTEGER Integer
integer number
std::size_t Size
size of a container
matrix used in linear algebra.
boost::numeric::ublas::compressed_matrix< Real > SparseMatrix
Preconditioner using the Incomplete LU algorithm and sparse matrices.