/usr/share/doc/octave/octave.html/Iterative-Techniques.html

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- Created by GNU Texinfo 6.5, http://www.gnu.org/software/texinfo/ -->
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<title>Iterative Techniques (GNU Octave)</title>

<meta name="description" content="Iterative Techniques (GNU Octave)">
<meta name="keywords" content="Iterative Techniques (GNU Octave)">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<link href="index.html#Top" rel="start" title="Top">
<link href="Concept-Index.html#Concept-Index" rel="index" title="Concept Index">
<link href="index.html#SEC_Contents" rel="contents" title="Table of Contents">
<link href="Sparse-Matrices.html#Sparse-Matrices" rel="up" title="Sparse Matrices">
<link href="Real-Life-Example.html#Real-Life-Example" rel="next" title="Real Life Example">
<link href="Sparse-Linear-Algebra.html#Sparse-Linear-Algebra" rel="prev" title="Sparse Linear Algebra">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.indentedblock {margin-right: 0em}
blockquote.smallindentedblock {margin-right: 0em; font-size: smaller}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smalllisp {margin-left: 3.2em}
kbd {font-style: oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nolinebreak {white-space: nowrap}
span.roman {font-family: initial; font-weight: normal}
span.sansserif {font-family: sans-serif; font-weight: normal}
ul.no-bullet {list-style: none}
-->
</style>
<link rel="stylesheet" type="text/css" href="octave.css">


</head>

<body lang="en">
<a name="Iterative-Techniques"></a>
<div class="header">
<p>
Next: <a href="Real-Life-Example.html#Real-Life-Example" accesskey="n" rel="next">Real Life Example</a>, Previous: <a href="Sparse-Linear-Algebra.html#Sparse-Linear-Algebra" accesskey="p" rel="prev">Sparse Linear Algebra</a>, Up: <a href="Sparse-Matrices.html#Sparse-Matrices" accesskey="u" rel="up">Sparse Matrices</a> &nbsp; [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html#Concept-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="Iterative-Techniques-Applied-to-Sparse-Matrices"></a>
<h3 class="section">22.3 Iterative Techniques Applied to Sparse Matrices</h3>

<p>The left division <code>\</code> and right division <code>/</code> operators,
discussed in the previous section, use direct solvers to resolve a
linear equation of the form <code><var>x</var> = <var>A</var> \ <var>b</var></code> or
<code><var>x</var> = <var>b</var> / <var>A</var></code>.  Octave also includes a number of
functions to solve sparse linear equations using iterative techniques.
</p>
<a name="XREFpcg"></a><dl>
<dt><a name="index-pcg"></a>: <em><var>x</var> =</em> <strong>pcg</strong> <em>(<var>A</var>, <var>b</var>, <var>tol</var>, <var>maxit</var>, <var>m1</var>, <var>m2</var>, <var>x0</var>, &hellip;)</em></dt>
<dt><a name="index-pcg-1"></a>: <em>[<var>x</var>, <var>flag</var>, <var>relres</var>, <var>iter</var>, <var>resvec</var>, <var>eigest</var>] =</em> <strong>pcg</strong> <em>(&hellip;)</em></dt>
<dd>
<p>Solve the linear system of equations <code><var>A</var>&nbsp;*&nbsp;<var>x</var>&nbsp;=&nbsp;<var>b</var></code><!-- /@w -->
by means of the Preconditioned Conjugate Gradient iterative method.
</p>
<p>The input arguments are
</p>
<ul>
<li> <var>A</var> can be either a square (preferably sparse) matrix or a function
handle, inline function or string containing the name of a function which
computes <code><var>A</var>&nbsp;*&nbsp;<var>x</var></code><!-- /@w -->.  In principle, <var>A</var> should be
symmetric and positive definite; if <code>pcg</code> finds <var>A</var> not to be
positive definite, a warning is printed and the <var>flag</var> output will be
set.

</li><li> <var>b</var> is the right-hand side vector.

</li><li> <var>tol</var> is the required relative tolerance for the residual error,
<code><var>b</var>&nbsp;<span class="nolinebreak">-</span>&nbsp;<var>A</var>&nbsp;*&nbsp;<var>x</var></code><!-- /@w -->.  The iteration stops if
<code>norm&nbsp;(<var>b</var>&nbsp;<span class="nolinebreak">-</span>&nbsp;<var>A</var>&nbsp;*&nbsp;<var>x</var>)</code>&nbsp;&le;&nbsp;<var>tol</var>&nbsp;*&nbsp;norm&nbsp;(<var>b</var>)<!-- /@w --><!-- /@w -->.
If <var>tol</var> is omitted or empty then a tolerance of 1e-6 is used.

</li><li> <var>maxit</var> is the maximum allowable number of iterations; if <var>maxit</var>
is omitted or empty then a value of 20 is used.

</li><li> <var>m</var> = <var>m1</var> * <var>m2</var> is the (left) preconditioning matrix, so that
the iteration is (theoretically) equivalent to solving by <code>pcg</code>
<code><var>P</var>&nbsp;*&nbsp;<var>x</var>&nbsp;=&nbsp;<var>m</var>&nbsp;\&nbsp;<var>b</var></code><!-- /@w -->, with
<code><var>P</var>&nbsp;=&nbsp;<var>m</var>&nbsp;\&nbsp;<var>A</var></code><!-- /@w -->.
Note that a proper choice of the preconditioner may dramatically improve
the overall performance of the method.  Instead of matrices <var>m1</var> and
<var>m2</var>, the user may pass two functions which return the results of
applying the inverse of <var>m1</var> and <var>m2</var> to a vector (usually this is
the preferred way of using the preconditioner).  If <var>m1</var> is omitted or
empty <code>[]</code> then no preconditioning is applied.  If <var>m2</var> is
omitted, <var>m</var> = <var>m1</var> will be used as a preconditioner.

</li><li> <var>x0</var> is the initial guess.  If <var>x0</var> is omitted or empty then the
function sets <var>x0</var> to a zero vector by default.
</li></ul>

<p>The arguments which follow <var>x0</var> are treated as parameters, and passed in
a proper way to any of the functions (<var>A</var> or <var>m</var>) which are passed
to <code>pcg</code>.  See the examples below for further details.  The output
arguments are
</p>
<ul>
<li> <var>x</var> is the computed approximation to the solution of
<code><var>A</var>&nbsp;*&nbsp;<var>x</var>&nbsp;=&nbsp;<var>b</var></code><!-- /@w -->.

</li><li> <var>flag</var> reports on the convergence.  A value of 0 means the solution
converged and the tolerance criterion given by <var>tol</var> is satisfied.
A value of 1 means that the <var>maxit</var> limit for the iteration count was
reached.  A value of 3 indicates that the (preconditioned) matrix was found
not to be positive definite.

</li><li> <var>relres</var> is the ratio of the final residual to its initial value,
measured in the Euclidean norm.

</li><li> <var>iter</var> is the actual number of iterations performed.

</li><li> <var>resvec</var> describes the convergence history of the method.
<code><var>resvec</var>(i,1)</code> is the Euclidean norm of the residual, and
<code><var>resvec</var>(i,2)</code> is the preconditioned residual norm, after the
(<var>i</var>-1)-th iteration, <code><var>i</var> = 1, 2, &hellip;, <var>iter</var>+1</code>.
The preconditioned residual norm is defined as
<code>norm (<var>r</var>) ^ 2 = <var>r</var>' * (<var>m</var> \ <var>r</var>)</code> where
<code><var>r</var> = <var>b</var> - <var>A</var> * <var>x</var></code>, see also the
description of <var>m</var>.  If <var>eigest</var> is not required, only
<code><var>resvec</var>(:,1)</code> is returned.

</li><li> <var>eigest</var> returns the estimate for the smallest <code><var>eigest</var>(1)</code>
and largest <code><var>eigest</var>(2)</code> eigenvalues of the preconditioned matrix
<code><var>P</var>&nbsp;=&nbsp;<var>m</var>&nbsp;\&nbsp;<var>A</var></code><!-- /@w -->.  In particular, if no
preconditioning is used, the estimates for the extreme eigenvalues of
<var>A</var> are returned.  <code><var>eigest</var>(1)</code> is an overestimate and
<code><var>eigest</var>(2)</code> is an underestimate, so that
<code><var>eigest</var>(2) / <var>eigest</var>(1)</code> is a lower bound for
<code>cond (<var>P</var>, 2)</code>, which nevertheless in the limit should
theoretically be equal to the actual value of the condition number.
The method which computes <var>eigest</var> works only for symmetric positive
definite <var>A</var> and <var>m</var>, and the user is responsible for verifying this
assumption.
</li></ul>

<p>Let us consider a trivial problem with a diagonal matrix (we exploit the
sparsity of A)
</p>
<div class="example">
<pre class="example">n = 10;
A = diag (sparse (1:n));
b = rand (n, 1);
[l, u, p] = ilu (A, struct (&quot;droptol&quot;, 1.e-3));
</pre></div>

<p><small>EXAMPLE 1:</small> Simplest use of <code>pcg</code>
</p>
<div class="example">
<pre class="example">x = pcg (A, b)
</pre></div>

<p><small>EXAMPLE 2:</small> <code>pcg</code> with a function which computes
<code><var>A</var> * <var>x</var></code>
</p>
<div class="example">
<pre class="example">function y = apply_a (x)
  y = [1:N]' .* x;
endfunction

x = pcg (&quot;apply_a&quot;, b)
</pre></div>

<p><small>EXAMPLE 3:</small> <code>pcg</code> with a preconditioner: <var>l</var> * <var>u</var>
</p>
<div class="example">
<pre class="example">x = pcg (A, b, 1.e-6, 500, l*u)
</pre></div>

<p><small>EXAMPLE 4:</small> <code>pcg</code> with a preconditioner: <var>l</var> * <var>u</var>.
Faster than <small>EXAMPLE 3</small> since lower and upper triangular matrices are
easier to invert
</p>
<div class="example">
<pre class="example">x = pcg (A, b, 1.e-6, 500, l, u)
</pre></div>

<p><small>EXAMPLE 5:</small> Preconditioned iteration, with full diagnostics.  The
preconditioner (quite strange, because even the original matrix <var>A</var> is
trivial) is defined as a function
</p>
<div class="example">
<pre class="example">function y = apply_m (x)
  k = floor (length (x) - 2);
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec, eigest] = ...
                   pcg (A, b, [], [], &quot;apply_m&quot;);
semilogy (1:iter+1, resvec);
</pre></div>

<p><small>EXAMPLE 6:</small> Finally, a preconditioner which depends on a parameter
<var>k</var>.
</p>
<div class="example">
<pre class="example">function y = apply_M (x, varargin)
  K = varargin{1};
  y = x;
  y(1:K) = x(1:K) ./ [1:K]';
endfunction

[x, flag, relres, iter, resvec, eigest] = ...
     pcg (A, b, [], [], &quot;apply_m&quot;, [], [], 3)
</pre></div>

<p>References:
</p>
<ol>
<li> C.T. Kelley, <cite>Iterative Methods for Linear and Nonlinear Equations</cite>,
SIAM, 1995. (the base PCG algorithm)

</li><li> Y. Saad, <cite>Iterative Methods for Sparse Linear Systems</cite>,
PWS 1996. (condition number estimate from PCG)
Revised version of this book is available online at
<a href="http://www-users.cs.umn.edu/~saad/books.html">http://www-users.cs.umn.edu/~saad/books.html</a>
</li></ol>


<p><strong>See also:</strong> <a href="Creating-Sparse-Matrices.html#XREFsparse">sparse</a>, <a href="#XREFpcr">pcr</a>.
</p></dd></dl>


<a name="XREFpcr"></a><dl>
<dt><a name="index-pcr"></a>: <em><var>x</var> =</em> <strong>pcr</strong> <em>(<var>A</var>, <var>b</var>, <var>tol</var>, <var>maxit</var>, <var>m</var>, <var>x0</var>, &hellip;)</em></dt>
<dt><a name="index-pcr-1"></a>: <em>[<var>x</var>, <var>flag</var>, <var>relres</var>, <var>iter</var>, <var>resvec</var>] =</em> <strong>pcr</strong> <em>(&hellip;)</em></dt>
<dd>
<p>Solve the linear system of equations <code><var>A</var> * <var>x</var> = <var>b</var></code> by
means of the Preconditioned Conjugate Residuals iterative method.
</p>
<p>The input arguments are
</p>
<ul>
<li> <var>A</var> can be either a square (preferably sparse) matrix or a function
handle, inline function or string containing the name of a function which
computes <code><var>A</var> * <var>x</var></code>.  In principle <var>A</var> should be
symmetric and non-singular; if <code>pcr</code> finds <var>A</var> to be numerically
singular, you will get a warning message and the <var>flag</var> output
parameter will be set.

</li><li> <var>b</var> is the right hand side vector.

</li><li> <var>tol</var> is the required relative tolerance for the residual error,
<code><var>b</var> - <var>A</var> * <var>x</var></code>.  The iteration stops if
<code>norm (<var>b</var> - <var>A</var> * <var>x</var>) &lt;=
<var>tol</var> * norm (<var>b</var> - <var>A</var> * <var>x0</var>)</code>.
If <var>tol</var> is empty or is omitted, the function sets
<code><var>tol</var> = 1e-6</code> by default.

</li><li> <var>maxit</var> is the maximum allowable number of iterations; if <code>[]</code> is
supplied for <code>maxit</code>, or <code>pcr</code> has less arguments, a default
value equal to 20 is used.

</li><li> <var>m</var> is the (left) preconditioning matrix, so that the iteration is
(theoretically) equivalent to solving by
<code>pcr</code> <code><var>P</var> * <var>x</var> = <var>m</var> \ <var>b</var></code>, with
<code><var>P</var> = <var>m</var> \ <var>A</var></code>.  Note that a proper choice of the
preconditioner may dramatically improve the overall performance of the
method.  Instead of matrix <var>m</var>, the user may pass a function which
returns the results of applying the inverse of <var>m</var> to a vector
(usually this is the preferred way of using the preconditioner).  If
<code>[]</code> is supplied for <var>m</var>, or <var>m</var> is omitted, no
preconditioning is applied.

</li><li> <var>x0</var> is the initial guess.  If <var>x0</var> is empty or omitted, the
function sets <var>x0</var> to a zero vector by default.
</li></ul>

<p>The arguments which follow <var>x0</var> are treated as parameters, and passed
in a proper way to any of the functions (<var>A</var> or <var>m</var>) which are
passed to <code>pcr</code>.  See the examples below for further details.
</p>
<p>The output arguments are
</p>
<ul>
<li> <var>x</var> is the computed approximation to the solution of
<code><var>A</var> * <var>x</var> = <var>b</var></code>.

</li><li> <var>flag</var> reports on the convergence.  <code><var>flag</var> = 0</code> means the
solution converged and the tolerance criterion given by <var>tol</var> is
satisfied.  <code><var>flag</var> = 1</code> means that the <var>maxit</var> limit for the
iteration count was reached.  <code><var>flag</var> = 3</code> reports a <code>pcr</code>
breakdown, see [1] for details.

</li><li> <var>relres</var> is the ratio of the final residual to its initial value,
measured in the Euclidean norm.

</li><li> <var>iter</var> is the actual number of iterations performed.

</li><li> <var>resvec</var> describes the convergence history of the method, so that
<code><var>resvec</var> (i)</code> contains the Euclidean norms of the residual after
the (<var>i</var>-1)-th iteration, <code><var>i</var> = 1,2, &hellip;, <var>iter</var>+1</code>.
</li></ul>

<p>Let us consider a trivial problem with a diagonal matrix (we exploit the
sparsity of A)
</p>
<div class="example">
<pre class="example">n = 10;
A = sparse (diag (1:n));
b = rand (N, 1);
</pre></div>

<p><small>EXAMPLE 1:</small> Simplest use of <code>pcr</code>
</p>
<div class="example">
<pre class="example">x = pcr (A, b)
</pre></div>

<p><small>EXAMPLE 2:</small> <code>pcr</code> with a function which computes
<code><var>A</var> * <var>x</var></code>.
</p>
<div class="example">
<pre class="example">function y = apply_a (x)
  y = [1:10]' .* x;
endfunction

x = pcr (&quot;apply_a&quot;, b)
</pre></div>

<p><small>EXAMPLE 3:</small>  Preconditioned iteration, with full diagnostics.  The
preconditioner (quite strange, because even the original matrix
<var>A</var> is trivial) is defined as a function
</p>
<div class="example">
<pre class="example">function y = apply_m (x)
  k = floor (length (x) - 2);
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], &quot;apply_m&quot;)
semilogy ([1:iter+1], resvec);
</pre></div>

<p><small>EXAMPLE 4:</small> Finally, a preconditioner which depends on a
parameter <var>k</var>.
</p>
<div class="example">
<pre class="example">function y = apply_m (x, varargin)
  k = varargin{1};
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], &quot;apply_m&quot;', [], 3)
</pre></div>

<p>References:
</p>
<p>[1] W. Hackbusch, <cite>Iterative Solution of Large Sparse
Systems of Equations</cite>, section 9.5.4; Springer, 1994
</p>

<p><strong>See also:</strong> <a href="Creating-Sparse-Matrices.html#XREFsparse">sparse</a>, <a href="#XREFpcg">pcg</a>.
</p></dd></dl>


<p>The speed with which an iterative solver converges to a solution can be
accelerated with the use of a pre-conditioning matrix <var>M</var>.  In this
case the linear equation <code><var>M</var>^-1 * <var>x</var> = <var>M</var>^-1 *
<var>A</var> \ <var>b</var></code> is solved instead.  Typical pre-conditioning matrices
are partial factorizations of the original matrix.
</p>
<a name="XREFichol"></a><dl>
<dt><a name="index-ichol"></a>: <em><var>L</var> =</em> <strong>ichol</strong> <em>(<var>A</var>)</em></dt>
<dt><a name="index-ichol-1"></a>: <em><var>L</var> =</em> <strong>ichol</strong> <em>(<var>A</var>, <var>opts</var>)</em></dt>
<dd>
<p>Compute the incomplete Cholesky factorization of the sparse square matrix
<var>A</var>.
</p>
<p>By default, <code>ichol</code> uses only the lower triangle of <var>A</var> and
produces a lower triangular factor <var>L</var> such that <code>L*L'</code>
approximates <var>A</var>.
</p>
<p>The factor given by this routine may be useful as a preconditioner for a
system of linear equations being solved by iterative methods such as
PCG (Preconditioned Conjugate Gradient).
</p>
<p>The factorization may be modified by passing options in a structure
<var>opts</var>.  The option name is a field of the structure and the setting
is the value of field.  Names and specifiers are case sensitive.
</p>
<dl compact="compact">
<dt>type</dt>
<dd><p>Type of factorization.
</p>
<dl compact="compact">
<dt><code>&quot;nofill&quot;</code> (default)</dt>
<dd><p>Incomplete Cholesky factorization with no fill-in (IC(0)).
</p>
</dd>
<dt><code>&quot;ict&quot;</code></dt>
<dd><p>Incomplete Cholesky factorization with threshold dropping (ICT).
</p></dd>
</dl>

</dd>
<dt>diagcomp</dt>
<dd><p>A non-negative scalar <var>alpha</var> for incomplete Cholesky factorization of
<code><var>A</var> + <var>alpha</var> * diag (diag (<var>A</var>))</code> instead of <var>A</var>.
This can be useful when <var>A</var> is not positive definite.  The default value
is 0.
</p>
</dd>
<dt>droptol</dt>
<dd><p>A non-negative scalar specifying the drop tolerance for factorization if
performing ICT.  The default value is 0 which produces the
complete Cholesky factorization.
</p>
<p>Non-diagonal entries of <var>L</var> are set to 0 unless
</p>
<p><code>abs (<var>L</var>(i,j)) &gt;= droptol * norm (<var>A</var>(j:end, j), 1)</code>.
</p>
</dd>
<dt>michol</dt>
<dd><p>Modified incomplete Cholesky factorization:
</p>
<dl compact="compact">
<dt><code>&quot;off&quot;</code> (default)</dt>
<dd><p>Row and column sums are not necessarily preserved.
</p>
</dd>
<dt><code>&quot;on&quot;</code></dt>
<dd><p>The diagonal of <var>L</var> is modified so that row (and column) sums are
preserved even when elements have been dropped during the factorization.
The relationship preserved is: <code><var>A</var> * e = <var>L</var> * <var>L</var>' * e</code>,
where e is a vector of ones.
</p></dd>
</dl>

