1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4Copyright (C) 2009 Ralph Schreyer
5
6This file is part of QuantLib, a free-software/open-source library
7for financial quantitative analysts and developers - http://quantlib.org/
8
9QuantLib is free software: you can redistribute it and/or modify it
10under the terms of the QuantLib license. You should have received a
11copy of the license along with this program; if not, please email
12<quantlib-dev@lists.sf.net>. The license is also available online at
13<http://quantlib.org/license.shtml>.
14
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#include <ql/math/matrixutilities/sparseilupreconditioner.hpp>
21#include <ql/math/matrix.hpp>
22
23#include <set>
24
25namespace QuantLib {
26
27 using namespace boost::numeric::ublas;
28
29 SparseILUPreconditioner::SparseILUPreconditioner(const SparseMatrix& A,
30 Integer lfil)
31 : L_(A.size1(),A.size2()),
32 U_(A.size1(),A.size2()) {
33
34 QL_REQUIRE(A.size1() == A.size2(),
35 "sparse ILU preconditioner works only with square matrices");
36
37 for (SparseMatrix::size_type i=0; i < L_.size1(); ++i)
38 L_(i,i) = 1.0;
39
40 const Integer n = A.size1();
41 std::set<Integer> lBandSet, uBandSet;
42
43 compressed_matrix<Integer> levs(n,n);
44 Integer lfilp = lfil + 1;
45
46 for (Integer ii=0; ii<n; ++ii) {
47 Array w(n);
48 for(Integer k=0; k<n; ++k) {
49 w[k] = A(ii,k);
50 }
51
52 std::vector<Integer> levii(n, 0);
53 for (Integer i=0; i<n; ++i) {
54 if( w[i] > QL_EPSILON
55 || w[i] < -1.0*QL_EPSILON) levii[i] = 1;
56 }
57 Integer jj = -1;
58 while (jj < ii) {
59 for (Integer k=jj+1; k<n; ++k) {
60 if (levii[k] != 0) {
61 jj = k;
62 break;
63 }
64 }
65 if (jj >= ii) {
66 break;
67 }
68 Integer jlev = levii[jj];
69 if (jlev <= lfilp) {
70 std::vector<Integer> nonZeros;
71 std::vector<Real> nonZeroEntries;
72 nonZeros.reserve(n: uBandSet.size()+1);
73 nonZeroEntries.reserve(n: uBandSet.size()+1);
74 const Real entry = U_(jj,jj);
75 if(entry > QL_EPSILON || entry < -1.0*QL_EPSILON) {
76 nonZeros.push_back(x: jj);
77 nonZeroEntries.push_back(x: entry);
78 }
79 auto iter = uBandSet.begin();
80 auto end = uBandSet.end();
81 for (; iter != end; ++iter) {
82 const Real entry = U_(jj,jj+*iter);
83 if(entry > QL_EPSILON || entry < -1.0*QL_EPSILON) {
84 nonZeros.push_back(x: jj+*iter);
85 nonZeroEntries.push_back(x: entry);
86 }
87 }
88 Real fact = w[jj];
89 if(!nonZeroEntries.empty()) {
90 fact /= nonZeroEntries[0];
91 }
92 for (Size k=0; k<nonZeros.size(); ++k) {
93 const Integer j = nonZeros[k] ;
94 const Integer temp = levs(jj,j) + jlev ;
95 if (levii[j] == 0) {
96 if (temp <= lfilp) {
97 w[j] = - fact*nonZeroEntries[k];
98 levii[j] = temp;
99 }
100 }
101 else {
102 w[j] -= fact*nonZeroEntries[k];
103 levii[j] = std::min(a: levii[j],b: temp);
104 }
105 }
106 w[jj] = fact;
107 }
108 }
109 std::vector<Integer> wNonZeros;
110 std::vector<Real> wNonZeroEntries;
111 wNonZeros.reserve(n: w.size());
112 wNonZeroEntries.reserve(n: w.size());
113 for (Size i=0; i<w.size(); ++i) {
114 const Real entry = w[i];
115 if(entry > QL_EPSILON || entry < -1.0*QL_EPSILON) {
116 wNonZeros.push_back(x: i);
117 wNonZeroEntries.push_back(x: entry);
118 }
119 }
120 std::vector<Integer> leviiNonZeroEntries;
121 leviiNonZeroEntries.reserve(n: levii.size());
122 for (int entry : levii) {
123 if (entry > QL_EPSILON || entry < -1.0 * QL_EPSILON) {
124 leviiNonZeroEntries.push_back(x: entry);
125 }
126 }
127 for (Size k=0; k<wNonZeros.size(); ++k) {
128 Integer j = wNonZeros[k];
129 if (j < ii) {
130 L_(ii,j) = wNonZeroEntries[k];
131 lBandSet.insert(x: ii-j);
132 }
133 else {
134 U_(ii,j) = wNonZeroEntries[k];
135 levs(ii,j) = leviiNonZeroEntries[k];
136 if(j-ii > 0) {
137 uBandSet.insert(x: j-ii);
138 }
139 }
140 }
141 }
142 lBands_.resize(new_size: lBandSet.size());
143 uBands_.resize(new_size: uBandSet.size());
144 std::copy(first: lBandSet.begin(), last: lBandSet.end(), result: lBands_.begin());
145 std::copy(first: uBandSet.begin(), last: uBandSet.end(), result: uBands_.begin());
146 }
147
148 const SparseMatrix& SparseILUPreconditioner::L() const {
149 return L_;
150 }
151
152 const SparseMatrix& SparseILUPreconditioner::U() const {
153 return U_;
154 }
155
156 Array SparseILUPreconditioner::apply(const Array& b) const {
157 return backwardSolve(y: forwardSolve(b));
158 }
159
160 Array SparseILUPreconditioner::forwardSolve(const Array& b) const {
161 Integer n = b.size();
162 Array y(n, 0.0);
163 y[0]=b[0]/L_(0,0);
164 for (Integer i=1; i<=n-1; ++i) {
165 y[i] = b[i]/L_(i,i);
166 for (Integer j=lBands_.size()-1;
167 j>=0 && i-Integer(lBands_[j]) <= i-1; --j) {
168 const Integer k = i-Integer(lBands_[j]);
169 if (k >= 0)
170 y[i]-=L_(i,k)*y[k]/L_(i,i);
171 }
172 }
173 return y;
174 }
175
176 Array SparseILUPreconditioner::backwardSolve(const Array& y) const {
177 Size n = y.size();
178 Array x(n, 0.0);
179 x[n-1] = y[n-1]/U_(n-1,n-1);
180 for (Integer i=n-2; i>=0; --i) {
181 x[i] = y[i]/U_(i,i);
182 for (Size j=0; j<uBands_.size() && i+uBands_[j] <= n-1; ++j) {
183 x[i] -= U_(i,i+uBands_[j])*x[i+uBands_[j]]/U_(i,i);
184 }
185 }
186 return x;
187 }
188
189}
190
191

source code of quantlib/ql/math/matrixutilities/sparseilupreconditioner.cpp