SCM

SCM Repository

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

Annotation of /pkg/src/ldl.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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

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