</dd>
<dt>shape</dt>
<dd>
<dl compact="compact">
<dt><code>&quot;lower&quot;</code> (default)</dt>
<dd><p>Use only the lower triangle of <var>A</var> and return a lower triangular factor
<var>L</var> such that <code>L*L'</code> approximates <var>A</var>.
</p>
</dd>
<dt><code>&quot;upper&quot;</code></dt>
<dd><p>Use only the upper triangle of <var>A</var> and return an upper triangular factor
<var>U</var> such that <code>U'*U</code> approximates <var>A</var>.
</p></dd>
</dl>
</dd>
</dl>

<p>EXAMPLES
</p>
<p>The following problem demonstrates how to factorize a sample symmetric
positive definite matrix with the full Cholesky decomposition and with the
incomplete one.
</p>
<div class="example">
<pre class="example">A = [ 0.37, -0.05,  -0.05,  -0.07;
     -0.05,  0.116,  0.0,   -0.05;
     -0.05,  0.0,    0.116, -0.05;
     -0.07, -0.05,  -0.05,   0.202];
A = sparse (A);
nnz (tril (A))
ans =  9
L = chol (A, &quot;lower&quot;);
nnz (L)
ans =  10
norm (A - L * L', &quot;fro&quot;) / norm (A, &quot;fro&quot;)
ans =  1.1993e-16
opts.type = &quot;nofill&quot;;
L = ichol (A, opts);
nnz (L)
ans =  9
norm (A - L * L', &quot;fro&quot;) / norm (A, &quot;fro&quot;)
ans =  0.019736
</pre></div>

<p>Another example for decomposition is a finite difference matrix used to
solve a boundary value problem on the unit square.
</p>
<div class="example">
<pre class="example">nx = 400; ny = 200;
hx = 1 / (nx + 1); hy = 1 / (ny + 1);
Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
               [-1 0 1 ], nx, nx) / (hx ^ 2);
Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
               [-1 0 1 ], ny, ny) / (hy ^ 2);
