1 : /* glpscf.c (Schur complement factorization) */
2 :
3 : /***********************************************************************
4 : * This code is part of GLPK (GNU Linear Programming Kit).
5 : *
6 : * Copyright (C) 2000, 01, 02, 03, 04, 05, 06, 07 Andrew Makhorin,
7 : * Department for Applied Informatics, Moscow Aviation Institute,
8 : * Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
9 : *
10 : * GLPK is free software: you can redistribute it and/or modify it
11 : * under the terms of the GNU General Public License as published by
12 : * the Free Software Foundation, either version 3 of the License, or
13 : * (at your option) any later version.
14 : *
15 : * GLPK is distributed in the hope that it will be useful, but WITHOUT
16 : * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17 : * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
18 : * License for more details.
19 : *
20 : * You should have received a copy of the GNU General Public License
21 : * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
22 : ***********************************************************************/
23 :
24 : #include "glplib.h"
25 : #include "glpscf.h"
26 :
27 : #define _GLPSCF_DEBUG 0
28 :
29 : #define eps 1e-10
30 :
31 : /***********************************************************************
32 : * NAME
33 : *
34 : * scf_create_it - create Schur complement factorization
35 : *
36 : * SYNOPSIS
37 : *
38 : * #include "glpscf.h"
39 : * SCF *scf_create_it(int n_max);
40 : *
41 : * DESCRIPTION
42 : *
43 : * The routine scf_create_it creates the factorization of matrix C,
44 : * which initially has no rows and columns.
45 : *
46 : * The parameter n_max specifies the maximal order of matrix C to be
47 : * factorized, 1 <= n_max <= 32767.
48 : *
49 : * RETURNS
50 : *
51 : * The routine scf_create_it returns a pointer to the structure SCF,
52 : * which defines the factorization. */
53 :
54 : SCF *scf_create_it(int n_max)
55 0 : { SCF *scf;
56 : #if _GLPSCF_DEBUG
57 : xprintf("scf_create_it: warning: debug mode enabled\n");
58 : #endif
59 0 : if (!(1 <= n_max && n_max <= 32767))
60 0 : xfault("scf_create_it: n_max = %d; invalid parameter\n",
61 : n_max);
62 0 : scf = xmalloc(sizeof(SCF));
63 0 : scf->n_max = n_max;
64 0 : scf->n = 0;
65 0 : scf->f = xcalloc(1 + n_max * n_max, sizeof(double));
66 0 : scf->u = xcalloc(1 + n_max * (n_max + 1) / 2, sizeof(double));
67 0 : scf->p = xcalloc(1 + n_max, sizeof(int));
68 0 : scf->t_opt = SCF_TBG;
69 0 : scf->rank = 0;
70 : #if _GLPSCF_DEBUG
71 : scf->c = xcalloc(1 + n_max * n_max, sizeof(double));
72 : #else
73 0 : scf->c = NULL;
74 : #endif
75 0 : scf->w = xcalloc(1 + n_max, sizeof(double));
76 0 : return scf;
77 : }
78 :
79 : /***********************************************************************
80 : * The routine f_loc determines location of matrix element F[i,j] in
81 : * the one-dimensional array f. */
82 :
83 : static int f_loc(SCF *scf, int i, int j)
84 0 : { int n_max = scf->n_max;
85 0 : int n = scf->n;
86 0 : xassert(1 <= i && i <= n);
87 0 : xassert(1 <= j && j <= n);
88 0 : return (i - 1) * n_max + j;
89 : }
90 :
91 : /***********************************************************************
92 : * The routine u_loc determines location of matrix element U[i,j] in
93 : * the one-dimensional array u. */
94 :
95 : static int u_loc(SCF *scf, int i, int j)
96 0 : { int n_max = scf->n_max;
97 0 : int n = scf->n;
98 0 : xassert(1 <= i && i <= n);
99 0 : xassert(i <= j && j <= n);
100 0 : return (i - 1) * n_max + j - i * (i - 1) / 2;
101 : }
102 :
103 : /***********************************************************************
104 : * The routine bg_transform applies Bartels-Golub version of gaussian
105 : * elimination to restore triangular structure of matrix U.
106 : *
107 : * On entry matrix U has the following structure:
108 : *
109 : * 1 k n
110 : * 1 * * * * * * * * * *
111 : * . * * * * * * * * *
112 : * . . * * * * * * * *
113 : * . . . * * * * * * *
114 : * k . . . . * * * * * *
115 : * . . . . . * * * * *
116 : * . . . . . . * * * *
117 : * . . . . . . . * * *
118 : * . . . . . . . . * *
119 : * n . . . . # # # # # #
120 : *
121 : * where '#' is a row spike to be eliminated.
122 : *
123 : * Elements of n-th row are passed separately in locations un[k], ...,
124 : * un[n]. On exit the content of the array un is destroyed.
125 : *
126 : * REFERENCES
127 : *
128 : * R.H.Bartels, G.H.Golub, "The Simplex Method of Linear Programming
129 : * Using LU-decomposition", Comm. ACM, 12, pp. 266-68, 1969. */
130 :
131 : static void bg_transform(SCF *scf, int k, double un[])
132 0 : { int n = scf->n;
133 0 : double *f = scf->f;
134 0 : double *u = scf->u;
135 : int j, k1, kj, kk, n1, nj;
136 : double t;
137 0 : xassert(1 <= k && k <= n);
138 : /* main elimination loop */
139 0 : for (k = k; k < n; k++)
140 : { /* determine location of U[k,k] */
141 0 : kk = u_loc(scf, k, k);
142 : /* determine location of F[k,1] */
143 0 : k1 = f_loc(scf, k, 1);
144 : /* determine location of F[n,1] */
145 0 : n1 = f_loc(scf, n, 1);
146 : /* if |U[k,k]| < |U[n,k]|, interchange k-th and n-th rows to
147 : provide |U[k,k]| >= |U[n,k]| */
148 0 : if (fabs(u[kk]) < fabs(un[k]))
149 : { /* interchange k-th and n-th rows of matrix U */
150 0 : for (j = k, kj = kk; j <= n; j++, kj++)
151 0 : t = u[kj], u[kj] = un[j], un[j] = t;
152 : /* interchange k-th and n-th rows of matrix F to keep the
153 : main equality F * C = U * P */
154 0 : for (j = 1, kj = k1, nj = n1; j <= n; j++, kj++, nj++)
155 0 : t = f[kj], f[kj] = f[nj], f[nj] = t;
156 : }
157 : /* now |U[k,k]| >= |U[n,k]| */
158 : /* if U[k,k] is too small in the magnitude, replace U[k,k] and
159 : U[n,k] by exact zero */
160 0 : if (fabs(u[kk]) < eps) u[kk] = un[k] = 0.0;
161 : /* if U[n,k] is already zero, elimination is not needed */
162 0 : if (un[k] == 0.0) continue;
163 : /* compute gaussian multiplier t = U[n,k] / U[k,k] */
164 0 : t = un[k] / u[kk];
165 : /* apply gaussian elimination to nullify U[n,k] */
166 : /* (n-th row of U) := (n-th row of U) - t * (k-th row of U) */
167 0 : for (j = k+1, kj = kk+1; j <= n; j++, kj++)
168 0 : un[j] -= t * u[kj];
169 : /* (n-th row of F) := (n-th row of F) - t * (k-th row of F)
170 : to keep the main equality F * C = U * P */
171 0 : for (j = 1, kj = k1, nj = n1; j <= n; j++, kj++, nj++)
172 0 : f[nj] -= t * f[kj];
173 : }
174 : /* if U[n,n] is too small in the magnitude, replace it by exact
175 : zero */
176 0 : if (fabs(un[n]) < eps) un[n] = 0.0;
177 : /* store U[n,n] in a proper location */
178 0 : u[u_loc(scf, n, n)] = un[n];
179 : return;
180 : }
181 :
182 : /***********************************************************************
183 : * The routine givens computes the parameters of Givens plane rotation
184 : * c = cos(teta) and s = sin(teta) such that:
185 : *
186 : * ( c -s ) ( a ) ( r )
187 : * ( ) ( ) = ( ) ,
188 : * ( s c ) ( b ) ( 0 )
189 : *
190 : * where a and b are given scalars.
191 : *
192 : * REFERENCES
193 : *
194 : * G.H.Golub, C.F.Van Loan, "Matrix Computations", 2nd ed. */
195 :
196 : static void givens(double a, double b, double *c, double *s)
197 0 : { double t;
198 0 : if (b == 0.0)
199 0 : (*c) = 1.0, (*s) = 0.0;
200 0 : else if (fabs(a) <= fabs(b))
201 0 : t = - a / b, (*s) = 1.0 / sqrt(1.0 + t * t), (*c) = (*s) * t;
202 : else
203 0 : t = - b / a, (*c) = 1.0 / sqrt(1.0 + t * t), (*s) = (*c) * t;
204 : return;
205 : }
206 :
207 : /*----------------------------------------------------------------------
208 : * The routine gr_transform applies Givens plane rotations to restore
209 : * triangular structure of matrix U.
210 : *
211 : * On entry matrix U has the following structure:
212 : *
213 : * 1 k n
214 : * 1 * * * * * * * * * *
215 : * . * * * * * * * * *
216 : * . . * * * * * * * *
217 : * . . . * * * * * * *
218 : * k . . . . * * * * * *
219 : * . . . . . * * * * *
220 : * . . . . . . * * * *
221 : * . . . . . . . * * *
222 : * . . . . . . . . * *
223 : * n . . . . # # # # # #
224 : *
225 : * where '#' is a row spike to be eliminated.
226 : *
227 : * Elements of n-th row are passed separately in locations un[k], ...,
228 : * un[n]. On exit the content of the array un is destroyed.
229 : *
230 : * REFERENCES
231 : *
232 : * R.H.Bartels, G.H.Golub, "The Simplex Method of Linear Programming
233 : * Using LU-decomposition", Comm. ACM, 12, pp. 266-68, 1969. */
234 :
235 : static void gr_transform(SCF *scf, int k, double un[])
236 0 : { int n = scf->n;
237 0 : double *f = scf->f;
238 0 : double *u = scf->u;
239 : int j, k1, kj, kk, n1, nj;
240 : double c, s;
241 0 : xassert(1 <= k && k <= n);
242 : /* main elimination loop */
243 0 : for (k = k; k < n; k++)
244 : { /* determine location of U[k,k] */
245 0 : kk = u_loc(scf, k, k);
246 : /* determine location of F[k,1] */
247 0 : k1 = f_loc(scf, k, 1);
248 : /* determine location of F[n,1] */
249 0 : n1 = f_loc(scf, n, 1);
250 : /* if both U[k,k] and U[n,k] are too small in the magnitude,
251 : replace them by exact zero */
252 0 : if (fabs(u[kk]) < eps && fabs(un[k]) < eps)
253 0 : u[kk] = un[k] = 0.0;
254 : /* if U[n,k] is already zero, elimination is not needed */
255 0 : if (un[k] == 0.0) continue;
256 : /* compute the parameters of Givens plane rotation */
257 0 : givens(u[kk], un[k], &c, &s);
258 : /* apply Givens rotation to k-th and n-th rows of matrix U */
259 0 : for (j = k, kj = kk; j <= n; j++, kj++)
260 0 : { double ukj = u[kj], unj = un[j];
261 0 : u[kj] = c * ukj - s * unj;
262 0 : un[j] = s * ukj + c * unj;
263 : }
264 : /* apply Givens rotation to k-th and n-th rows of matrix F
265 : to keep the main equality F * C = U * P */
266 0 : for (j = 1, kj = k1, nj = n1; j <= n; j++, kj++, nj++)
267 0 : { double fkj = f[kj], fnj = f[nj];
268 0 : f[kj] = c * fkj - s * fnj;
269 0 : f[nj] = s * fkj + c * fnj;
270 : }
271 : }
272 : /* if U[n,n] is too small in the magnitude, replace it by exact
273 : zero */
274 0 : if (fabs(un[n]) < eps) un[n] = 0.0;
275 : /* store U[n,n] in a proper location */
276 0 : u[u_loc(scf, n, n)] = un[n];
277 : return;
278 : }
279 :
280 : /***********************************************************************
281 : * The routine transform restores triangular structure of matrix U.
282 : * It is a driver to the routines bg_transform and gr_transform (see
283 : * comments to these routines above). */
284 :
285 : static void transform(SCF *scf, int k, double un[])
286 0 : { switch (scf->t_opt)
287 : { case SCF_TBG:
288 0 : bg_transform(scf, k, un);
289 0 : break;
290 : case SCF_TGR:
291 0 : gr_transform(scf, k, un);
292 0 : break;
293 : default:
294 0 : xassert(scf != scf);
295 : }
296 : return;
297 : }
298 :
299 : /***********************************************************************
300 : * The routine estimate_rank estimates the rank of matrix C.
301 : *
302 : * Since all transformations applied to matrix F are non-singular,
303 : * and F is assumed to be well conditioned, from the main equaility
304 : * F * C = U * P it follows that rank(C) = rank(U), where rank(U) is
305 : * estimated as the number of non-zero diagonal elements of U. */
306 :
307 : static int estimate_rank(SCF *scf)
308 0 : { int n_max = scf->n_max;
309 0 : int n = scf->n;
310 0 : double *u = scf->u;
311 0 : int i, ii, inc, rank = 0;
312 0 : for (i = 1, ii = u_loc(scf, i, i), inc = n_max; i <= n;
313 0 : i++, ii += inc, inc--)
314 0 : if (u[ii] != 0.0) rank++;
315 0 : return rank;
316 : }
317 :
318 : #if _GLPSCF_DEBUG
319 : /***********************************************************************
320 : * The routine check_error computes the maximal relative error between
321 : * left- and right-hand sides of the main equality F * C = U * P. (This
322 : * routine is intended only for debugging.) */
323 :
324 : static void check_error(SCF *scf, const char *func)
325 : { int n = scf->n;
326 : double *f = scf->f;
327 : double *u = scf->u;
328 : int *p = scf->p;
329 : double *c = scf->c;
330 : int i, j, k;
331 : double d, dmax = 0.0, s, t;
332 : xassert(c != NULL);
333 : for (i = 1; i <= n; i++)
334 : { for (j = 1; j <= n; j++)
335 : { /* compute element (i,j) of product F * C */
336 : s = 0.0;
337 : for (k = 1; k <= n; k++)
338 : s += f[f_loc(scf, i, k)] * c[f_loc(scf, k, j)];
339 : /* compute element (i,j) of product U * P */
340 : k = p[j];
341 : t = (i <= k ? u[u_loc(scf, i, k)] : 0.0);
342 : /* compute the maximal relative error */
343 : d = fabs(s - t) / (1.0 + fabs(t));
344 : if (dmax < d) dmax = d;
345 : }
346 : }
347 : if (dmax > 1e-8)
348 : xprintf("%s: dmax = %g; relative error too large\n", func,
349 : dmax);
350 : return;
351 : }
352 : #endif
353 :
354 : /***********************************************************************
355 : * NAME
356 : *
357 : * scf_update_exp - update factorization on expanding C
358 : *
359 : * SYNOPSIS
360 : *
361 : * #include "glpscf.h"
362 : * int scf_update_exp(SCF *scf, const double x[], const double y[],
363 : * double z);
364 : *
365 : * DESCRIPTION
366 : *
367 : * The routine scf_update_exp updates the factorization of matrix C on
368 : * expanding it by adding a new row and column as follows:
369 : *
370 : * ( C x )
371 : * new C = ( )
372 : * ( y' z )
373 : *
374 : * where x[1,...,n] is a new column, y[1,...,n] is a new row, and z is
375 : * a new diagonal element.
376 : *
377 : * If on entry the factorization is empty, the parameters x and y can
378 : * be specified as NULL.
379 : *
380 : * RETURNS
381 : *
382 : * 0 The factorization has been successfully updated.
383 : *
384 : * SCF_ESING
385 : * The factorization has been successfully updated, however, new
386 : * matrix C is singular within working precision. Note that the new
387 : * factorization remains valid.
388 : *
389 : * SCF_ELIMIT
390 : * There is not enough room to expand the factorization, because
391 : * n = n_max. The factorization remains unchanged.
392 : *
393 : * ALGORITHM
394 : *
395 : * We can see that:
396 : *
397 : * ( F 0 ) ( C x ) ( FC Fx ) ( UP Fx )
398 : * ( ) ( ) = ( ) = ( ) =
399 : * ( 0 1 ) ( y' z ) ( y' z ) ( y' z )
400 : *
401 : * ( U Fx ) ( P 0 )
402 : * = ( ) ( ),
403 : * ( y'P' z ) ( 0 1 )
404 : *
405 : * therefore to keep the main equality F * C = U * P we can take:
406 : *
407 : * ( F 0 ) ( U Fx ) ( P 0 )
408 : * new F = ( ), new U = ( ), new P = ( ),
409 : * ( 0 1 ) ( y'P' z ) ( 0 1 )
410 : *
411 : * and eliminate the row spike y'P' in the last row of new U to restore
412 : * its upper triangular structure. */
413 :
414 : int scf_update_exp(SCF *scf, const double x[], const double y[],
415 : double z)
416 0 : { int n_max = scf->n_max;
417 0 : int n = scf->n;
418 0 : double *f = scf->f;
419 0 : double *u = scf->u;
420 0 : int *p = scf->p;
421 : #if _GLPSCF_DEBUG
422 : double *c = scf->c;
423 : #endif
424 0 : double *un = scf->w;
425 0 : int i, ij, in, j, k, nj, ret = 0;
426 : double t;
427 : /* check if the factorization can be expanded */
428 0 : if (n == n_max)
429 : { /* there is not enough room */
430 0 : ret = SCF_ELIMIT;
431 0 : goto done;
432 : }
433 : /* increase the order of the factorization */
434 0 : scf->n = ++n;
435 : /* fill new zero column of matrix F */
436 0 : for (i = 1, in = f_loc(scf, i, n); i < n; i++, in += n_max)
437 0 : f[in] = 0.0;
438 : /* fill new zero row of matrix F */
439 0 : for (j = 1, nj = f_loc(scf, n, j); j < n; j++, nj++)
440 0 : f[nj] = 0.0;
441 : /* fill new unity diagonal element of matrix F */
442 0 : f[f_loc(scf, n, n)] = 1.0;
443 : /* compute new column of matrix U, which is (old F) * x */
444 0 : for (i = 1; i < n; i++)
445 : { /* u[i,n] := (i-th row of old F) * x */
446 0 : t = 0.0;
447 0 : for (j = 1, ij = f_loc(scf, i, 1); j < n; j++, ij++)
448 0 : t += f[ij] * x[j];
449 0 : u[u_loc(scf, i, n)] = t;
450 : }
451 : /* compute new (spiked) row of matrix U, which is (old P) * y */
452 0 : for (j = 1; j < n; j++) un[j] = y[p[j]];
453 : /* store new diagonal element of matrix U, which is z */
454 0 : un[n] = z;
455 : /* expand matrix P */
456 0 : p[n] = n;
457 : #if _GLPSCF_DEBUG
458 : /* expand matrix C */
459 : /* fill its new column, which is x */
460 : for (i = 1, in = f_loc(scf, i, n); i < n; i++, in += n_max)
461 : c[in] = x[i];
462 : /* fill its new row, which is y */
463 : for (j = 1, nj = f_loc(scf, n, j); j < n; j++, nj++)
464 : c[nj] = y[j];
465 : /* fill its new diagonal element, which is z */
466 : c[f_loc(scf, n, n)] = z;
467 : #endif
468 : /* restore upper triangular structure of matrix U */
469 0 : for (k = 1; k < n; k++)
470 0 : if (un[k] != 0.0) break;
471 0 : transform(scf, k, un);
472 : /* estimate the rank of matrices C and U */
473 0 : scf->rank = estimate_rank(scf);
474 0 : if (scf->rank != n) ret = SCF_ESING;
475 : #if _GLPSCF_DEBUG
476 : /* check that the factorization is accurate enough */
477 : check_error(scf, "scf_update_exp");
478 : #endif
479 0 : done: return ret;
480 : }
481 :
482 : /***********************************************************************
483 : * The routine solve solves the system C * x = b.
484 : *
485 : * From the main equation F * C = U * P it follows that:
486 : *
487 : * C * x = b => F * C * x = F * b => U * P * x = F * b =>
488 : *
489 : * P * x = inv(U) * F * b => x = P' * inv(U) * F * b.
490 : *
491 : * On entry the array x contains right-hand side vector b. On exit this
492 : * array contains solution vector x. */
493 :
494 : static void solve(SCF *scf, double x[])
495 0 : { int n = scf->n;
496 0 : double *f = scf->f;
497 0 : double *u = scf->u;
498 0 : int *p = scf->p;
499 0 : double *y = scf->w;
500 : int i, j, ij;
501 : double t;
502 : /* y := F * b */
503 0 : for (i = 1; i <= n; i++)
504 : { /* y[i] = (i-th row of F) * b */
505 0 : t = 0.0;
506 0 : for (j = 1, ij = f_loc(scf, i, 1); j <= n; j++, ij++)
507 0 : t += f[ij] * x[j];
508 0 : y[i] = t;
509 : }
510 : /* y := inv(U) * y */
511 0 : for (i = n; i >= 1; i--)
512 0 : { t = y[i];
513 0 : for (j = n, ij = u_loc(scf, i, n); j > i; j--, ij--)
514 0 : t -= u[ij] * y[j];
515 0 : y[i] = t / u[ij];
516 : }
517 : /* x := P' * y */
518 0 : for (i = 1; i <= n; i++) x[p[i]] = y[i];
519 : return;
520 : }
521 :
522 : /***********************************************************************
523 : * The routine tsolve solves the transposed system C' * x = b.
524 : *
525 : * From the main equation F * C = U * P it follows that:
526 : *
527 : * C' * F' = P' * U',
528 : *
529 : * therefore:
530 : *
531 : * C' * x = b => C' * F' * inv(F') * x = b =>
532 : *
533 : * P' * U' * inv(F') * x = b => U' * inv(F') * x = P * b =>
534 : *
535 : * inv(F') * x = inv(U') * P * b => x = F' * inv(U') * P * b.
536 : *
537 : * On entry the array x contains right-hand side vector b. On exit this
538 : * array contains solution vector x. */
539 :
540 : static void tsolve(SCF *scf, double x[])
541 0 : { int n = scf->n;
542 0 : double *f = scf->f;
543 0 : double *u = scf->u;
544 0 : int *p = scf->p;
545 0 : double *y = scf->w;
546 : int i, j, ij;
547 : double t;
548 : /* y := P * b */
549 0 : for (i = 1; i <= n; i++) y[i] = x[p[i]];
550 : /* y := inv(U') * y */
551 0 : for (i = 1; i <= n; i++)
552 : { /* compute y[i] */
553 0 : ij = u_loc(scf, i, i);
554 0 : t = (y[i] /= u[ij]);
555 : /* substitute y[i] in other equations */
556 0 : for (j = i+1, ij++; j <= n; j++, ij++)
557 0 : y[j] -= u[ij] * t;
558 : }
559 : /* x := F' * y (computed as linear combination of rows of F) */
560 0 : for (j = 1; j <= n; j++) x[j] = 0.0;
561 0 : for (i = 1; i <= n; i++)
562 0 : { t = y[i]; /* coefficient of linear combination */
563 0 : for (j = 1, ij = f_loc(scf, i, 1); j <= n; j++, ij++)
564 0 : x[j] += f[ij] * t;
565 : }
566 : return;
567 : }
568 :
569 : /***********************************************************************
570 : * NAME
571 : *
572 : * scf_solve_it - solve either system C * x = b or C' * x = b
573 : *
574 : * SYNOPSIS
575 : *
576 : * #include "glpscf.h"
577 : * void scf_solve_it(SCF *scf, int tr, double x[]);
578 : *
579 : * DESCRIPTION
580 : *
581 : * The routine scf_solve_it solves either the system C * x = b (if tr
582 : * is zero) or the system C' * x = b, where C' is a matrix transposed
583 : * to C (if tr is non-zero). C is assumed to be non-singular.
584 : *
585 : * On entry the array x should contain the right-hand side vector b in
586 : * locations x[1], ..., x[n], where n is the order of matrix C. On exit
587 : * the array x contains the solution vector x in the same locations. */
588 :
589 : void scf_solve_it(SCF *scf, int tr, double x[])
590 0 : { if (scf->rank < scf->n)
591 0 : xfault("scf_solve_it: singular matrix\n");
592 0 : if (!tr)
593 0 : solve(scf, x);
594 : else
595 0 : tsolve(scf, x);
596 : return;
597 : }
598 :
599 : void scf_reset_it(SCF *scf)
600 0 : { /* reset factorization for empty matrix C */
601 0 : scf->n = scf->rank = 0;
602 : return;
603 : }
604 :
605 : /***********************************************************************
606 : * NAME
607 : *
608 : * scf_delete_it - delete Schur complement factorization
609 : *
610 : * SYNOPSIS
611 : *
612 : * #include "glpscf.h"
613 : * void scf_delete_it(SCF *scf);
614 : *
615 : * DESCRIPTION
616 : *
617 : * The routine scf_delete_it deletes the specified factorization and
618 : * frees all the memory allocated to this object. */
619 :
620 : void scf_delete_it(SCF *scf)
621 0 : { xfree(scf->f);
622 0 : xfree(scf->u);
623 0 : xfree(scf->p);
624 : #if _GLPSCF_DEBUG
625 : xfree(scf->c);
626 : #endif
627 0 : xfree(scf->w);
628 0 : xfree(scf);
629 : return;
630 : }
631 :
632 : /* eof */
|
