1
\documentclass[12pt]{article}
3
\newcommand{\m}[1]{{\bf{#1}}} % for matrices and vectors
4
\newcommand{\tr}{^{\sf T}} % transpose
12
%-------------------------------------------------------------------------------
14
%-------------------------------------------------------------------------------
16
\title{User Guide for LDL, a concise sparse Cholesky package}
17
\author{Timothy A. Davis\thanks{
18
Dept.~of Computer and Information Science and Engineering,
19
Univ.~of Florida, Gainesville, FL, USA.
20
email: davis@cise.ufl.edu.
21
http://www.cise.ufl.edu/$\sim$davis.
22
This work was supported by the National
23
Science Foundation, under grant CCR-0203270.
24
Portions of the work were done while on sabbatical at Stanford University
25
and Lawrence Berkeley National Laboratory (with funding from Stanford
26
University and the SciDAC program).
33
%-------------------------------------------------------------------------------
35
The {\tt LDL} software package is a set of short, concise routines for
36
factorizing symmetric positive-definite sparse matrices, with some
37
applicability to symmetric indefinite matrices. Its primary purpose is
38
to illustrate much of the basic theory of sparse matrix algorithms in as
39
concise a code as possible, including an elegant method
40
of sparse symmetric factorization that computes the factorization row-by-row
41
but stores it column-by-column. The entire symbolic and numeric factorization
42
consists of less than 50 lines of code. The package is written in C,
43
and includes a MATLAB interface.
45
%-------------------------------------------------------------------------------
47
%-------------------------------------------------------------------------------
49
%-------------------------------------------------------------------------------
51
{\tt LDL} is a set of short, concise routines that compute the $\m{LDL}\tr$
52
factorization of a sparse symmetric matrix $\m{A}$. Its primary purpose is
53
to illustrate much of the basic theory of sparse matrix algorithms in as
54
compact a code as possible, including an elegant method of
55
sparse symmetric factorization (related to \cite{Liu86c,Liu91}).
56
The lower triangular factor $\m{L}$ is computed row-by-row, in contrast to the
57
conventional column-by-column method.
58
Although it does not achieve the same level of performance
59
as methods based on dense matrix kernels
60
(such as \cite{NgPeyton93,RothbergGupta91}),
61
its performance is competitive with column-by-column methods that do not
62
use dense kernels \cite{GeorgeLiu79, GeorgeLiu, GilbertMolerSchreiber}.
64
Section~\ref{Algorithm} gives a brief description of the algorithm
65
used in the symbolic and numeric factorization. A more detailed tutorial-level
66
discussion may be found in \cite{Stewart03}. Details
67
of the concise implementation of this method are given in
68
Section~\ref{Implementation}. Sections~\ref{MATLAB}~and~\ref{C} give an
69
overview of how to use the package in MATLAB and in a stand-alone C program.
71
%-------------------------------------------------------------------------------
74
%-------------------------------------------------------------------------------
76
The underlying numerical algorithm is described below. The $k$th
77
step solves a lower triangular system of dimension $k-1$ to compute the
78
$k$th row of $\m{L}$ and the $d_{kk}$ entry of the diagonal matrix $\m{D}$.
79
Colon notation is used for submatrices. For example,
80
$\m{L}_{k,1:k-1}$ refers to the first $k-1$ columns of
81
the $k$th row of $\m{L}$. Similarly, $\m{L}_{1:k-1,1:k-1}$ refers to
82
the leading $(k-1)$-by-$(k-1)$ submatrix of $\m{L}$.
86
\hspace{2em} \= \hspace{2em} \= \hspace{2em} \= \\
88
($\m{LDL}\tr$ factorization of a $n$-by-$n$ symmetric matrix $\m{A}$)} \\
89
\> {\bf for} $k = 1$ {\bf to} $n$ \\
90
\>\> (step 1) Solve $\m{L}_{1:k-1,1:k-1}\m{y} = \m{A}_{1:k-1,k}$ for $\m{y}$ \\
91
\>\> (step 2) $\m{L}_{k,1:k-1} = (\m{D}_{1:k-1,1:k-1}^{-1} \m{y})\tr$ \\
92
\>\> (step 3) $l_{kk} = 1$ \\
93
\>\> (step 4) $d_{kk} = a_{kk} - \m{L}_{k,1:k-1}\m{y}$ \\
98
The algorithm computes an $\m{LDL}\tr$ factorization without numerical pivoting.
99
It can thus factorize any symmetric positive definite matrix, and any
100
symmetric indefinite matrix whose leading minors are all well-conditioned.
102
When $\m{A}$ and $\m{L}$ are sparse, step 1 of Algorithm~1 requires a
103
triangular solve of the form $\m{Lx}=\m{b}$, where all three terms in
104
the equation are sparse. This is the most costly step of the Algorithm.
105
Steps 2 through 4 are fairly straightforward.
107
Let ${\cal X}$ and ${\cal B}$ refer to the set of indices of nonzero entries
108
in $\m{x}$ and $\m{b}$, respectively, in the lower triangular system
109
$\m{Lx}=\m{b}$. To compute $\m{x}$ efficiently
110
the nonzero pattern ${\cal X}$ must be found first.
111
In the general case when $\m{L}$ is arbitrary \cite{GilbertPeierls88},
113
pattern ${\cal X}$ is the set of nodes reachable via paths in the graph $G_L$
114
from all nodes in the set ${\cal B}$, and where the graph $G_L$ has
115
$n$ nodes and a directed edge $(j,i)$ if and only if $l_{ij}$ is nonzero.
116
To compute the numerical solution to $\m{Lx}=\m{b}$ by accessing the columns of
117
$\m{L}$ one at a time, ${\cal X}$ can be traversed
118
in any topological order of the subgraph of $G_L$ consisting of nodes in
119
${\cal X}$. That is, $x_j$ must be computed before $x_i$ if there is a path
120
from $j$ to $i$ in $G_L$. The natural order ($1, 2, \ldots, n$) is one such
121
ordering, but that requires a costly sort of ${\cal X}$.
122
With a graph traversal and topological sort, the solution of $\m{Lx}=\m{b}$
123
can be computed using Algorithm~2 below.
124
The computation of ${\cal X}$ and $\m{x}$ both take
125
time proportional to the floating-point operation count.
129
\hspace{2em} \= \hspace{2em} \= \hspace{2em} \= \\
131
(Solve $\m{Lx}=\m{b}$, where $\m{L}$ is lower triangular with unit diagonal)} \\
132
\> ${\cal X} = \mbox{Reach}_{G_L} ({\cal B})$ \\
133
\> $\m{x} = \m{b}$ \\
134
\> {\bf for} $i \in {\cal X}$ in any topological order \\
135
\>\> $\m{x}_{i+1:n} = \m{x}_{i+1:n} - \m{L}_{i+1:n,i} x_i$ \\
140
The general result also governs the pattern of $\m{y}$ in Algorithm~1.
141
However, in this case $\m{L}$ arises from a sparse Cholesky factorization,
142
and is governed by the elimination tree \cite{Liu90a}.
143
A general graph traversal is not required.
144
In the elimination tree, the parent of node $i$ is the smallest $j > i$
145
such that $l_{ji}$ is nonzero. Node $i$ has no parent if column $i$ of
146
$\m{L}$ is completely zero below the diagonal; $i$ is a root of the
147
elimination tree in this case. The nonzero pattern of $\m{x}$ is the
148
union of all the nodes on the paths from any node $i$ (where $b_i$ is nonzero) to the
149
root of the elimination tree \cite[Thm 2.4]{Liu86c}. It is referred to here as a tree,
150
but in general it can be a forest.
152
Rather than a general topological sort of the subgraph of $G_L$ consisting
153
nodes reachable from nodes in ${\cal B}$, a simpler
154
tree traversal can be used. First, select any nonzero entry $b_i$
155
and follow the path from $i$ to the root of tree.
156
Nodes along this path are marked and placed in a stack,
157
with $i$ at the top of the
158
stack and the root at the bottom.
159
Repeat for every other nonzero entry in $b_i$, in arbitrary order, but stop
160
just before reaching a marked node (the result can be empty if $i$ is already
161
in the stack). The stack now contains ${\cal X}$, a topological ordering of
162
the nonzero pattern of $\m{x}$, which can be used in Algorithm~2 to solve
163
$\m{Lx}=\m{b}$. The time to compute ${\cal X}$
164
using an elimination tree traversal is much faster than the general graph
165
traversal, taking time proportional to the size of ${\cal X}$ rather than the
166
number of floating-point operations required to compute $\m{x}$.
168
In the $k$th step of the factorization, the set ${\cal X}$ becomes the
169
nonzero pattern of row $k$ of $\m{L}$. This step requires the elimination
170
tree of $\m{L}_{1:k-1,1:k-1}$, and must construct the elimination tree of
171
$\m{L}_{1:k,1:k}$ for step $k+1$. Recall that the parent of $i$ in the
172
tree is the smallest $j$ such that $i < j$ and $l_{ji} \ne 0$.
173
Thus, if any node $i$ already has a parent $j$, then $j$ will remain the
174
parent of $i$ in the elimination trees of all other larger leading submatrices
175
of $\m{L}$, and in the elimination tree of $\m{L}$ itself.
176
If $l_{ki} \ne 0$ and $i$ does not have a parent in the elimination tree of
177
$\m{L}_{1:k-1,1:k-1}$, then the parent of $i$ is $k$
178
in the elimination tree of $\m{L}_{1:k,1:k}$.
179
Node $k$ becomes the parent of any node $i \in {\cal X}$ that does not yet
182
Since Algorithm~2 traverses $\m{L}$ in column order, $\m{L}$ is stored in a
183
conventional sparse column representation. Each column $j$ is stored as a list
184
of nonzero values and their corresponding row indices. When row $k$ is
185
computed, the new entries can be placed at the end of each list. As
186
a by-product of computing $\m{L}$ one row at a time,
187
the columns of $\m{L}$ are computed in a sorted manner. This is a convenient
189
MATLAB requires the columns of its sparse matrices to be sorted, for example.
190
Sorted columns improve the speed of Algorithm~2, since the memory access
191
pattern is more regular. The conventional column-by-column algorithm
192
\cite{GeorgeLiu79,GeorgeLiu} does not produce columns of $\m{L}$ with
195
A simple symbolic pre-analysis can be obtained by repeating the subtree traversals.
196
All that is required to compute the nonzero pattern of
197
the $k$th row of $\m{L}$ is the partially constructed elimination tree
198
and the nonzero pattern of the $k$th column of $\m{A}$. This is computed
199
in time proportional to the size of this set, using the elimination tree
200
traversal. Once constructed, the number of nonzeros in each column of
201
$\m{L}$ is incremented, for each entry in ${\cal X}$, and then ${\cal X}$
202
is discarded. The set ${\cal X}$ need not be constructed in topological
203
order, so no stack is required. The run time of the symbolic analysis
204
algorithm is thus proportional to the number of nonzeros in $\m{L}$.
205
This is more costly than the optimal algorithm \cite{GilbertNgPeyton94},
206
which takes time essentially proportional to the number of nonzeros in $\m{A}$.
207
The memory requirements are just the matrix $\m{A}$ and a few size-$n$ integer
208
arrays. The result of the algorithm is the elimination tree, a count
209
of the number of nonzeros in each column of $\m{L}$, and
210
the cumulative sum of the column counts.
212
%-------------------------------------------------------------------------------
213
\section{Implementation}
214
\label{Implementation}
215
%-------------------------------------------------------------------------------
217
Because of its simplicity, the implementation of this algorithm leads to
218
a very short, concise code. The symbolic analysis routine {\tt ldl\_symbolic}
219
shown in Figure~\ref{ldlsymbolic}
220
consists of only 18 lines of executable C code.
221
This includes 5 lines of code to allow for a
222
sparsity-preserving ordering $\m{P}$ so that either $\m{A}$ or $\m{PAP}\tr$
223
can be analyzed, 3 lines of code to compute the cumulative sum of
224
the column counts, and one line of code to speed up a {\tt for} loop.
225
An additional line of code allows for a more general form of the input
226
sparse matrix $\m{A}$.
228
The {\tt n}-by-{\tt n} sparse matrix $\m{A}$ is provided in compressed column
229
form as an {\tt int} array {\tt Ap} of length {\tt n+1},
230
an {\tt int} array {\tt Ai} of length {\tt nz},
231
and a {\tt double} array {\tt Ax} also of length {\tt nz},
232
where {\tt nz} is the number of entries in the matrix.
233
The numerical values of entries in column $j$ are stored in
234
{\tt Ax[Ap[j]} $\ldots$ {\tt Ap[j+1]-1]}
235
and the corresponding row indices are in
236
{\tt Ai[Ap[j]} $\ldots$ {\tt Ap[j+1]-1]}.
237
With {\tt Ap[0] = 0}, the number of entries in the matrix is {\tt nz = Ap[n]}.
238
If no fill-reducing ordering {\tt P} is provided,
239
only entries in the upper triangular part of $\m{A}$ are considered.
240
If {\tt P} is provided and row/column {\tt i} of the
241
matrix $\m{A}$ is the {\tt k}-th row/column of $\m{PAP}\tr$, then {\tt P[k]=i}.
242
Only entries in the upper
243
triangular part of $\m{PAP}\tr$ are considered. These entries may be
244
in the lower triangular part of $\m{A}$, so to ensure that the correct matrix
245
is factorized, all entries of $\m{A}$ should be provided when using the
246
permutation input {\tt P}.
248
The outputs of {\tt ldl\_symbolic} are three size-{\tt n} arrays:
249
{\tt Parent} holds the elimination tree,
250
{\tt Lnz} holds the counts of the number of entries in each column of
252
{\tt Lp} holds the cumulative sum of {\tt Lnz}.
253
The size-{\tt n} array {\tt Flag} is used as workspace.
254
None of the output or workspace arrays need to be initialized.
257
\caption{{\tt ldl\_symbolic:} finding the elimination tree and column counts}
263
int n, /* A and L are n-by-n, where n >= 0 */
264
int Ap [ ], /* input of size n+1, not modified */
265
int Ai [ ], /* input of size nz=Ap[n], not modified */
266
int Lp [ ], /* output of size n+1, not defined on input */
267
int Parent [ ], /* output of size n, not defined on input */
268
int Lnz [ ], /* output of size n, not defined on input */
269
int Flag [ ], /* workspace of size n, not defn. on input or output */
270
int P [ ], /* optional input of size n */
271
int Pinv [ ] /* optional output of size n (used if P is not NULL) */
274
int i, k, p, kk, p2 ;
277
/* If P is present then compute Pinv, the inverse of P */
278
for (k = 0 ; k < n ; k++)
283
for (k = 0 ; k < n ; k++)
285
/* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
286
Parent [k] = -1 ; /* parent of k is not yet known */
287
Flag [k] = k ; /* mark node k as visited */
288
Lnz [k] = 0 ; /* count of nonzeros in column k of L */
289
kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */
291
for (p = Ap [kk] ; p < p2 ; p++)
293
/* A (i,k) is nonzero (original or permuted A) */
294
i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ;
297
/* follow path from i to root of etree, stop at flagged node */
298
for ( ; Flag [i] != k ; i = Parent [i])
300
/* find parent of i if not yet determined */
301
if (Parent [i] == -1) Parent [i] = k ;
302
Lnz [i]++ ; /* L (k,i) is nonzero */
303
Flag [i] = k ; /* mark i as visited */
308
/* construct Lp index array from Lnz column counts */
310
for (k = 0 ; k < n ; k++)
312
Lp [k+1] = Lp [k] + Lnz [k] ;
319
The {\tt ldl\_numeric} numeric factorization routine shown
320
in Figure~\ref{ldlnumeric} consists of only 31 lines of
321
executable code. It includes this same subtree traversal algorithm
322
as {\tt ldl\_symbolic},
323
except that each path is placed on a stack that holds
324
nonzero pattern of the $k$th row of $\m{L}$.
325
This traversal is followed by a sparse forward solve
326
using this pattern, and all of the nonzero entries in
327
the resulting $k$th row of $\m{L}$ are appended to their respective columns
328
in the data structure of $\m{L}$.
331
\caption{{\tt ldl\_numeric:} numeric factorization}
335
int ldl_numeric /* returns n if successful, k if D (k,k) is zero */
337
int n, /* A and L are n-by-n, where n >= 0 */
338
int Ap [ ], /* input of size n+1, not modified */
339
int Ai [ ], /* input of size nz=Ap[n], not modified */
340
double Ax [ ], /* input of size nz=Ap[n], not modified */
341
int Lp [ ], /* input of size n+1, not modified */
342
int Parent [ ], /* input of size n, not modified */
343
int Lnz [ ], /* output of size n, not defn. on input */
344
int Li [ ], /* output of size lnz=Lp[n], not defined on input */
345
double Lx [ ], /* output of size lnz=Lp[n], not defined on input */
346
double D [ ], /* output of size n, not defined on input */
347
double Y [ ], /* workspace of size n, not defn. on input or output */
348
int Pattern [ ], /* workspace of size n, not defn. on input or output */
349
int Flag [ ], /* workspace of size n, not defn. on input or output */
350
int P [ ], /* optional input of size n */
351
int Pinv [ ] /* optional input of size n */
355
int i, k, p, kk, p2, len, top ;
356
for (k = 0 ; k < n ; k++)
358
/* compute nonzero Pattern of kth row of L, in topological order */
359
Y [k] = 0.0 ; /* Y(0:k) is now all zero */
360
top = n ; /* stack for pattern is empty */
361
Flag [k] = k ; /* mark node k as visited */
362
Lnz [k] = 0 ; /* count of nonzeros in column k of L */
363
kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */
365
for (p = Ap [kk] ; p < p2 ; p++)
367
i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ; /* get A(i,k) */
370
Y [i] += Ax [p] ; /* scatter A(i,k) into Y (sum duplicates) */
371
for (len = 0 ; Flag [i] != k ; i = Parent [i])
373
Pattern [len++] = i ; /* L(k,i) is nonzero */
374
Flag [i] = k ; /* mark i as visited */
376
while (len > 0) Pattern [--top] = Pattern [--len] ;
379
/* compute numerical values kth row of L (a sparse triangular solve) */
380
D [k] = Y [k] ; /* get D(k,k) and clear Y(k) */
382
for ( ; top < n ; top++)
384
i = Pattern [top] ; /* Pattern [top:n-1] is pattern of L(:,k) */
385
yi = Y [i] ; /* get and clear Y(i) */
387
p2 = Lp [i] + Lnz [i] ;
388
for (p = Lp [i] ; p < p2 ; p++)
390
Y [Li [p]] -= Lx [p] * yi ;
392
l_ki = yi / D [i] ; /* the nonzero entry L(k,i) */
394
Li [p] = k ; /* store L(k,i) in column form of L */
396
Lnz [i]++ ; /* increment count of nonzeros in col i */
398
if (D [k] == 0.0) return (k) ; /* failure, D(k,k) is zero */
400
return (n) ; /* success, diagonal of D is all nonzero */
406
After the matrix is factorized, the {\tt ldl\_lsolve}, {\tt ldl\_dsolve},
407
and {\tt ldl\_ltsolve} routines shown in Figure~\ref{ldlsolve}
408
are provided to solve
409
$\m{Lx}=\m{b}$, $\m{Dx}=\m{b}$, and $\m{L}\tr\m{x}=\m{b}$, respectively.
410
Together, they solve $\m{Ax}=\m{b}$, and consist of only 10 lines of executable
411
code. If a fill-reducing permutation is used,
412
{\tt ldl\_perm} and {\tt ldl\_permt} must be used to permute $\m{b}$ and
416
\caption{Solve routines}
422
int n, /* L is n-by-n, where n >= 0 */
423
double X [ ], /* size n. right-hand-side on input, soln. on output */
424
int Lp [ ], /* input of size n+1, not modified */
425
int Li [ ], /* input of size lnz=Lp[n], not modified */
426
double Lx [ ] /* input of size lnz=Lp[n], not modified */
430
for (j = 0 ; j < n ; j++)
433
for (p = Lp [j] ; p < p2 ; p++)
435
X [Li [p]] -= Lx [p] * X [j] ;
442
int n, /* D is n-by-n, where n >= 0 */
443
double X [ ], /* size n. right-hand-side on input, soln. on output */
444
double D [ ] /* input of size n, not modified */
448
for (j = 0 ; j < n ; j++)
456
int n, /* L is n-by-n, where n >= 0 */
457
double X [ ], /* size n. right-hand-side on input, soln. on output */
458
int Lp [ ], /* input of size n+1, not modified */
459
int Li [ ], /* input of size lnz=Lp[n], not modified */
460
double Lx [ ] /* input of size lnz=Lp[n], not modified */
464
for (j = n-1 ; j >= 0 ; j--)
467
for (p = Lp [j] ; p < p2 ; p++)
469
X [j] -= Lx [p] * X [Li [p]] ;
477
In addition to appearing as a Collected Algorithm of the ACM \cite{Davis05},
478
{\tt LDL} is available at http://www.cise.ufl.edu/research/sparse.
480
%-------------------------------------------------------------------------------
481
\section{Using LDL in MATLAB}
483
%-------------------------------------------------------------------------------
485
The simplest way to use {\tt LDL} is within MATLAB. Once the {\tt LDL}
486
mexFunction is compiled and installed, the MATLAB statement
487
{\tt [L, D, Parent, fl] = ldl (A)} returns the sparse factorization
488
{\tt A = (L+I)*D*(L+I)'}, where {\tt L} is lower triangular, {\tt D} is a
489
diagonal matrix, and {\tt I} is the {\tt n}-by-{\tt n}
490
identity matrix ({\tt LDL} does not return the unit diagonal of {\tt L}).
491
The elimination tree is returned in {\tt Parent}.
492
If no zero on the diagonal of {\tt D} is encountered, {\tt fl} is the
493
floating-point operation count. Otherwise, {\tt D(-fl,-fl)} is the first
494
zero entry encountered. Let {\tt d=-fl}. The function returns the
495
factorization of {\tt A (1:d,1:d)}, where rows {\tt d+1} to {\tt n} of {\tt L}
496
and {\tt D} are all zero. If a sparsity preserving permutation {\tt P} is
497
passed, {\tt [L, D, Parent, fl] = ldl (A,P)} operates on {\tt A(P,P)} without
498
forming it explicitly.
500
The statement {\tt x = ldl (A, [ ], b)} is roughly equivalent to
501
{\tt x = A}$\backslash${\tt b}, when {\tt A} is sparse, real, and symmetric.
502
The $\m{LDL}\tr$ factorization of {\tt A} is performed. If {\tt P} is
503
provided, {\tt x = ldl (A, P, b)} still performs
504
{\tt x = A}$\backslash${\tt b}, except that {\tt A(P,P)} is factorized
507
%-------------------------------------------------------------------------------
508
\section{Using LDL in a C program}
510
%-------------------------------------------------------------------------------
512
The C-callable {\tt LDL} library consists of nine user-callable routines
513
and one include file.
516
\item {\tt ldl\_symbolic}: given the nonzero pattern of a sparse symmetric
517
matrix $\m{A}$ and an optional permutation $\m{P}$, analyzes either
518
$\m{A}$ or $\m{PAP}\tr$, and returns the elimination tree, the
519
number of nonzeros in each column of $\m{L}$, and the {\tt Lp} array
520
for the sparse matrix data structure for $\m{L}$.
521
Duplicate entries are allowed in the columns of $\m{A}$, and the
522
row indices in each column need not be sorted.
523
Providing a sparsity-preserving ordering is critical for obtaining
524
good performance. A minimum degree ordering
525
(such as AMD \cite{AmestoyDavisDuff96,AmestoyDavisDuff03})
526
or a graph-partitioning based ordering are appropriate.
527
\item {\tt ldl\_numeric}: given {\tt Lp} and the elimination tree computed
528
by {\tt ldl\_symbolic}, and an optional permutation $\m{P}$,
529
returns the numerical factorization of $\m{A}$ or $\m{PAP}\tr$.
530
Duplicate entries are allowed in the columns of $\m{A}$
531
(any duplicate entries are summed), and the
532
row indices in each column need not be sorted.
533
The data structure for $\m{L}$ is the same as $\m{A}$, except that
534
no duplicates appear, and each column has sorted row indices.
535
\item {\tt ldl\_lsolve}: given the factor $\m{L}$ computed by
536
{\tt ldl\_numeric}, solves the linear system $\m{Lx}=\m{b}$, where
537
$\m{x}$ and $\m{b}$ are full $n$-by-1 vectors.
538
\item {\tt ldl\_dsolve}: given the factor $\m{D}$ computed by
539
{\tt ldl\_numeric}, solves the linear system $\m{Dx}=\m{b}$.
540
\item {\tt ldl\_ltsolve}: given the factor $\m{L}$ computed by
541
{\tt ldl\_numeric}, solves the linear system $\m{L}\tr\m{x}=\m{b}$.
542
\item {\tt ldl\_perm}: given a vector $\m{b}$ and a permutation $\m{P}$,
543
returns $\m{x}=\m{Pb}$.
544
\item {\tt ldl\_permt}: given a vector $\m{b}$ and a permutation $\m{P}$,
545
returns $\m{x}=\m{P}\tr\m{b}$.
546
\item {\tt ldl\_valid\_perm}: Except for checking if the diagonal of
547
$\m{D}$ is zero, none of the above routines check their inputs for errors.
548
This routine checks the validity of a permutation $\m{P}$.
549
\item {\tt ldl\_valid\_matrix}: checks if a matrix $\m{A}$ is valid as input
550
to {\tt ldl\_symbolic} and {\tt ldl\_numeric}.
553
Note that the primary input to the {\tt ldl\_symbolic} and
554
{\tt ldl\_numeric} is the sparse matrix $\m{A}$. It is provided in
555
column-oriented form, and only the upper triangular part is accessed.
556
This is slightly different than the primary output: the matrix $\m{L}$,
557
which is lower triangular in column-oriented form.
558
If you wish to factorize a symmetric matrix $\m{A}$ for which only the lower
559
triangular part is supplied, you would need to transpose $\m{A}$ before
560
passing it {\tt ldl\_symbolic} and {\tt ldl\_numeric}.
563
\caption{Example of use}
569
#define N 10 /* A is 10-by-10 */
570
#define ANZ 19 /* # of nonzeros on diagonal and upper triangular part of A */
571
#define LNZ 13 /* # of nonzeros below the diagonal of L */
573
int main (int argc, int **argv)
575
/* only the upper triangular part of A is required */
576
int Ap [N+1] = {0, 1, 2, 3, 4, 6, 7, 9, 11, 15, ANZ},
577
Ai [ANZ] = {0, 1, 2, 3, 1,4, 5, 4,6, 4,7, 0,4,7,8, 1,4,6,9 } ;
578
double Ax [ANZ] = {1.7, 1., 1.5, 1.1, .02,2.6, 1.2, .16,1.3, .09,1.6,
579
.13,.52,.11,1.4, .01,.53,.56,3.1},
580
b [N] = {.287, .22, .45, .44, 2.486, .72, 1.55, 1.424, 1.621, 3.759};
581
double Lx [LNZ], D [N], Y [N] ;
582
int Li [LNZ], Lp [N+1], Parent [N], Lnz [N], Flag [N], Pattern [N], d, i ;
584
/* factorize A into LDL' (P and Pinv not used) */
585
ldl_symbolic (N, Ap, Ai, Lp, Parent, Lnz, Flag, NULL, NULL) ;
586
printf ("Nonzeros in L, excluding diagonal: %d\n", Lp [N]) ;
587
d = ldl_numeric (N, Ap, Ai, Ax, Lp, Parent, Lnz, Li, Lx, D, Y, Pattern,
592
/* solve Ax=b, overwriting b with the solution x */
593
ldl_lsolve (N, b, Lp, Li, Lx) ;
594
ldl_dsolve (N, b, D) ;
595
ldl_ltsolve (N, b, Lp, Li, Lx) ;
596
for (i = 0 ; i < N ; i++) printf ("x [%d] = %g\n", i, b [i]) ;
600
printf ("ldl_numeric failed, D (%d,%d) is zero\n", d, d) ;
608
The program in Figure~\ref{ldlsimple}
609
illustrates the basic usage of the {\tt LDL} routines.
610
It analyzes and factorizes the sparse symmetric positive-definite matrix
614
\begin{array}{cccccccccc}
615
1.7 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & .13 & 0 \\
616
0 & 1. & 0 & 0 & .02 & 0 & 0 & 0 & 0 & .01 \\
617
0 & 0 & 1.5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
618
0 & 0 & 0 & 1.1 & 0 & 0 & 0 & 0 & 0 & 0 \\
619
0 & .02 & 0 & 0 & 2.6 & 0 & .16 & .09 & .52 & .53 \\
620
0 & 0 & 0 & 0 & 0 & 1.2 & 0 & 0 & 0 & 0 \\
621
0 & 0 & 0 & 0 & .16 & 0 & 1.3 & 0 & 0 & .56 \\
622
0 & 0 & 0 & 0 & .09 & 0 & 0 & 1.6 & .11 & 0 \\
623
.13 & 0 & 0 & 0 & .52 & 0 & 0 & .11 & 1.4 & 0 \\
624
0 & .01 & 0 & 0 & .53 & 0 & .56 & 0 & 0 & 3.1 \\
629
and then solves a system $\m{Ax}=\m{b}$ whose true solution is
630
$x_i = i/10$. Note that {\tt Li} and {\tt Lx} are statically allocated.
631
Normally they would be allocated after their size, {\tt Lp[n]},
632
is determined by {\tt ldl\_symbolic}.
633
More example programs are included with the {\tt LDL} package.
635
\section{Acknowledgments}
637
I would like to thank Pete Stewart for his comments on an earlier draft
638
of this software and its accompanying paper.
641
\bibliographystyle{plain}
added
removed
LDL/ldl_userguide.tex