A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
nnz (tril (A))
ans =  239400
opts.type = &quot;nofill&quot;;
L = ichol (A, opts);
nnz (tril (A))
ans =  239400
norm (A - L * L', &quot;fro&quot;) / norm (A, &quot;fro&quot;)
ans =  0.062327
</pre></div>

<p>References for implemented algorithms:
</p>
<p>[1] Y. Saad. &quot;Preconditioning Techniques.&quot; <cite>Iterative
Methods for Sparse Linear Systems</cite>, PWS Publishing Company, 1996.
</p>
<p>[2] M. Jones, P. Plassmann: <cite>An Improved Incomplete
Cholesky Factorization</cite>, 1992.
</p>
<p><strong>See also:</strong> <a href="Matrix-Factorizations.html#XREFchol">chol</a>, <a href="#XREFilu">ilu</a>, <a href="#XREFpcg">pcg</a>.
</p></dd></dl>


<a name="XREFilu"></a><dl>
<dt><a name="index-ilu"></a>: <em></em> <strong>ilu</strong> <em>(<var>A</var>)</em></dt>
<dt><a name="index-ilu-1"></a>: <em></em> <strong>ilu</strong> <em>(<var>A</var>, <var>opts</var>)</em></dt>
<dt><a name="index-ilu-2"></a>: <em>[<var>L</var>, <var>U</var>] =</em> <strong>ilu</strong> <em>(&hellip;)</em></dt>
<dt><a name="index-ilu-3"></a>: <em>[<var>L</var>, <var>U</var>, <var>P</var>] =</em> <strong>ilu</strong> <em>(&hellip;)</em></dt>
<dd>
<p>Compute the incomplete LU factorization of the sparse square matrix <var>A</var>.
</p>
<p><code>ilu</code> returns a unit lower triangular matrix <var>L</var>, an upper
triangular matrix <var>U</var>, and optionally a permutation matrix <var>P</var>, such
that <code><var>L</var>*<var>U</var></code> approximates <code><var>P</var>*<var>A</var></code>.
</p>
<p>The factors given by this routine may be useful as preconditioners for a
system of linear equations being solved by iterative methods such as BICG
(BiConjugate Gradients) or GMRES (Generalized Minimum Residual Method).
</p>
<p>The factorization may be modified by passing options in a structure
<var>opts</var>.  The option name is a field of the structure and the setting
is the value of field.  Names and specifiers are case sensitive.
</p>
<dl compact="compact">
<dt><code>type</code></dt>
<dd><p>Type of factorization.
</p>
<dl compact="compact">
<dt><code>&quot;nofill&quot;</code></dt>
<dd><p>ILU factorization with no fill-in (ILU(0)).
</p>
<p>Additional supported options: <code>milu</code>.
</p>
</dd>
<dt><code>&quot;crout&quot;</code></dt>
<dd><p>Crout version of ILU factorization (ILUC).
</p>
<p>Additional supported options: <code>milu</code>, <code>droptol</code>.
</p>
</dd>
<dt><code>&quot;ilutp&quot;</code> (default)</dt>
<dd><p>ILU factorization with threshold and pivoting.
</p>
<p>Additional supported options: <code>milu</code>, <code>droptol</code>, <code>udiag</code>,
<code>thresh</code>.
</p></dd>
</dl>

</dd>
<dt><code>droptol</code></dt>
<dd><p>A non-negative scalar specifying the drop tolerance for factorization.  The
default value is 0 which produces the complete LU factorization.
</p>
<p>Non-diagonal entries of <var>U</var> are set to 0 unless
</p>
<p><code>abs (<var>U</var>(i,j)) &gt;= droptol * norm (<var>A</var>(:,j))</code>.
</p>
<p>Non-diagonal entries of <var>L</var> are set to 0 unless
</p>
<p><code>abs (<var>L</var>(i,j)) &gt;= droptol * norm (<var>A</var>(:,j))/<var>U</var>(j,j)</code>.
</p>
</dd>
<dt><code>milu</code></dt>
<dd><p>Modified incomplete LU factorization:
</p>
<dl compact="compact">
<dt><code>&quot;row&quot;</code></dt>
<dd><p>Row-sum modified incomplete LU factorization.
The factorization preserves row sums:
<code><var>A</var> * e = <var>L</var> * <var>U</var> * e</code>, where e is a vector of ones.
</p>
</dd>
<dt><code>&quot;col&quot;</code></dt>
<dd><p>Column-sum modified incomplete LU factorization.
The factorization preserves column sums:
<code>e' * <var>A</var> = e' * <var>L</var> * <var>U</var></code>.
</p>
</dd>
<dt><code>&quot;off&quot;</code> (default)</dt>
<dd><p>Row and column sums are not necessarily preserved.
</p></dd>
</dl>

</dd>
<dt><code>udiag</code></dt>
<dd><p>If true, any zeros on the diagonal of the upper triangular factor are
replaced by the local drop tolerance
<code>droptol * norm (<var>A</var>(:,j))/<var>U</var>(j,j)</code>.  The default is false.
</p>
</dd>
<dt><code>thresh</code></dt>
<dd><p>Pivot threshold for factorization.  It can range between 0 (diagonal
pivoting) and 1 (default), where the maximum magnitude entry in the column
is chosen to be the pivot.
</p></dd>
</dl>

<p>If <code>ilu</code> is called with just one output, the returned matrix is
<code><var>L</var> + <var>U</var> - speye (size (<var>A</var>))</code>, where <var>L</var> is unit
lower triangular and <var>U</var> is upper triangular.
</p>
<p>With two outputs, <code>ilu</code> returns a unit lower triangular matrix <var>L</var>
and an upper triangular matrix <var>U</var>.  For <var>opts</var>.type ==
<code>&quot;ilutp&quot;</code>, one of the factors is permuted based on the value of
<var>opts</var>.milu.  When <var>opts</var>.milu == <code>&quot;row&quot;</code>, <var>U</var> is a
column permuted upper triangular factor.  Otherwise, <var>L</var> is a
row-permuted unit lower triangular factor.
</p>
<p>If there are three named outputs and <var>opts</var>.milu != <code>&quot;row&quot;</code>,
<var>P</var> is returned such that <var>L</var> and <var>U</var> are incomplete factors
of <code><var>P</var>*<var>A</var></code>.  When <var>opts</var>.milu == <code>&quot;row&quot;</code>, <var>P</var>
is returned such that <var>L</var> and <var>U</var> are incomplete factors of
<code><var>A</var>*<var>P</var></code>.
</p>
<p>EXAMPLES
</p>
<div class="example">
<pre class="example">A = gallery (&quot;neumann&quot;, 1600) + speye (1600);
opts.type = &quot;nofill&quot;;
nnz (A)
ans = 7840

nnz (lu (A))
ans = 126478

nnz (ilu (A, opts))
ans = 7840
</pre></div>

<p>This shows that <var>A</var> has 7,840 nonzeros, the complete LU factorization
has 126,478 nonzeros, and the incomplete LU factorization, with 0 level of
fill-in, has 7,840 nonzeros, the same amount as <var>A</var>.  Taken from:
http://www.mathworks.com/help/matlab/ref/ilu.html
</p>
<div class="example">
<pre class="example">A = gallery (&quot;wathen&quot;, 10, 10);
b = sum (A, 2);
tol = 1e-8;
maxit = 50;
opts.type = &quot;crout&quot;;
opts.droptol = 1e-4;
[L, U] = ilu (A, opts);
x = bicg (A, b, tol, maxit, L, U);
norm (A * x - b, inf)
</pre></div>

<p>This example uses ILU as preconditioner for a random FEM-Matrix, which has a
large condition number.  Without <var>L</var> and <var>U</var> BICG would not
converge.
</p>

<p><strong>See also:</strong> <a href="Matrix-Factorizations.html#XREFlu">lu</a>, <a href="#XREFichol">ichol</a>, <a href="Specialized-Solvers.html#XREFbicg">bicg</a>, <a href="Specialized-Solvers.html#XREFgmres">gmres</a>.
</p></dd></dl>


<hr>
<div class="header">
<p>
Next: <a href="Real-Life-Example.html#Real-Life-Example" accesskey="n" rel="next">Real Life Example</a>, Previous: <a href="Sparse-Linear-Algebra.html#Sparse-Linear-Algebra" accesskey="p" rel="prev">Sparse Linear Algebra</a>, Up: <a href="Sparse-Matrices.html#Sparse-Matrices" accesskey="u" rel="up">Sparse Matrices</a> &nbsp; [<a href="index.html#SEC_Contents" title="Table of contents" rel="contents">Contents</a>][<a href="Concept-Index.html#Concept-Index" title="Index" rel="index">Index</a>]</p>
</div>



</body>
</html>
octave-doc 4.2.2-1ubuntu1 / usr / share / doc / octave / octave.html / Iterative-Techniques.html
