SCM

SCM Repository

[matrix] Annotation of /pkg/src/R_ldl.c
ViewVC logotype

Annotation of /pkg/src/R_ldl.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 555 - (view) (download) (as text)

1 : bates 369 /* This is a modified version of the ldl.c file released by Timothy A. Davis */
2 :     /* in the LDL package and carrying the copyright shown below. The */
3 :     /* modifications are to replace scratch arrays passed as arguments by */
4 :     /* dynamically allocated arrays. */
5 :     /* Douglas Bates (Nov., 2004) */
6 :    
7 :     /* ========================================================================== */
8 :     /* === ldl.c: sparse LDL' factorization and solve package =================== */
9 :     /* ========================================================================== */
10 :    
11 :     /* LDL: a simple set of routines for sparse LDL' factorization. These routines
12 :     * are not terrifically fast (they do not use dense matrix kernels), but the
13 :     * code is very short. The purpose is to illustrate the algorithms in a very
14 :     * concise manner, primarily for educational purposes. Although the code is
15 :     * very concise, this package is slightly faster than the built-in sparse
16 :     * Cholesky factorization in MATLAB 6.5 (chol), when using the same input
17 :     * permutation.
18 :     *
19 :     * The routines compute the LDL' factorization of a real sparse symmetric
20 :     * matrix A (or PAP' if a permutation P is supplied), and solve upper
21 :     * and lower triangular systems with the resulting L and D factors. If A is
22 :     * positive definite then the factorization will be accurate. A can be
23 :     * indefinite (with negative values on the diagonal D), but in this case no
24 :     * guarantee of accuracy is provided, since no numeric pivoting is performed.
25 :     *
26 :     * The n-by-n sparse matrix A is in compressed-column form. The nonzero values
27 :     * in column j are stored in Ax [Ap [j] ... Ap [j+1]-1], with corresponding row
28 :     * indices in Ai [Ap [j] ... Ap [j+1]-1]. Ap [0] = 0 is required, and thus
29 :     * nz = Ap [n] is the number of nonzeros in A. Ap is an int array of size n+1.
30 :     * The int array Ai and the double array Ax are of size nz. This data structure
31 :     * is identical to the one used by MATLAB, except for the following
32 :     * generalizations. The row indices in each column of A need not be in any
33 :     * particular order, although they must be in the range 0 to n-1. Duplicate
34 :     * entries can be present; any duplicates are summed. That is, if row index i
35 :     * appears twice in a column j, then the value of A (i,j) is the sum of the two
36 :     * entries. The data structure used here for the input matrix A is more
37 :     * flexible than MATLAB's, which requires sorted columns with no duplicate
38 :     * entries.
39 :     *
40 :     * Only the diagonal and upper triangular part of A (or PAP' if a permutation
41 :     * P is provided) is accessed. The lower triangular parts of the matrix A or
42 :     * PAP' can be present, but they are ignored.
43 :     *
44 :     * The optional input permutation is provided as an array P of length n. If
45 :     * P [k] = j, the row and column j of A is the kth row and column of PAP'.
46 :     * If P is present then the factorization is LDL' = PAP' or L*D*L' = A(P,P) in
47 :     * 0-based MATLAB notation. If P is not present (a null pointer) then no
48 :     * permutation is performed, and the factorization is LDL' = A.
49 :     *
50 :     * The lower triangular matrix L is stored in the same compressed-column
51 :     * form (an int Lp array of size n+1, an int Li array of size Lp [n], and a
52 :     * double array Lx of the same size as Li). It has a unit diagonal, which is
53 :     * not stored. The row indices in each column of L are always returned in
54 :     * ascending order, with no duplicate entries. This format is compatible with
55 :     * MATLAB, except that it would be more convenient for MATLAB to include the
56 :     * unit diagonal of L. Doing so here would add additional complexity to the
57 :     * code, and is thus omitted in the interest of keeping this code short and
58 :     * readable.
59 :     *
60 :     * The elimination tree is held in the Parent [0..n-1] array. It is normally
61 :     * not required by the user, but it is required by R_ldl_numeric. The diagonal
62 :     * matrix D is held as an array D [0..n-1] of size n.
63 :     *
64 :     * --------------------
65 :     * C-callable routines:
66 :     * --------------------
67 :     *
68 : bates 555 * R_ldl_symbolic: Given the pattern of A, computes the Lp and Parent
69 :     * arrays required by R_ldl_numeric. Takes time proportional to the
70 :     * number of nonzeros in L. Computes the inverse Pinv of P if P is
71 :     * provided. Also returns Lnz, the count of nonzeros in each column
72 :     * of L below the diagonal (this is not required by R_ldl_numeric).
73 :     * R_ldl_numeric: Given the pattern and numerical values of A, the Lp
74 :     * array, the Parent array, and P and Pinv if applicable, computes
75 :     * the pattern and numerical values of L and D.
76 : bates 369 * R_ldl_lsolve: Solves Lx=b for a dense vector b.
77 :     * R_ldl_dsolve: Solves Dx=b for a dense vector b.
78 :     * R_ldl_ltsolve: Solves L'x=b for a dense vector b.
79 :     * R_ldl_perm: Computes x=Pb for a dense vector b.
80 :     * R-ldl_permt: Computes x=P'b for a dense vector b.
81 :     * R_ldl_valid_perm: checks the validity of a permutation vector
82 :     * R_ldl_valid_matrix: checks the validity of the sparse matrix A
83 :     *
84 :     * ----------------------------
85 :     * Limitations of this package:
86 :     * ----------------------------
87 :     *
88 : bates 555 * In the interest of keeping this code simple and readable,
89 :     * R_ldl_symbolic and R_ldl_numeric assume their inputs are valid.
90 :     * You can check your own inputs prior to calling these routines with
91 :     * the R_ldl_valid_perm and R_ldl_valid_matrix routines. Except for
92 :     * the two R_ldl_valid_* routines, no routine checks to see if the
93 :     * array arguments are present (non-NULL). Like all C routines, no
94 :     * routine can determine if the arrays are long enough and don't
95 :     * overlap.
96 : bates 369 *
97 : bates 555 * The R_ldl_numeric does check the numerical factorization, however.
98 :     * It returns n if the factorization is successful. If D (k,k) is
99 :     * zero, then k is returned, and L is only partially computed.
100 : bates 369 *
101 : bates 555 * No pivoting to control fill-in is performed, which is often
102 :     * critical for obtaining good performance. I recommend that you
103 :     * compute the permutation P using AMD or SYMAMD (approximate minimum
104 :     * degree ordering routines), or an appropriate graph-partitioning
105 :     * based ordering. See the ldldemo.m routine for an example in
106 :     * MATLAB, and the ldlmain.c stand-alone C program for examples of how
107 :     * to find P. Routines for manipulating compressed-column matrices
108 :     * are available in UMFPACK. AMD, SYMAMD, UMFPACK, and this LDL
109 :     * package are all available at
110 :     * http://www.cise.ufl.edu/research/sparse.
111 : bates 369 *
112 :     * -------------------------
113 :     * Possible simplifications:
114 :     * -------------------------
115 :     *
116 : bates 555 * These routines could be made even simpler with a few additional
117 :     * assumptions. If no input permutation were performed, the caller
118 :     * would have to permute the matrix first, but the computation of
119 :     * Pinv, and the use of P and Pinv could be removed. If only the
120 :     * diagonal and upper triangular part of A or PAP' are present, then
121 :     * the tests in the "if (i < k)" statement in R_ldl_symbolic and "if
122 :     * (i <= k)" in R_ldl_numeric, are always true, and could be removed
123 :     * (i can equal k in R_ldl_symbolic, but then the body of the if
124 :     * statement would correctly do no work since Flag [k] == k). If we
125 :     * could assume that no duplicate entries are present, then the
126 :     * statement Y [i] += Ax [p] could be replaced with Y [i] = Ax [p] in
127 :     * R_ldl_numeric.
128 : bates 369 *
129 :     * --------------------------
130 :     * Description of the method:
131 :     * --------------------------
132 :     *
133 : bates 555 * LDL computes the symbolic factorization by finding the pattern of L
134 :     * one row at a time. It does this based on the following theory.
135 :     * Consider a sparse system Lx=b, where L, x, and b, are all sparse,
136 :     * and where L comes from a Cholesky (or LDL') factorization. The
137 :     * elimination tree (etree) of L is defined as follows. The parent of
138 :     * node j is the smallest k > j such that L (k,j) is nonzero. Node j
139 :     * has no parent if column j of L is completely zero below the
140 :     * diagonal (j is a root of the etree in this case). The nonzero
141 :     * pattern of x is the union of the paths from each node i to the
142 :     * root, for each nonzero b (i). To compute the numerical solution to
143 :     * Lx=b, we can traverse the columns of L corresponding to nonzero
144 :     * values of x. This traversal does not need to be done in the order
145 :     * 0 to n-1. It can be done in any "topological" order, such that x
146 :     * (i) is computed before x (j) if i is a descendant of j in the
147 :     * elimination tree.
148 : bates 369 *
149 : bates 555 * The row-form of the LDL' factorization is shown in the MATLAB
150 :     * function ldlrow.m in this LDL package. Note that row k of L is
151 :     * found via a sparse triangular solve of L (1:k-1, 1:k-1) \ A (1:k-1,
152 :     * k), to use 1-based MATLAB notation. Thus, we can start with the
153 :     * nonzero pattern of the kth column of A (above the diagonal), follow
154 :     * the paths up to the root of the etree of the (k-1)-by-(k-1) leading
155 :     * submatrix of L, and obtain the pattern of the kth row of L. Note
156 :     * that we only need the leading (k-1)-by-(k-1) submatrix of L to do
157 :     * this. The elimination tree can be constructed as we go.
158 : bates 369 *
159 : bates 555 * The symbolic factorization does the same thing, except that it
160 :     * discards the pattern of L as it is computed. It simply counts the
161 :     * number of nonzeros in each column of L and then constructs the Lp
162 :     * index array when it's done. The symbolic factorization does not
163 :     * need to do this in topological order. Compare R_ldl_symbolic with
164 :     * the first part of R_ldl_numeric, and note that the while (len > 0)
165 :     * loop is not present in R_ldl_symbolic.
166 : bates 369 *
167 : bates 555 * LDL Version 1.0 (Dec. 31, 2003), Copyright (c) 2003 by Timothy A
168 :     * Davis, University of Florida. All Rights Reserved. Developed
169 :     * while on sabbatical at Stanford University and Lawrence Berkeley
170 :     * National Laboratory. Refer to the README file for the License.
171 :     * Available at http://www.cise.ufl.edu/research/sparse.
172 : bates 369 */
173 :    
174 :     #include "R_ldl.h"
175 :    
176 :     /* ========================================================================== */
177 : bates 555 /* === R_ldl_symbolic ======================================================= */
178 : bates 369 /* ========================================================================== */
179 :    
180 :     /**
181 : bates 555 * The input to this routine is a sparse matrix A, stored in column
182 :     * form, and an optional permutation P. The output is the elimination
183 :     * tree and the number of nonzeros in each column of L. Parent [i] =
184 :     * k if k is the parent of i in the tree. The Parent array is
185 :     * required by R_ldl_numeric. Lnz [k] gives the number of nonzeros in
186 :     * the kth column of L, excluding the diagonal.
187 : bates 369 *
188 : bates 555 * If P is NULL, then it is ignored. The factorization will be LDL' =
189 :     * A. Pinv is not computed. In this case, neither P nor Pinv are
190 :     * required by R_ldl_numeric.
191 : bates 369 *
192 :     * If P is not NULL, then it is assumed to be a valid permutation. If
193 : bates 555 * row and column j of A is the kth pivot, the P [k] = j. The
194 :     * factorization will be LDL' = PAP', or A (p,p) in MATLAB notation.
195 :     * The inverse permutation Pinv is computed, where Pinv [j] = k if P
196 :     * [k] = j. In this case, both P and Pinv are required as inputs to
197 :     * R_ldl_numeric.
198 : bates 369 *
199 : bates 555 * The floating-point operation count of the subsequent call to
200 :     * R_ldl_numeric is not returned, but could be computed after
201 :     * R_ldl_symbolic is done. It is the sum of (Lnz [k]) * (Lnz [k] + 2)
202 :     * for k = 0 to n-1.
203 : bates 369 *
204 :     * @param n A and L are n-by-n, where n >= 0
205 :     * @param Ap column pointers of size n+1
206 :     * @param Ai row indices of size nz=Ap[n]
207 :     * @param Lp column pointers of size n+1
208 :     * @param Parent elimination tree of size n
209 :     * @param P optional permutation vector of size n [use (int *) NULL for none]
210 :     * @param Pinv optional inverse permutation [not used if P is NULL]
211 :     */
212 :     void
213 :     R_ldl_symbolic(int n, const int Ap[], const int Ai[], int Lp[],
214 :     int Parent[], const int P[], int Pinv[])
215 :     {
216 :     int i, k, p, kk, p2;
217 :     int *Flag = Calloc(n, int);
218 :     int *Lnz = Calloc(n, int);
219 :     if (P)
220 :     {
221 :     /* If P is present then compute Pinv, the inverse of P */
222 :     for (k = 0 ; k < n ; k++)
223 :     {
224 :     Pinv [P [k]] = k ;
225 :     }
226 :     }
227 :     for (k = 0 ; k < n ; k++)
228 :     {
229 :     /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */
230 : bates 555 Parent [k] = -1 ; /* parent of k is not yet known */
231 :     Flag [k] = k ; /* mark node k as visited */
232 :     Lnz [k] = 0 ; /* count of nonzeros in column k of L */
233 :     kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */
234 : bates 369 p2 = Ap [kk+1] ;
235 :     for (p = Ap [kk] ; p < p2 ; p++)
236 :     {
237 :     /* A (i,k) is nonzero (original or permuted A) */
238 :     i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ;
239 :     if (i < k)
240 :     {
241 :     /* follow path from i to root of etree, stop at flagged node */
242 :     for ( ; Flag [i] != k ; i = Parent [i])
243 :     {
244 :     /* find parent of i if not yet determined */
245 :     if (Parent [i] == -1)
246 :     {
247 :     Parent [i] = k ;
248 :     }
249 : bates 555 Lnz [i]++ ; /* L (k,i) is nonzero */
250 :     Flag [i] = k ; /* mark i as visited */
251 : bates 369 }
252 :     }
253 :     }
254 :     }
255 :     /* construct Lp index array from Lnz column counts */
256 :     Lp [0] = 0 ;
257 :     for (k = 0 ; k < n ; k++)
258 :     {
259 :     Lp [k+1] = Lp [k] + Lnz [k] ;
260 :     }
261 :     Free(Flag); Free(Lnz);
262 :     }
263 :    
264 :     /**
265 : bates 555 * Given a sparse matrix A (the arguments n, Ap, Ai, and Ax) and its
266 :     * symbolic analysis (Lp and Parent, and optionally P and Pinv),
267 :     * compute the numeric LDL' factorization of A or PAP'. The outputs
268 :     * of this routine are arguments Li, Lx, and D.
269 : bates 369 *
270 :     * @param n A and L are n-by-n, where n >= 0
271 :     * @param Ap column pointer array of size n+1
272 :     * @param Ai row index array of size nz=Ap[n] (upper triangle only)
273 :     * @param Ax array of non-zero matrix elements of size nz=Ap[n]
274 :     * @param Lp column pointer array of size n+1
275 :     * @param Parent elimination tree of size n
276 :     * @param Li row index array of size lnz=Lp[n]
277 :     * @param Lx non-zero off-diagonal elements of L (size lnz=Lp[n])
278 :     * @param D vector of diagonal elements (size n)
279 :     * @param P optional permutation vector of size n [use (int *) NULL for none]
280 :     * @param Pinv optional inverse permutation [use (int *) NULL for none]
281 :     *
282 :     * @return n if successful, k if D (k,k) is zero
283 :     */
284 :     int
285 :     R_ldl_numeric(int n,
286 :     const int Ap[], const int Ai[], const double Ax[],
287 :     const int Lp[], const int Parent[],
288 :     int Li[], double Lx[], double D[],
289 :     const int P[], const int Pinv[])
290 :     {
291 :     double yi, l_ki ;
292 :     int i, k, p, kk, p2, len, top ;
293 :     int *Lnz = Calloc(n, int),
294 :     *Pattern = Calloc(n, int),
295 :     *Flag = Calloc(n, int);
296 :     double *Y = Calloc(n, double);
297 :    
298 :     for (k = 0 ; k < n ; k++)
299 :     {
300 :     /* compute nonzero Pattern of kth row of L, in topological order */
301 :     Y [k] = 0.0 ; /* Y (0:k) is now all zero */
302 : bates 555 top = n ; /* stack for pattern is empty */
303 :     Flag [k] = k ; /* mark node k as visited */
304 :     Lnz [k] = 0 ; /* count of nonzeros in column k of L */
305 :     kk = (P) ? (P [k]) : (k) ; /* kth original, or permuted, column */
306 : bates 369 p2 = Ap [kk+1] ;
307 :     for (p = Ap [kk] ; p < p2 ; p++)
308 :     {
309 : bates 555 i = (Pinv) ? (Pinv [Ai [p]]) : (Ai [p]) ; /* get A(i,k) */
310 : bates 369 if (i <= k)
311 :     {
312 : bates 555 Y [i] += Ax [p] ; /* scatter A(i,k) into Y (sum duplicates) */
313 : bates 369 /* follow path from i to root of etree, stop at flagged node */
314 :     for (len = 0 ; Flag [i] != k ; i = Parent [i])
315 :     {
316 : bates 555 Pattern [len++] = i ; /* L (k,i) is nonzero */
317 :     Flag [i] = k ; /* mark i as visited */
318 : bates 369 }
319 : bates 555 while (len > 0) /* push path on top of stack */
320 : bates 369 {
321 :     Pattern [--top] = Pattern [--len] ;
322 :     }
323 :     }
324 :     }
325 :     /* Pattern [top ... n-1] now contains nonzero pattern of L (:,k) */
326 :     /* compute numerical values kth row of L (a sparse triangular solve) */
327 : bates 555 D [k] = Y [k] ; /* get D (k,k) and clear Y (k) */
328 : bates 369 Y [k] = 0.0 ;
329 :     for ( ; top < n ; top++)
330 :     {
331 :     i = Pattern [top] ;
332 : bates 555 yi = Y [i] ; /* get and clear Y (i) */
333 : bates 369 Y [i] = 0.0 ;
334 :     p2 = Lp [i] + Lnz [i] ;
335 :     for (p = Lp [i] ; p < p2 ; p++)
336 :     {
337 :     Y [Li [p]] -= Lx [p] * yi ;
338 :     }
339 : bates 555 l_ki = yi / D [i] ; /* the nonzero entry L (k,i) */
340 : bates 369 D [k] -= l_ki * yi ;
341 : bates 555 Li [p] = k ; /* store L(k,k )in column form of L */
342 : bates 369 Lx [p] = l_ki ;
343 : bates 555 Lnz [i]++ ; /* increment count of nonzeros in col i */
344 : bates 369 }
345 :     if (D [k] == 0.0)
346 :     {
347 :     Free(Y); Free(Pattern); Free(Flag); Free(Lnz);
348 : bates 555 return (k) ; /* failure, D (k,k) is zero */
349 : bates 369 }
350 :     }
351 :     Free(Y); Free(Pattern); Free(Flag); Free(Lnz);
352 : bates 555 return (n) ; /* success, diagonal of D is all nonzero */
353 : bates 369 }
354 :    
355 :     /**
356 :     * solve Lx=b
357 :     *
358 :     * @param n L is n-by-n, where n >= 0
359 :     * @param X size n. right-hand-side on input, soln. on output
360 :     * @param Lp column pointer array of size n+1
361 :     * @param Li row index array of size lnz=Lp[n]
362 :     * @param Lx non-zero off-diagonal elements (size lnz=Lp[n])
363 :     */
364 :     void
365 :     R_ldl_lsolve (int n, double X[],
366 :     const int Lp[], const int Li[], const double Lx[])
367 :     {
368 :     int j, p, p2 ;
369 :     for (j = 0 ; j < n ; j++)
370 :     {
371 :     p2 = Lp [j+1] ;
372 :     for (p = Lp [j] ; p < p2 ; p++)
373 :     {
374 :     X [Li [p]] -= Lx [p] * X [j] ;
375 :     }
376 :     }
377 :     }
378 :    
379 :     /**
380 :     * solve Dx=b
381 :     *
382 :     * @param n L is n-by-n, where n >= 0
383 :     * @param X size n. right-hand-side on input, soln. on output
384 :     * @param D diagonal elements of size n
385 :     */
386 :     void
387 :     R_ldl_dsolve(int n, double X[], const double D[])
388 :     {
389 :     int j ;
390 :     for (j = 0 ; j < n ; j++)
391 :     {
392 :     X [j] /= D [j] ;
393 :     }
394 :     }
395 :    
396 :     /**
397 :     * solve L'x=b
398 :     *
399 :     * @param n L is n-by-n, where n >= 0
400 :     * @param X size n. right-hand-side on input, soln. on output
401 :     * @param Lp column pointer array of size n+1
402 :     * @param Li row index array of size lnz=Lp[n]
403 :     * @param Lx non-zero off-diagonal elements (size lnz=Lp[n])
404 :     */
405 :     void
406 :     R_ldl_ltsolve (int n, double X[],
407 :     const int Lp[], const int Li[], const double Lx[])
408 :     {
409 :     int j, p, p2 ;
410 :     for (j = n-1 ; j >= 0 ; j--)
411 :     {
412 :     p2 = Lp [j+1] ;
413 :     for (p = Lp [j] ; p < p2 ; p++)
414 :     {
415 :     X [j] -= Lx [p] * X [Li [p]] ;
416 :     }
417 :     }
418 :     }
419 :    
420 :     /**
421 :     * permute a vector, x=Pb
422 :     *
423 :     * @param n size of X, B, and P
424 :     * @param X output of size n
425 :     * @param B input of size n
426 :     * @param P input permutation array of size n
427 :     */
428 :     void
429 :     R_ldl_perm(int n, double X[], const double B[], const int P[])
430 :     {
431 :     int j ;
432 :     for (j = 0 ; j < n ; j++)
433 :     {
434 :     X [j] = B [P [j]] ;
435 :     }
436 :     }
437 :    
438 :     /**
439 :     * permute a vector, x=P'b
440 :     *
441 :     * @param n size of X, B, and P
442 :     * @param X output of size n
443 :     * @param B input of size n
444 :     * @param P input permutation array of size n
445 :     */
446 :     void
447 :     R_ldl_permt(int n, double X[], const double B[], const int P[])
448 :     {
449 :     int j ;
450 :     for (j = 0 ; j < n ; j++)
451 :     {
452 :     X [P [j]] = B [j] ;
453 :     }
454 :     }
455 :    
456 :     /**
457 :     * Check if a permutation vector is valid
458 :     *
459 :     * @param n size of permutation
460 :     * @param P input of size n, a permutation of 0:n-1
461 :     *
462 :     * @return 1 if valid, otherwise 0
463 :     */
464 :     int
465 :     R_ldl_valid_perm (int n, const int P[])
466 :     {
467 :    
468 :     int j, k ;
469 :     int *Flag = (int *) R_alloc(n, sizeof(int));
470 :    
471 :     if (n < 0 || !Flag)
472 :     {
473 : bates 555 return (0) ; /* n must be >= 0, and Flag must be present */
474 : bates 369 }
475 :     if (!P)
476 :     {
477 : bates 555 return (1) ; /* If NULL, P is assumed to be the identity perm. */
478 : bates 369 }
479 :     for (j = 0 ; j < n ; j++)
480 :     {
481 : bates 555 Flag [j] = 0 ; /* clear the Flag array */
482 : bates 369 }
483 :     for (k = 0 ; k < n ; k++)
484 :     {
485 :     j = P [k] ;
486 :     if (j < 0 || j >= n || Flag [j] != 0)
487 :     {
488 : bates 555 return (0) ; /* P is not valid */
489 : bates 369 }
490 :     Flag [j] = 1 ;
491 :     }
492 : bates 555 return (1) ; /* P is valid */
493 : bates 369 }
494 :    
495 :     /**
496 : bates 555 * This routine checks to see if a sparse matrix A is valid for input
497 :     * to R_ldl_symbolic and R_ldl_numeric. It returns 1 if the matrix is
498 :     * valid, 0 otherwise. A is in sparse column form. The numerical
499 :     * values in column j are stored in Ax [Ap [j] ... Ap [j+1]-1], with
500 :     * row indices in Ai [Ap [j] ... Ap [j+1]-1]. The Ax array is not
501 :     * checked.
502 : bates 369 *
503 :     * @param n A is n by n (n >= 0)
504 :     * @param Ap column pointer array of size n+1
505 :     * @param Ai row index array of size nz=Ap[n]
506 :     *
507 :     * @return 1 if valid sparse matrix, otherwise 0
508 :     */
509 :     int
510 :     R_ldl_valid_matrix (int n, const int Ap[], const int Ai[])
511 :     {
512 :     int j, p ;
513 :     if (n < 0 || !Ap || !Ai || Ap [0] != 0)
514 :     {
515 : bates 555 return (0) ; /* n must be >= 0, and Ap and Ai must be present */
516 : bates 369 }
517 :     for (j = 0 ; j < n ; j++)
518 :     {
519 :     if (Ap [j] > Ap [j+1])
520 :     {
521 : bates 555 return (0) ; /* Ap must be monotonically nondecreasing */
522 : bates 369 }
523 :     }
524 :     for (p = 0 ; p < Ap [n] ; p++)
525 :     {
526 :     if (Ai [p] < 0 || Ai [p] >= n)
527 :     {
528 : bates 555 return (0) ; /* row indices must be in the range 0 to n-1 */
529 : bates 369 }
530 :     }
531 : bates 555 return (1) ; /* matrix is valid */
532 : bates 369 }

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge