1 | /* glpmat.c */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
9 | * E-mail: <mao@gnu.org>. |
---|
10 | * |
---|
11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
12 | * under the terms of the GNU General Public License as published by |
---|
13 | * the Free Software Foundation, either version 3 of the License, or |
---|
14 | * (at your option) any later version. |
---|
15 | * |
---|
16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
19 | * License for more details. |
---|
20 | * |
---|
21 | * You should have received a copy of the GNU General Public License |
---|
22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
23 | ***********************************************************************/ |
---|
24 | |
---|
25 | #include "glpenv.h" |
---|
26 | #include "glpmat.h" |
---|
27 | #include "glpqmd.h" |
---|
28 | #include "amd/amd.h" |
---|
29 | #include "colamd/colamd.h" |
---|
30 | |
---|
31 | /*---------------------------------------------------------------------- |
---|
32 | -- check_fvs - check sparse vector in full-vector storage format. |
---|
33 | -- |
---|
34 | -- SYNOPSIS |
---|
35 | -- |
---|
36 | -- #include "glpmat.h" |
---|
37 | -- int check_fvs(int n, int nnz, int ind[], double vec[]); |
---|
38 | -- |
---|
39 | -- DESCRIPTION |
---|
40 | -- |
---|
41 | -- The routine check_fvs checks if a given vector of dimension n in |
---|
42 | -- full-vector storage format has correct representation. |
---|
43 | -- |
---|
44 | -- RETURNS |
---|
45 | -- |
---|
46 | -- The routine returns one of the following codes: |
---|
47 | -- |
---|
48 | -- 0 - the vector is correct; |
---|
49 | -- 1 - the number of elements (n) is negative; |
---|
50 | -- 2 - the number of non-zero elements (nnz) is negative; |
---|
51 | -- 3 - some element index is out of range; |
---|
52 | -- 4 - some element index is duplicate; |
---|
53 | -- 5 - some non-zero element is out of pattern. */ |
---|
54 | |
---|
55 | int check_fvs(int n, int nnz, int ind[], double vec[]) |
---|
56 | { int i, t, ret, *flag = NULL; |
---|
57 | /* check the number of elements */ |
---|
58 | if (n < 0) |
---|
59 | { ret = 1; |
---|
60 | goto done; |
---|
61 | } |
---|
62 | /* check the number of non-zero elements */ |
---|
63 | if (nnz < 0) |
---|
64 | { ret = 2; |
---|
65 | goto done; |
---|
66 | } |
---|
67 | /* check vector indices */ |
---|
68 | flag = xcalloc(1+n, sizeof(int)); |
---|
69 | for (i = 1; i <= n; i++) flag[i] = 0; |
---|
70 | for (t = 1; t <= nnz; t++) |
---|
71 | { i = ind[t]; |
---|
72 | if (!(1 <= i && i <= n)) |
---|
73 | { ret = 3; |
---|
74 | goto done; |
---|
75 | } |
---|
76 | if (flag[i]) |
---|
77 | { ret = 4; |
---|
78 | goto done; |
---|
79 | } |
---|
80 | flag[i] = 1; |
---|
81 | } |
---|
82 | /* check vector elements */ |
---|
83 | for (i = 1; i <= n; i++) |
---|
84 | { if (!flag[i] && vec[i] != 0.0) |
---|
85 | { ret = 5; |
---|
86 | goto done; |
---|
87 | } |
---|
88 | } |
---|
89 | /* the vector is ok */ |
---|
90 | ret = 0; |
---|
91 | done: if (flag != NULL) xfree(flag); |
---|
92 | return ret; |
---|
93 | } |
---|
94 | |
---|
95 | /*---------------------------------------------------------------------- |
---|
96 | -- check_pattern - check pattern of sparse matrix. |
---|
97 | -- |
---|
98 | -- SYNOPSIS |
---|
99 | -- |
---|
100 | -- #include "glpmat.h" |
---|
101 | -- int check_pattern(int m, int n, int A_ptr[], int A_ind[]); |
---|
102 | -- |
---|
103 | -- DESCRIPTION |
---|
104 | -- |
---|
105 | -- The routine check_pattern checks the pattern of a given mxn matrix |
---|
106 | -- in storage-by-rows format. |
---|
107 | -- |
---|
108 | -- RETURNS |
---|
109 | -- |
---|
110 | -- The routine returns one of the following codes: |
---|
111 | -- |
---|
112 | -- 0 - the pattern is correct; |
---|
113 | -- 1 - the number of rows (m) is negative; |
---|
114 | -- 2 - the number of columns (n) is negative; |
---|
115 | -- 3 - A_ptr[1] is not 1; |
---|
116 | -- 4 - some column index is out of range; |
---|
117 | -- 5 - some column indices are duplicate. */ |
---|
118 | |
---|
119 | int check_pattern(int m, int n, int A_ptr[], int A_ind[]) |
---|
120 | { int i, j, ptr, ret, *flag = NULL; |
---|
121 | /* check the number of rows */ |
---|
122 | if (m < 0) |
---|
123 | { ret = 1; |
---|
124 | goto done; |
---|
125 | } |
---|
126 | /* check the number of columns */ |
---|
127 | if (n < 0) |
---|
128 | { ret = 2; |
---|
129 | goto done; |
---|
130 | } |
---|
131 | /* check location A_ptr[1] */ |
---|
132 | if (A_ptr[1] != 1) |
---|
133 | { ret = 3; |
---|
134 | goto done; |
---|
135 | } |
---|
136 | /* check row patterns */ |
---|
137 | flag = xcalloc(1+n, sizeof(int)); |
---|
138 | for (j = 1; j <= n; j++) flag[j] = 0; |
---|
139 | for (i = 1; i <= m; i++) |
---|
140 | { /* check pattern of row i */ |
---|
141 | for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++) |
---|
142 | { j = A_ind[ptr]; |
---|
143 | /* check column index */ |
---|
144 | if (!(1 <= j && j <= n)) |
---|
145 | { ret = 4; |
---|
146 | goto done; |
---|
147 | } |
---|
148 | /* check for duplication */ |
---|
149 | if (flag[j]) |
---|
150 | { ret = 5; |
---|
151 | goto done; |
---|
152 | } |
---|
153 | flag[j] = 1; |
---|
154 | } |
---|
155 | /* clear flags */ |
---|
156 | for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++) |
---|
157 | { j = A_ind[ptr]; |
---|
158 | flag[j] = 0; |
---|
159 | } |
---|
160 | } |
---|
161 | /* the pattern is ok */ |
---|
162 | ret = 0; |
---|
163 | done: if (flag != NULL) xfree(flag); |
---|
164 | return ret; |
---|
165 | } |
---|
166 | |
---|
167 | /*---------------------------------------------------------------------- |
---|
168 | -- transpose - transpose sparse matrix. |
---|
169 | -- |
---|
170 | -- *Synopsis* |
---|
171 | -- |
---|
172 | -- #include "glpmat.h" |
---|
173 | -- void transpose(int m, int n, int A_ptr[], int A_ind[], |
---|
174 | -- double A_val[], int AT_ptr[], int AT_ind[], double AT_val[]); |
---|
175 | -- |
---|
176 | -- *Description* |
---|
177 | -- |
---|
178 | -- For a given mxn sparse matrix A the routine transpose builds a nxm |
---|
179 | -- sparse matrix A' which is a matrix transposed to A. |
---|
180 | -- |
---|
181 | -- The arrays A_ptr, A_ind, and A_val specify a given mxn matrix A to |
---|
182 | -- be transposed in storage-by-rows format. The parameter A_val can be |
---|
183 | -- NULL, in which case numeric values are not copied. The arrays A_ptr, |
---|
184 | -- A_ind, and A_val are not changed on exit. |
---|
185 | -- |
---|
186 | -- On entry the arrays AT_ptr, AT_ind, and AT_val must be allocated, |
---|
187 | -- but their content is ignored. On exit the routine stores a resultant |
---|
188 | -- nxm matrix A' in these arrays in storage-by-rows format. Note that |
---|
189 | -- if the parameter A_val is NULL, the array AT_val is not used. |
---|
190 | -- |
---|
191 | -- The routine transpose has a side effect that elements in rows of the |
---|
192 | -- resultant matrix A' follow in ascending their column indices. */ |
---|
193 | |
---|
194 | void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[], |
---|
195 | int AT_ptr[], int AT_ind[], double AT_val[]) |
---|
196 | { int i, j, t, beg, end, pos, len; |
---|
197 | /* determine row lengths of resultant matrix */ |
---|
198 | for (j = 1; j <= n; j++) AT_ptr[j] = 0; |
---|
199 | for (i = 1; i <= m; i++) |
---|
200 | { beg = A_ptr[i], end = A_ptr[i+1]; |
---|
201 | for (t = beg; t < end; t++) AT_ptr[A_ind[t]]++; |
---|
202 | } |
---|
203 | /* set up row pointers of resultant matrix */ |
---|
204 | pos = 1; |
---|
205 | for (j = 1; j <= n; j++) |
---|
206 | len = AT_ptr[j], pos += len, AT_ptr[j] = pos; |
---|
207 | AT_ptr[n+1] = pos; |
---|
208 | /* build resultant matrix */ |
---|
209 | for (i = m; i >= 1; i--) |
---|
210 | { beg = A_ptr[i], end = A_ptr[i+1]; |
---|
211 | for (t = beg; t < end; t++) |
---|
212 | { pos = --AT_ptr[A_ind[t]]; |
---|
213 | AT_ind[pos] = i; |
---|
214 | if (A_val != NULL) AT_val[pos] = A_val[t]; |
---|
215 | } |
---|
216 | } |
---|
217 | return; |
---|
218 | } |
---|
219 | |
---|
220 | /*---------------------------------------------------------------------- |
---|
221 | -- adat_symbolic - compute S = P*A*D*A'*P' (symbolic phase). |
---|
222 | -- |
---|
223 | -- *Synopsis* |
---|
224 | -- |
---|
225 | -- #include "glpmat.h" |
---|
226 | -- int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], |
---|
227 | -- int A_ind[], int S_ptr[]); |
---|
228 | -- |
---|
229 | -- *Description* |
---|
230 | -- |
---|
231 | -- The routine adat_symbolic implements the symbolic phase to compute |
---|
232 | -- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix, |
---|
233 | -- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix |
---|
234 | -- transposed to A, P' is an inverse of P. |
---|
235 | -- |
---|
236 | -- The parameter m is the number of rows in A and the order of P. |
---|
237 | -- |
---|
238 | -- The parameter n is the number of columns in A and the order of D. |
---|
239 | -- |
---|
240 | -- The array P_per specifies permutation matrix P. It is not changed on |
---|
241 | -- exit. |
---|
242 | -- |
---|
243 | -- The arrays A_ptr and A_ind specify the pattern of matrix A. They are |
---|
244 | -- not changed on exit. |
---|
245 | -- |
---|
246 | -- On exit the routine stores the pattern of upper triangular part of |
---|
247 | -- matrix S without diagonal elements in the arrays S_ptr and S_ind in |
---|
248 | -- storage-by-rows format. The array S_ptr should be allocated on entry, |
---|
249 | -- however, its content is ignored. The array S_ind is allocated by the |
---|
250 | -- routine itself which returns a pointer to it. |
---|
251 | -- |
---|
252 | -- *Returns* |
---|
253 | -- |
---|
254 | -- The routine returns a pointer to the array S_ind. */ |
---|
255 | |
---|
256 | int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[], |
---|
257 | int S_ptr[]) |
---|
258 | { int i, j, t, ii, jj, tt, k, size, len; |
---|
259 | int *S_ind, *AT_ptr, *AT_ind, *ind, *map, *temp; |
---|
260 | /* build the pattern of A', which is a matrix transposed to A, to |
---|
261 | efficiently access A in column-wise manner */ |
---|
262 | AT_ptr = xcalloc(1+n+1, sizeof(int)); |
---|
263 | AT_ind = xcalloc(A_ptr[m+1], sizeof(int)); |
---|
264 | transpose(m, n, A_ptr, A_ind, NULL, AT_ptr, AT_ind, NULL); |
---|
265 | /* allocate the array S_ind */ |
---|
266 | size = A_ptr[m+1] - 1; |
---|
267 | if (size < m) size = m; |
---|
268 | S_ind = xcalloc(1+size, sizeof(int)); |
---|
269 | /* allocate and initialize working arrays */ |
---|
270 | ind = xcalloc(1+m, sizeof(int)); |
---|
271 | map = xcalloc(1+m, sizeof(int)); |
---|
272 | for (jj = 1; jj <= m; jj++) map[jj] = 0; |
---|
273 | /* compute pattern of S; note that symbolically S = B*B', where |
---|
274 | B = P*A, B' is matrix transposed to B */ |
---|
275 | S_ptr[1] = 1; |
---|
276 | for (ii = 1; ii <= m; ii++) |
---|
277 | { /* compute pattern of ii-th row of S */ |
---|
278 | len = 0; |
---|
279 | i = P_per[ii]; /* i-th row of A = ii-th row of B */ |
---|
280 | for (t = A_ptr[i]; t < A_ptr[i+1]; t++) |
---|
281 | { k = A_ind[t]; |
---|
282 | /* walk through k-th column of A */ |
---|
283 | for (tt = AT_ptr[k]; tt < AT_ptr[k+1]; tt++) |
---|
284 | { j = AT_ind[tt]; |
---|
285 | jj = P_per[m+j]; /* j-th row of A = jj-th row of B */ |
---|
286 | /* a[i,k] != 0 and a[j,k] != 0 ergo s[ii,jj] != 0 */ |
---|
287 | if (ii < jj && !map[jj]) ind[++len] = jj, map[jj] = 1; |
---|
288 | } |
---|
289 | } |
---|
290 | /* now (ind) is pattern of ii-th row of S */ |
---|
291 | S_ptr[ii+1] = S_ptr[ii] + len; |
---|
292 | /* at least (S_ptr[ii+1] - 1) locations should be available in |
---|
293 | the array S_ind */ |
---|
294 | if (S_ptr[ii+1] - 1 > size) |
---|
295 | { temp = S_ind; |
---|
296 | size += size; |
---|
297 | S_ind = xcalloc(1+size, sizeof(int)); |
---|
298 | memcpy(&S_ind[1], &temp[1], (S_ptr[ii] - 1) * sizeof(int)); |
---|
299 | xfree(temp); |
---|
300 | } |
---|
301 | xassert(S_ptr[ii+1] - 1 <= size); |
---|
302 | /* (ii-th row of S) := (ind) */ |
---|
303 | memcpy(&S_ind[S_ptr[ii]], &ind[1], len * sizeof(int)); |
---|
304 | /* clear the row pattern map */ |
---|
305 | for (t = 1; t <= len; t++) map[ind[t]] = 0; |
---|
306 | } |
---|
307 | /* free working arrays */ |
---|
308 | xfree(AT_ptr); |
---|
309 | xfree(AT_ind); |
---|
310 | xfree(ind); |
---|
311 | xfree(map); |
---|
312 | /* reallocate the array S_ind to free unused locations */ |
---|
313 | temp = S_ind; |
---|
314 | size = S_ptr[m+1] - 1; |
---|
315 | S_ind = xcalloc(1+size, sizeof(int)); |
---|
316 | memcpy(&S_ind[1], &temp[1], size * sizeof(int)); |
---|
317 | xfree(temp); |
---|
318 | return S_ind; |
---|
319 | } |
---|
320 | |
---|
321 | /*---------------------------------------------------------------------- |
---|
322 | -- adat_numeric - compute S = P*A*D*A'*P' (numeric phase). |
---|
323 | -- |
---|
324 | -- *Synopsis* |
---|
325 | -- |
---|
326 | -- #include "glpmat.h" |
---|
327 | -- void adat_numeric(int m, int n, int P_per[], |
---|
328 | -- int A_ptr[], int A_ind[], double A_val[], double D_diag[], |
---|
329 | -- int S_ptr[], int S_ind[], double S_val[], double S_diag[]); |
---|
330 | -- |
---|
331 | -- *Description* |
---|
332 | -- |
---|
333 | -- The routine adat_numeric implements the numeric phase to compute |
---|
334 | -- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix, |
---|
335 | -- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix |
---|
336 | -- transposed to A, P' is an inverse of P. |
---|
337 | -- |
---|
338 | -- The parameter m is the number of rows in A and the order of P. |
---|
339 | -- |
---|
340 | -- The parameter n is the number of columns in A and the order of D. |
---|
341 | -- |
---|
342 | -- The matrix P is specified in the array P_per, which is not changed |
---|
343 | -- on exit. |
---|
344 | -- |
---|
345 | -- The matrix A is specified in the arrays A_ptr, A_ind, and A_val in |
---|
346 | -- storage-by-rows format. These arrays are not changed on exit. |
---|
347 | -- |
---|
348 | -- Diagonal elements of the matrix D are specified in the array D_diag, |
---|
349 | -- where D_diag[0] is not used, D_diag[i] = d[i,i] for i = 1, ..., n. |
---|
350 | -- The array D_diag is not changed on exit. |
---|
351 | -- |
---|
352 | -- The pattern of the upper triangular part of the matrix S without |
---|
353 | -- diagonal elements (previously computed by the routine adat_symbolic) |
---|
354 | -- is specified in the arrays S_ptr and S_ind, which are not changed on |
---|
355 | -- exit. Numeric values of non-diagonal elements of S are stored in |
---|
356 | -- corresponding locations of the array S_val, and values of diagonal |
---|
357 | -- elements of S are stored in locations S_diag[1], ..., S_diag[n]. */ |
---|
358 | |
---|
359 | void adat_numeric(int m, int n, int P_per[], |
---|
360 | int A_ptr[], int A_ind[], double A_val[], double D_diag[], |
---|
361 | int S_ptr[], int S_ind[], double S_val[], double S_diag[]) |
---|
362 | { int i, j, t, ii, jj, tt, beg, end, beg1, end1, k; |
---|
363 | double sum, *work; |
---|
364 | work = xcalloc(1+n, sizeof(double)); |
---|
365 | for (j = 1; j <= n; j++) work[j] = 0.0; |
---|
366 | /* compute S = B*D*B', where B = P*A, B' is a matrix transposed |
---|
367 | to B */ |
---|
368 | for (ii = 1; ii <= m; ii++) |
---|
369 | { i = P_per[ii]; /* i-th row of A = ii-th row of B */ |
---|
370 | /* (work) := (i-th row of A) */ |
---|
371 | beg = A_ptr[i], end = A_ptr[i+1]; |
---|
372 | for (t = beg; t < end; t++) |
---|
373 | work[A_ind[t]] = A_val[t]; |
---|
374 | /* compute ii-th row of S */ |
---|
375 | beg = S_ptr[ii], end = S_ptr[ii+1]; |
---|
376 | for (t = beg; t < end; t++) |
---|
377 | { jj = S_ind[t]; |
---|
378 | j = P_per[jj]; /* j-th row of A = jj-th row of B */ |
---|
379 | /* s[ii,jj] := sum a[i,k] * d[k,k] * a[j,k] */ |
---|
380 | sum = 0.0; |
---|
381 | beg1 = A_ptr[j], end1 = A_ptr[j+1]; |
---|
382 | for (tt = beg1; tt < end1; tt++) |
---|
383 | { k = A_ind[tt]; |
---|
384 | sum += work[k] * D_diag[k] * A_val[tt]; |
---|
385 | } |
---|
386 | S_val[t] = sum; |
---|
387 | } |
---|
388 | /* s[ii,ii] := sum a[i,k] * d[k,k] * a[i,k] */ |
---|
389 | sum = 0.0; |
---|
390 | beg = A_ptr[i], end = A_ptr[i+1]; |
---|
391 | for (t = beg; t < end; t++) |
---|
392 | { k = A_ind[t]; |
---|
393 | sum += A_val[t] * D_diag[k] * A_val[t]; |
---|
394 | work[k] = 0.0; |
---|
395 | } |
---|
396 | S_diag[ii] = sum; |
---|
397 | } |
---|
398 | xfree(work); |
---|
399 | return; |
---|
400 | } |
---|
401 | |
---|
402 | /*---------------------------------------------------------------------- |
---|
403 | -- min_degree - minimum degree ordering. |
---|
404 | -- |
---|
405 | -- *Synopsis* |
---|
406 | -- |
---|
407 | -- #include "glpmat.h" |
---|
408 | -- void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]); |
---|
409 | -- |
---|
410 | -- *Description* |
---|
411 | -- |
---|
412 | -- The routine min_degree uses the minimum degree ordering algorithm |
---|
413 | -- to find a permutation matrix P for a given sparse symmetric positive |
---|
414 | -- matrix A which minimizes the number of non-zeros in upper triangular |
---|
415 | -- factor U for Cholesky factorization P*A*P' = U'*U. |
---|
416 | -- |
---|
417 | -- The parameter n is the order of matrices A and P. |
---|
418 | -- |
---|
419 | -- The pattern of the given matrix A is specified on entry in the arrays |
---|
420 | -- A_ptr and A_ind in storage-by-rows format. Only the upper triangular |
---|
421 | -- part without diagonal elements (which all are assumed to be non-zero) |
---|
422 | -- should be specified as if A were upper triangular. The arrays A_ptr |
---|
423 | -- and A_ind are not changed on exit. |
---|
424 | -- |
---|
425 | -- The permutation matrix P is stored by the routine in the array P_per |
---|
426 | -- on exit. |
---|
427 | -- |
---|
428 | -- *Algorithm* |
---|
429 | -- |
---|
430 | -- The routine min_degree is based on some subroutines from the package |
---|
431 | -- SPARSPAK (see comments in the module glpqmd). */ |
---|
432 | |
---|
433 | void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]) |
---|
434 | { int i, j, ne, t, pos, len; |
---|
435 | int *xadj, *adjncy, *deg, *marker, *rchset, *nbrhd, *qsize, |
---|
436 | *qlink, nofsub; |
---|
437 | /* determine number of non-zeros in complete pattern */ |
---|
438 | ne = A_ptr[n+1] - 1; |
---|
439 | ne += ne; |
---|
440 | /* allocate working arrays */ |
---|
441 | xadj = xcalloc(1+n+1, sizeof(int)); |
---|
442 | adjncy = xcalloc(1+ne, sizeof(int)); |
---|
443 | deg = xcalloc(1+n, sizeof(int)); |
---|
444 | marker = xcalloc(1+n, sizeof(int)); |
---|
445 | rchset = xcalloc(1+n, sizeof(int)); |
---|
446 | nbrhd = xcalloc(1+n, sizeof(int)); |
---|
447 | qsize = xcalloc(1+n, sizeof(int)); |
---|
448 | qlink = xcalloc(1+n, sizeof(int)); |
---|
449 | /* determine row lengths in complete pattern */ |
---|
450 | for (i = 1; i <= n; i++) xadj[i] = 0; |
---|
451 | for (i = 1; i <= n; i++) |
---|
452 | { for (t = A_ptr[i]; t < A_ptr[i+1]; t++) |
---|
453 | { j = A_ind[t]; |
---|
454 | xassert(i < j && j <= n); |
---|
455 | xadj[i]++, xadj[j]++; |
---|
456 | } |
---|
457 | } |
---|
458 | /* set up row pointers for complete pattern */ |
---|
459 | pos = 1; |
---|
460 | for (i = 1; i <= n; i++) |
---|
461 | len = xadj[i], pos += len, xadj[i] = pos; |
---|
462 | xadj[n+1] = pos; |
---|
463 | xassert(pos - 1 == ne); |
---|
464 | /* construct complete pattern */ |
---|
465 | for (i = 1; i <= n; i++) |
---|
466 | { for (t = A_ptr[i]; t < A_ptr[i+1]; t++) |
---|
467 | { j = A_ind[t]; |
---|
468 | adjncy[--xadj[i]] = j, adjncy[--xadj[j]] = i; |
---|
469 | } |
---|
470 | } |
---|
471 | /* call the main minimimum degree ordering routine */ |
---|
472 | genqmd(&n, xadj, adjncy, P_per, P_per + n, deg, marker, rchset, |
---|
473 | nbrhd, qsize, qlink, &nofsub); |
---|
474 | /* make sure that permutation matrix P is correct */ |
---|
475 | for (i = 1; i <= n; i++) |
---|
476 | { j = P_per[i]; |
---|
477 | xassert(1 <= j && j <= n); |
---|
478 | xassert(P_per[n+j] == i); |
---|
479 | } |
---|
480 | /* free working arrays */ |
---|
481 | xfree(xadj); |
---|
482 | xfree(adjncy); |
---|
483 | xfree(deg); |
---|
484 | xfree(marker); |
---|
485 | xfree(rchset); |
---|
486 | xfree(nbrhd); |
---|
487 | xfree(qsize); |
---|
488 | xfree(qlink); |
---|
489 | return; |
---|
490 | } |
---|
491 | |
---|
492 | /**********************************************************************/ |
---|
493 | |
---|
494 | void amd_order1(int n, int A_ptr[], int A_ind[], int P_per[]) |
---|
495 | { /* approximate minimum degree ordering (AMD) */ |
---|
496 | int k, ret; |
---|
497 | double Control[AMD_CONTROL], Info[AMD_INFO]; |
---|
498 | /* get the default parameters */ |
---|
499 | amd_defaults(Control); |
---|
500 | #if 0 |
---|
501 | /* and print them */ |
---|
502 | amd_control(Control); |
---|
503 | #endif |
---|
504 | /* make all indices 0-based */ |
---|
505 | for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--; |
---|
506 | for (k = 1; k <= n+1; k++) A_ptr[k]--; |
---|
507 | /* call the ordering routine */ |
---|
508 | ret = amd_order(n, &A_ptr[1], &A_ind[1], &P_per[1], Control, Info) |
---|
509 | ; |
---|
510 | #if 0 |
---|
511 | amd_info(Info); |
---|
512 | #endif |
---|
513 | xassert(ret == AMD_OK || ret == AMD_OK_BUT_JUMBLED); |
---|
514 | /* retsore 1-based indices */ |
---|
515 | for (k = 1; k <= n+1; k++) A_ptr[k]++; |
---|
516 | for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++; |
---|
517 | /* patch up permutation matrix */ |
---|
518 | memset(&P_per[n+1], 0, n * sizeof(int)); |
---|
519 | for (k = 1; k <= n; k++) |
---|
520 | { P_per[k]++; |
---|
521 | xassert(1 <= P_per[k] && P_per[k] <= n); |
---|
522 | xassert(P_per[n+P_per[k]] == 0); |
---|
523 | P_per[n+P_per[k]] = k; |
---|
524 | } |
---|
525 | return; |
---|
526 | } |
---|
527 | |
---|
528 | /**********************************************************************/ |
---|
529 | |
---|
530 | static void *allocate(size_t n, size_t size) |
---|
531 | { void *ptr; |
---|
532 | ptr = xcalloc(n, size); |
---|
533 | memset(ptr, 0, n * size); |
---|
534 | return ptr; |
---|
535 | } |
---|
536 | |
---|
537 | static void release(void *ptr) |
---|
538 | { xfree(ptr); |
---|
539 | return; |
---|
540 | } |
---|
541 | |
---|
542 | void symamd_ord(int n, int A_ptr[], int A_ind[], int P_per[]) |
---|
543 | { /* approximate minimum degree ordering (SYMAMD) */ |
---|
544 | int k, ok; |
---|
545 | int stats[COLAMD_STATS]; |
---|
546 | /* make all indices 0-based */ |
---|
547 | for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--; |
---|
548 | for (k = 1; k <= n+1; k++) A_ptr[k]--; |
---|
549 | /* call the ordering routine */ |
---|
550 | ok = symamd(n, &A_ind[1], &A_ptr[1], &P_per[1], NULL, stats, |
---|
551 | allocate, release); |
---|
552 | #if 0 |
---|
553 | symamd_report(stats); |
---|
554 | #endif |
---|
555 | xassert(ok); |
---|
556 | /* restore 1-based indices */ |
---|
557 | for (k = 1; k <= n+1; k++) A_ptr[k]++; |
---|
558 | for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++; |
---|
559 | /* patch up permutation matrix */ |
---|
560 | memset(&P_per[n+1], 0, n * sizeof(int)); |
---|
561 | for (k = 1; k <= n; k++) |
---|
562 | { P_per[k]++; |
---|
563 | xassert(1 <= P_per[k] && P_per[k] <= n); |
---|
564 | xassert(P_per[n+P_per[k]] == 0); |
---|
565 | P_per[n+P_per[k]] = k; |
---|
566 | } |
---|
567 | return; |
---|
568 | } |
---|
569 | |
---|
570 | /*---------------------------------------------------------------------- |
---|
571 | -- chol_symbolic - compute Cholesky factorization (symbolic phase). |
---|
572 | -- |
---|
573 | -- *Synopsis* |
---|
574 | -- |
---|
575 | -- #include "glpmat.h" |
---|
576 | -- int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]); |
---|
577 | -- |
---|
578 | -- *Description* |
---|
579 | -- |
---|
580 | -- The routine chol_symbolic implements the symbolic phase of Cholesky |
---|
581 | -- factorization A = U'*U, where A is a given sparse symmetric positive |
---|
582 | -- definite matrix, U is a resultant upper triangular factor, U' is a |
---|
583 | -- matrix transposed to U. |
---|
584 | -- |
---|
585 | -- The parameter n is the order of matrices A and U. |
---|
586 | -- |
---|
587 | -- The pattern of the given matrix A is specified on entry in the arrays |
---|
588 | -- A_ptr and A_ind in storage-by-rows format. Only the upper triangular |
---|
589 | -- part without diagonal elements (which all are assumed to be non-zero) |
---|
590 | -- should be specified as if A were upper triangular. The arrays A_ptr |
---|
591 | -- and A_ind are not changed on exit. |
---|
592 | -- |
---|
593 | -- The pattern of the matrix U without diagonal elements (which all are |
---|
594 | -- assumed to be non-zero) is stored on exit from the routine in the |
---|
595 | -- arrays U_ptr and U_ind in storage-by-rows format. The array U_ptr |
---|
596 | -- should be allocated on entry, however, its content is ignored. The |
---|
597 | -- array U_ind is allocated by the routine which returns a pointer to it |
---|
598 | -- on exit. |
---|
599 | -- |
---|
600 | -- *Returns* |
---|
601 | -- |
---|
602 | -- The routine returns a pointer to the array U_ind. |
---|
603 | -- |
---|
604 | -- *Method* |
---|
605 | -- |
---|
606 | -- The routine chol_symbolic computes the pattern of the matrix U in a |
---|
607 | -- row-wise manner. No pivoting is used. |
---|
608 | -- |
---|
609 | -- It is known that to compute the pattern of row k of the matrix U we |
---|
610 | -- need to merge the pattern of row k of the matrix A and the patterns |
---|
611 | -- of each row i of U, where u[i,k] is non-zero (these rows are already |
---|
612 | -- computed and placed above row k). |
---|
613 | -- |
---|
614 | -- However, to reduce the number of rows to be merged the routine uses |
---|
615 | -- an advanced algorithm proposed in: |
---|
616 | -- |
---|
617 | -- D.J.Rose, R.E.Tarjan, and G.S.Lueker. Algorithmic aspects of vertex |
---|
618 | -- elimination on graphs. SIAM J. Comput. 5, 1976, 266-83. |
---|
619 | -- |
---|
620 | -- The authors of the cited paper show that we have the same result if |
---|
621 | -- we merge row k of the matrix A and such rows of the matrix U (among |
---|
622 | -- rows 1, ..., k-1) whose leftmost non-diagonal non-zero element is |
---|
623 | -- placed in k-th column. This feature signficantly reduces the number |
---|
624 | -- of rows to be merged, especially on the final steps, where rows of |
---|
625 | -- the matrix U become quite dense. |
---|
626 | -- |
---|
627 | -- To determine rows, which should be merged on k-th step, for a fixed |
---|
628 | -- time the routine uses linked lists of row numbers of the matrix U. |
---|
629 | -- Location head[k] contains the number of a first row, whose leftmost |
---|
630 | -- non-diagonal non-zero element is placed in column k, and location |
---|
631 | -- next[i] contains the number of a next row with the same property as |
---|
632 | -- row i. */ |
---|
633 | |
---|
634 | int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]) |
---|
635 | { int i, j, k, t, len, size, beg, end, min_j, *U_ind, *head, *next, |
---|
636 | *ind, *map, *temp; |
---|
637 | /* initially we assume that on computing the pattern of U fill-in |
---|
638 | will double the number of non-zeros in A */ |
---|
639 | size = A_ptr[n+1] - 1; |
---|
640 | if (size < n) size = n; |
---|
641 | size += size; |
---|
642 | U_ind = xcalloc(1+size, sizeof(int)); |
---|
643 | /* allocate and initialize working arrays */ |
---|
644 | head = xcalloc(1+n, sizeof(int)); |
---|
645 | for (i = 1; i <= n; i++) head[i] = 0; |
---|
646 | next = xcalloc(1+n, sizeof(int)); |
---|
647 | ind = xcalloc(1+n, sizeof(int)); |
---|
648 | map = xcalloc(1+n, sizeof(int)); |
---|
649 | for (j = 1; j <= n; j++) map[j] = 0; |
---|
650 | /* compute the pattern of matrix U */ |
---|
651 | U_ptr[1] = 1; |
---|
652 | for (k = 1; k <= n; k++) |
---|
653 | { /* compute the pattern of k-th row of U, which is the union of |
---|
654 | k-th row of A and those rows of U (among 1, ..., k-1) whose |
---|
655 | leftmost non-diagonal non-zero is placed in k-th column */ |
---|
656 | /* (ind) := (k-th row of A) */ |
---|
657 | len = A_ptr[k+1] - A_ptr[k]; |
---|
658 | memcpy(&ind[1], &A_ind[A_ptr[k]], len * sizeof(int)); |
---|
659 | for (t = 1; t <= len; t++) |
---|
660 | { j = ind[t]; |
---|
661 | xassert(k < j && j <= n); |
---|
662 | map[j] = 1; |
---|
663 | } |
---|
664 | /* walk through rows of U whose leftmost non-diagonal non-zero |
---|
665 | is placed in k-th column */ |
---|
666 | for (i = head[k]; i != 0; i = next[i]) |
---|
667 | { /* (ind) := (ind) union (i-th row of U) */ |
---|
668 | beg = U_ptr[i], end = U_ptr[i+1]; |
---|
669 | for (t = beg; t < end; t++) |
---|
670 | { j = U_ind[t]; |
---|
671 | if (j > k && !map[j]) ind[++len] = j, map[j] = 1; |
---|
672 | } |
---|
673 | } |
---|
674 | /* now (ind) is the pattern of k-th row of U */ |
---|
675 | U_ptr[k+1] = U_ptr[k] + len; |
---|
676 | /* at least (U_ptr[k+1] - 1) locations should be available in |
---|
677 | the array U_ind */ |
---|
678 | if (U_ptr[k+1] - 1 > size) |
---|
679 | { temp = U_ind; |
---|
680 | size += size; |
---|
681 | U_ind = xcalloc(1+size, sizeof(int)); |
---|
682 | memcpy(&U_ind[1], &temp[1], (U_ptr[k] - 1) * sizeof(int)); |
---|
683 | xfree(temp); |
---|
684 | } |
---|
685 | xassert(U_ptr[k+1] - 1 <= size); |
---|
686 | /* (k-th row of U) := (ind) */ |
---|
687 | memcpy(&U_ind[U_ptr[k]], &ind[1], len * sizeof(int)); |
---|
688 | /* determine column index of leftmost non-diagonal non-zero in |
---|
689 | k-th row of U and clear the row pattern map */ |
---|
690 | min_j = n + 1; |
---|
691 | for (t = 1; t <= len; t++) |
---|
692 | { j = ind[t], map[j] = 0; |
---|
693 | if (min_j > j) min_j = j; |
---|
694 | } |
---|
695 | /* include k-th row into corresponding linked list */ |
---|
696 | if (min_j <= n) next[k] = head[min_j], head[min_j] = k; |
---|
697 | } |
---|
698 | /* free working arrays */ |
---|
699 | xfree(head); |
---|
700 | xfree(next); |
---|
701 | xfree(ind); |
---|
702 | xfree(map); |
---|
703 | /* reallocate the array U_ind to free unused locations */ |
---|
704 | temp = U_ind; |
---|
705 | size = U_ptr[n+1] - 1; |
---|
706 | U_ind = xcalloc(1+size, sizeof(int)); |
---|
707 | memcpy(&U_ind[1], &temp[1], size * sizeof(int)); |
---|
708 | xfree(temp); |
---|
709 | return U_ind; |
---|
710 | } |
---|
711 | |
---|
712 | /*---------------------------------------------------------------------- |
---|
713 | -- chol_numeric - compute Cholesky factorization (numeric phase). |
---|
714 | -- |
---|
715 | -- *Synopsis* |
---|
716 | -- |
---|
717 | -- #include "glpmat.h" |
---|
718 | -- int chol_numeric(int n, |
---|
719 | -- int A_ptr[], int A_ind[], double A_val[], double A_diag[], |
---|
720 | -- int U_ptr[], int U_ind[], double U_val[], double U_diag[]); |
---|
721 | -- |
---|
722 | -- *Description* |
---|
723 | -- |
---|
724 | -- The routine chol_symbolic implements the numeric phase of Cholesky |
---|
725 | -- factorization A = U'*U, where A is a given sparse symmetric positive |
---|
726 | -- definite matrix, U is a resultant upper triangular factor, U' is a |
---|
727 | -- matrix transposed to U. |
---|
728 | -- |
---|
729 | -- The parameter n is the order of matrices A and U. |
---|
730 | -- |
---|
731 | -- Upper triangular part of the matrix A without diagonal elements is |
---|
732 | -- specified in the arrays A_ptr, A_ind, and A_val in storage-by-rows |
---|
733 | -- format. Diagonal elements of A are specified in the array A_diag, |
---|
734 | -- where A_diag[0] is not used, A_diag[i] = a[i,i] for i = 1, ..., n. |
---|
735 | -- The arrays A_ptr, A_ind, A_val, and A_diag are not changed on exit. |
---|
736 | -- |
---|
737 | -- The pattern of the matrix U without diagonal elements (previously |
---|
738 | -- computed with the routine chol_symbolic) is specified in the arrays |
---|
739 | -- U_ptr and U_ind, which are not changed on exit. Numeric values of |
---|
740 | -- non-diagonal elements of U are stored in corresponding locations of |
---|
741 | -- the array U_val, and values of diagonal elements of U are stored in |
---|
742 | -- locations U_diag[1], ..., U_diag[n]. |
---|
743 | -- |
---|
744 | -- *Returns* |
---|
745 | -- |
---|
746 | -- The routine returns the number of non-positive diagonal elements of |
---|
747 | -- the matrix U which have been replaced by a huge positive number (see |
---|
748 | -- the method description below). Zero return code means the matrix A |
---|
749 | -- has been successfully factorized. |
---|
750 | -- |
---|
751 | -- *Method* |
---|
752 | -- |
---|
753 | -- The routine chol_numeric computes the matrix U in a row-wise manner |
---|
754 | -- using standard gaussian elimination technique. No pivoting is used. |
---|
755 | -- |
---|
756 | -- Initially the routine sets U = A, and before k-th elimination step |
---|
757 | -- the matrix U is the following: |
---|
758 | -- |
---|
759 | -- 1 k n |
---|
760 | -- 1 x x x x x x x x x x |
---|
761 | -- . x x x x x x x x x |
---|
762 | -- . . x x x x x x x x |
---|
763 | -- . . . x x x x x x x |
---|
764 | -- k . . . . * * * * * * |
---|
765 | -- . . . . * * * * * * |
---|
766 | -- . . . . * * * * * * |
---|
767 | -- . . . . * * * * * * |
---|
768 | -- . . . . * * * * * * |
---|
769 | -- n . . . . * * * * * * |
---|
770 | -- |
---|
771 | -- where 'x' are elements of already computed rows, '*' are elements of |
---|
772 | -- the active submatrix. (Note that the lower triangular part of the |
---|
773 | -- active submatrix being symmetric is not stored and diagonal elements |
---|
774 | -- are stored separately in the array U_diag.) |
---|
775 | -- |
---|
776 | -- The matrix A is assumed to be positive definite. However, if it is |
---|
777 | -- close to semi-definite, on some elimination step a pivot u[k,k] may |
---|
778 | -- happen to be non-positive due to round-off errors. In this case the |
---|
779 | -- routine uses a technique proposed in: |
---|
780 | -- |
---|
781 | -- S.J.Wright. The Cholesky factorization in interior-point and barrier |
---|
782 | -- methods. Preprint MCS-P600-0596, Mathematics and Computer Science |
---|
783 | -- Division, Argonne National Laboratory, Argonne, Ill., May 1996. |
---|
784 | -- |
---|
785 | -- The routine just replaces non-positive u[k,k] by a huge positive |
---|
786 | -- number. This involves non-diagonal elements in k-th row of U to be |
---|
787 | -- close to zero that, in turn, involves k-th component of a solution |
---|
788 | -- vector to be close to zero. Note, however, that this technique works |
---|
789 | -- only if the system A*x = b is consistent. */ |
---|
790 | |
---|
791 | int chol_numeric(int n, |
---|
792 | int A_ptr[], int A_ind[], double A_val[], double A_diag[], |
---|
793 | int U_ptr[], int U_ind[], double U_val[], double U_diag[]) |
---|
794 | { int i, j, k, t, t1, beg, end, beg1, end1, count = 0; |
---|
795 | double ukk, uki, *work; |
---|
796 | work = xcalloc(1+n, sizeof(double)); |
---|
797 | for (j = 1; j <= n; j++) work[j] = 0.0; |
---|
798 | /* U := (upper triangle of A) */ |
---|
799 | /* note that the upper traingle of A is a subset of U */ |
---|
800 | for (i = 1; i <= n; i++) |
---|
801 | { beg = A_ptr[i], end = A_ptr[i+1]; |
---|
802 | for (t = beg; t < end; t++) |
---|
803 | j = A_ind[t], work[j] = A_val[t]; |
---|
804 | beg = U_ptr[i], end = U_ptr[i+1]; |
---|
805 | for (t = beg; t < end; t++) |
---|
806 | j = U_ind[t], U_val[t] = work[j], work[j] = 0.0; |
---|
807 | U_diag[i] = A_diag[i]; |
---|
808 | } |
---|
809 | /* main elimination loop */ |
---|
810 | for (k = 1; k <= n; k++) |
---|
811 | { /* transform k-th row of U */ |
---|
812 | ukk = U_diag[k]; |
---|
813 | if (ukk > 0.0) |
---|
814 | U_diag[k] = ukk = sqrt(ukk); |
---|
815 | else |
---|
816 | U_diag[k] = ukk = DBL_MAX, count++; |
---|
817 | /* (work) := (transformed k-th row) */ |
---|
818 | beg = U_ptr[k], end = U_ptr[k+1]; |
---|
819 | for (t = beg; t < end; t++) |
---|
820 | work[U_ind[t]] = (U_val[t] /= ukk); |
---|
821 | /* transform other rows of U */ |
---|
822 | for (t = beg; t < end; t++) |
---|
823 | { i = U_ind[t]; |
---|
824 | xassert(i > k); |
---|
825 | /* (i-th row) := (i-th row) - u[k,i] * (k-th row) */ |
---|
826 | uki = work[i]; |
---|
827 | beg1 = U_ptr[i], end1 = U_ptr[i+1]; |
---|
828 | for (t1 = beg1; t1 < end1; t1++) |
---|
829 | U_val[t1] -= uki * work[U_ind[t1]]; |
---|
830 | U_diag[i] -= uki * uki; |
---|
831 | } |
---|
832 | /* (work) := 0 */ |
---|
833 | for (t = beg; t < end; t++) |
---|
834 | work[U_ind[t]] = 0.0; |
---|
835 | } |
---|
836 | xfree(work); |
---|
837 | return count; |
---|
838 | } |
---|
839 | |
---|
840 | /*---------------------------------------------------------------------- |
---|
841 | -- u_solve - solve upper triangular system U*x = b. |
---|
842 | -- |
---|
843 | -- *Synopsis* |
---|
844 | -- |
---|
845 | -- #include "glpmat.h" |
---|
846 | -- void u_solve(int n, int U_ptr[], int U_ind[], double U_val[], |
---|
847 | -- double U_diag[], double x[]); |
---|
848 | -- |
---|
849 | -- *Description* |
---|
850 | -- |
---|
851 | -- The routine u_solve solves an linear system U*x = b, where U is an |
---|
852 | -- upper triangular matrix. |
---|
853 | -- |
---|
854 | -- The parameter n is the order of matrix U. |
---|
855 | -- |
---|
856 | -- The matrix U without diagonal elements is specified in the arrays |
---|
857 | -- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements |
---|
858 | -- of U are specified in the array U_diag, where U_diag[0] is not used, |
---|
859 | -- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not |
---|
860 | -- changed on exit. |
---|
861 | -- |
---|
862 | -- The right-hand side vector b is specified on entry in the array x, |
---|
863 | -- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit |
---|
864 | -- the routine stores computed components of the vector of unknowns x |
---|
865 | -- in the array x in the same manner. */ |
---|
866 | |
---|
867 | void u_solve(int n, int U_ptr[], int U_ind[], double U_val[], |
---|
868 | double U_diag[], double x[]) |
---|
869 | { int i, t, beg, end; |
---|
870 | double temp; |
---|
871 | for (i = n; i >= 1; i--) |
---|
872 | { temp = x[i]; |
---|
873 | beg = U_ptr[i], end = U_ptr[i+1]; |
---|
874 | for (t = beg; t < end; t++) |
---|
875 | temp -= U_val[t] * x[U_ind[t]]; |
---|
876 | xassert(U_diag[i] != 0.0); |
---|
877 | x[i] = temp / U_diag[i]; |
---|
878 | } |
---|
879 | return; |
---|
880 | } |
---|
881 | |
---|
882 | /*---------------------------------------------------------------------- |
---|
883 | -- ut_solve - solve lower triangular system U'*x = b. |
---|
884 | -- |
---|
885 | -- *Synopsis* |
---|
886 | -- |
---|
887 | -- #include "glpmat.h" |
---|
888 | -- void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[], |
---|
889 | -- double U_diag[], double x[]); |
---|
890 | -- |
---|
891 | -- *Description* |
---|
892 | -- |
---|
893 | -- The routine ut_solve solves an linear system U'*x = b, where U is a |
---|
894 | -- matrix transposed to an upper triangular matrix. |
---|
895 | -- |
---|
896 | -- The parameter n is the order of matrix U. |
---|
897 | -- |
---|
898 | -- The matrix U without diagonal elements is specified in the arrays |
---|
899 | -- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements |
---|
900 | -- of U are specified in the array U_diag, where U_diag[0] is not used, |
---|
901 | -- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not |
---|
902 | -- changed on exit. |
---|
903 | -- |
---|
904 | -- The right-hand side vector b is specified on entry in the array x, |
---|
905 | -- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit |
---|
906 | -- the routine stores computed components of the vector of unknowns x |
---|
907 | -- in the array x in the same manner. */ |
---|
908 | |
---|
909 | void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[], |
---|
910 | double U_diag[], double x[]) |
---|
911 | { int i, t, beg, end; |
---|
912 | double temp; |
---|
913 | for (i = 1; i <= n; i++) |
---|
914 | { xassert(U_diag[i] != 0.0); |
---|
915 | temp = (x[i] /= U_diag[i]); |
---|
916 | if (temp == 0.0) continue; |
---|
917 | beg = U_ptr[i], end = U_ptr[i+1]; |
---|
918 | for (t = beg; t < end; t++) |
---|
919 | x[U_ind[t]] -= U_val[t] * temp; |
---|
920 | } |
---|
921 | return; |
---|
922 | } |
---|
923 | |
---|
924 | /* eof */ |
---|