1 /* scf.c (sparse updatable Schur-complement-based factorization) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 * Copyright (C) 2013-2014 Free Software Foundation, Inc.
6 * Written by Andrew Makhorin <mao@gnu.org>.
7 *
8 * GLPK is free software: you can redistribute it and/or modify it
9 * under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * GLPK is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
20 ***********************************************************************/
21
22 #include "env.h"
23 #include "scf.h"
24
25 /***********************************************************************
26 * scf_r0_solve - solve system R0 * x = b or R0'* x = b
27 *
28 * This routine solves the system R0 * x = b (if tr is zero) or the
29 * system R0'* x = b (if tr is non-zero), where R0 is the left factor
30 * of the initial matrix A0 = R0 * S0.
31 *
32 * On entry the array x should contain elements of the right-hand side
33 * vector b in locations x[1], ..., x[n0], where n0 is the order of the
34 * matrix R0. On exit the array x will contain elements of the solution
35 * vector in the same locations. */
36
scf_r0_solve(SCF * scf,int tr,double x[])37 void scf_r0_solve(SCF *scf, int tr, double x[/*1+n0*/])
38 { switch (scf->type)
39 { case 1:
40 /* A0 = F0 * V0, so R0 = F0 */
41 if (!tr)
42 luf_f_solve(scf->a0.luf, x);
43 else
44 luf_ft_solve(scf->a0.luf, x);
45 break;
46 case 2:
47 /* A0 = I * A0, so R0 = I */
48 break;
49 default:
50 xassert(scf != scf);
51 }
52 return;
53 }
54
55 /***********************************************************************
56 * scf_s0_solve - solve system S0 * x = b or S0'* x = b
57 *
58 * This routine solves the system S0 * x = b (if tr is zero) or the
59 * system S0'* x = b (if tr is non-zero), where S0 is the right factor
60 * of the initial matrix A0 = R0 * S0.
61 *
62 * On entry the array x should contain elements of the right-hand side
63 * vector b in locations x[1], ..., x[n0], where n0 is the order of the
64 * matrix S0. On exit the array x will contain elements of the solution
65 * vector in the same locations.
66 *
67 * The routine uses locations [1], ..., [n0] of three working arrays
68 * w1, w2, and w3. (In case of type = 1 arrays w2 and w3 are not used
69 * and can be specified as NULL.) */
70
scf_s0_solve(SCF * scf,int tr,double x[],double w1[],double w2[],double w3[])71 void scf_s0_solve(SCF *scf, int tr, double x[/*1+n0*/],
72 double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/])
73 { int n0 = scf->n0;
74 switch (scf->type)
75 { case 1:
76 /* A0 = F0 * V0, so S0 = V0 */
77 if (!tr)
78 luf_v_solve(scf->a0.luf, x, w1);
79 else
80 luf_vt_solve(scf->a0.luf, x, w1);
81 break;
82 case 2:
83 /* A0 = I * A0, so S0 = A0 */
84 if (!tr)
85 btf_a_solve(scf->a0.btf, x, w1, w2, w3);
86 else
87 btf_at_solve(scf->a0.btf, x, w1, w2, w3);
88 break;
89 default:
90 xassert(scf != scf);
91 }
92 memcpy(&x[1], &w1[1], n0 * sizeof(double));
93 return;
94 }
95
96 /***********************************************************************
97 * scf_r_prod - compute product y := y + alpha * R * x
98 *
99 * This routine computes the product y := y + alpha * R * x, where
100 * x is a n0-vector, alpha is a scalar, y is a nn-vector.
101 *
102 * Since matrix R is available by rows, the product components are
103 * computed as inner products:
104 *
105 * y[i] = y[i] + alpha * (i-th row of R) * x
106 *
107 * for i = 1, 2, ..., nn. */
108
scf_r_prod(SCF * scf,double y[],double a,const double x[])109 void scf_r_prod(SCF *scf, double y[/*1+nn*/], double a, const double
110 x[/*1+n0*/])
111 { int nn = scf->nn;
112 SVA *sva = scf->sva;
113 int *sv_ind = sva->ind;
114 double *sv_val = sva->val;
115 int rr_ref = scf->rr_ref;
116 int *rr_ptr = &sva->ptr[rr_ref-1];
117 int *rr_len = &sva->len[rr_ref-1];
118 int i, ptr, end;
119 double t;
120 for (i = 1; i <= nn; i++)
121 { /* t := (i-th row of R) * x */
122 t = 0.0;
123 for (end = (ptr = rr_ptr[i]) + rr_len[i]; ptr < end; ptr++)
124 t += sv_val[ptr] * x[sv_ind[ptr]];
125 /* y[i] := y[i] + alpha * t */
126 y[i] += a * t;
127 }
128 return;
129 }
130
131 /***********************************************************************
132 * scf_rt_prod - compute product y := y + alpha * R'* x
133 *
134 * This routine computes the product y := y + alpha * R'* x, where
135 * R' is a matrix transposed to R, x is a nn-vector, alpha is a scalar,
136 * y is a n0-vector.
137 *
138 * Since matrix R is available by rows, the product is computed as a
139 * linear combination:
140 *
141 * y := y + alpha * (R'[1] * x[1] + ... + R'[nn] * x[nn]),
142 *
143 * where R'[i] is i-th row of R. */
144
scf_rt_prod(SCF * scf,double y[],double a,const double x[])145 void scf_rt_prod(SCF *scf, double y[/*1+n0*/], double a, const double
146 x[/*1+nn*/])
147 { int nn = scf->nn;
148 SVA *sva = scf->sva;
149 int *sv_ind = sva->ind;
150 double *sv_val = sva->val;
151 int rr_ref = scf->rr_ref;
152 int *rr_ptr = &sva->ptr[rr_ref-1];
153 int *rr_len = &sva->len[rr_ref-1];
154 int i, ptr, end;
155 double t;
156 for (i = 1; i <= nn; i++)
157 { if (x[i] == 0.0)
158 continue;
159 /* y := y + alpha * R'[i] * x[i] */
160 t = a * x[i];
161 for (end = (ptr = rr_ptr[i]) + rr_len[i]; ptr < end; ptr++)
162 y[sv_ind[ptr]] += sv_val[ptr] * t;
163 }
164 return;
165 }
166
167 /***********************************************************************
168 * scf_s_prod - compute product y := y + alpha * S * x
169 *
170 * This routine computes the product y := y + alpha * S * x, where
171 * x is a nn-vector, alpha is a scalar, y is a n0 vector.
172 *
173 * Since matrix S is available by columns, the product is computed as
174 * a linear combination:
175 *
176 * y := y + alpha * (S[1] * x[1] + ... + S[nn] * x[nn]),
177 *
178 * where S[j] is j-th column of S. */
179
scf_s_prod(SCF * scf,double y[],double a,const double x[])180 void scf_s_prod(SCF *scf, double y[/*1+n0*/], double a, const double
181 x[/*1+nn*/])
182 { int nn = scf->nn;
183 SVA *sva = scf->sva;
184 int *sv_ind = sva->ind;
185 double *sv_val = sva->val;
186 int ss_ref = scf->ss_ref;
187 int *ss_ptr = &sva->ptr[ss_ref-1];
188 int *ss_len = &sva->len[ss_ref-1];
189 int j, ptr, end;
190 double t;
191 for (j = 1; j <= nn; j++)
192 { if (x[j] == 0.0)
193 continue;
194 /* y := y + alpha * S[j] * x[j] */
195 t = a * x[j];
196 for (end = (ptr = ss_ptr[j]) + ss_len[j]; ptr < end; ptr++)
197 y[sv_ind[ptr]] += sv_val[ptr] * t;
198 }
199 return;
200 }
201
202 /***********************************************************************
203 * scf_st_prod - compute product y := y + alpha * S'* x
204 *
205 * This routine computes the product y := y + alpha * S'* x, where
206 * S' is a matrix transposed to S, x is a n0-vector, alpha is a scalar,
207 * y is a nn-vector.
208 *
209 * Since matrix S is available by columns, the product components are
210 * computed as inner products:
211 *
212 * y[j] := y[j] + alpha * (j-th column of S) * x
213 *
214 * for j = 1, 2, ..., nn. */
215
scf_st_prod(SCF * scf,double y[],double a,const double x[])216 void scf_st_prod(SCF *scf, double y[/*1+nn*/], double a, const double
217 x[/*1+n0*/])
218 { int nn = scf->nn;
219 SVA *sva = scf->sva;
220 int *sv_ind = sva->ind;
221 double *sv_val = sva->val;
222 int ss_ref = scf->ss_ref;
223 int *ss_ptr = &sva->ptr[ss_ref-1];
224 int *ss_len = &sva->len[ss_ref-1];
225 int j, ptr, end;
226 double t;
227 for (j = 1; j <= nn; j++)
228 { /* t := (j-th column of S) * x */
229 t = 0.0;
230 for (end = (ptr = ss_ptr[j]) + ss_len[j]; ptr < end; ptr++)
231 t += sv_val[ptr] * x[sv_ind[ptr]];
232 /* y[j] := y[j] + alpha * t */
233 y[j] += a * t;
234 }
235 return;
236 }
237
238 /***********************************************************************
239 * scf_a_solve - solve system A * x = b
240 *
241 * This routine solves the system A * x = b, where A is the current
242 * matrix.
243 *
244 * On entry the array x should contain elements of the right-hand side
245 * vector b in locations x[1], ..., x[n], where n is the order of the
246 * matrix A. On exit the array x will contain elements of the solution
247 * vector in the same locations.
248 *
249 * For details see the program documentation. */
250
scf_a_solve(SCF * scf,double x[],double w[],double work1[],double work2[],double work3[])251 void scf_a_solve(SCF *scf, double x[/*1+n*/],
252 double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/],
253 double work2[/*1+n*/], double work3[/*1+n*/])
254 { int n = scf->n;
255 int n0 = scf->n0;
256 int nn = scf->nn;
257 int *pp_ind = scf->pp_ind;
258 int *qq_inv = scf->qq_inv;
259 int i, ii;
260 /* (u1, u2) := inv(P) * (b, 0) */
261 for (ii = 1; ii <= n0+nn; ii++)
262 { i = pp_ind[ii];
263 #if 1 /* FIXME: currently P = I */
264 xassert(i == ii);
265 #endif
266 w[ii] = (i <= n ? x[i] : 0.0);
267 }
268 /* v1 := inv(R0) * u1 */
269 scf_r0_solve(scf, 0, &w[0]);
270 /* v2 := u2 - R * v1 */
271 scf_r_prod(scf, &w[n0], -1.0, &w[0]);
272 /* w2 := inv(C) * v2 */
273 ifu_a_solve(&scf->ifu, &w[n0], work1);
274 /* w1 := inv(S0) * (v1 - S * w2) */
275 scf_s_prod(scf, &w[0], -1.0, &w[n0]);
276 scf_s0_solve(scf, 0, &w[0], work1, work2, work3);
277 /* (x, x~) := inv(Q) * (w1, w2); x~ is not needed */
278 for (i = 1; i <= n; i++)
279 x[i] = w[qq_inv[i]];
280 return;
281 }
282
283 /***********************************************************************
284 * scf_at_solve - solve system A'* x = b
285 *
286 * This routine solves the system A'* x = b, where A' is a matrix
287 * transposed to the current matrix A.
288 *
289 * On entry the array x should contain elements of the right-hand side
290 * vector b in locations x[1], ..., x[n], where n is the order of the
291 * matrix A. On exit the array x will contain elements of the solution
292 * vector in the same locations.
293 *
294 * For details see the program documentation. */
295
scf_at_solve(SCF * scf,double x[],double w[],double work1[],double work2[],double work3[])296 void scf_at_solve(SCF *scf, double x[/*1+n*/],
297 double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/],
298 double work2[/*1+n*/], double work3[/*1+n*/])
299 { int n = scf->n;
300 int n0 = scf->n0;
301 int nn = scf->nn;
302 int *pp_inv = scf->pp_inv;
303 int *qq_ind = scf->qq_ind;
304 int i, ii;
305 /* (u1, u2) := Q * (b, 0) */
306 for (ii = 1; ii <= n0+nn; ii++)
307 { i = qq_ind[ii];
308 w[ii] = (i <= n ? x[i] : 0.0);
309 }
310 /* v1 := inv(S0') * u1 */
311 scf_s0_solve(scf, 1, &w[0], work1, work2, work3);
312 /* v2 := inv(C') * (u2 - S'* v1) */
313 scf_st_prod(scf, &w[n0], -1.0, &w[0]);
314 ifu_at_solve(&scf->ifu, &w[n0], work1);
315 /* w2 := v2 */
316 /* nop */
317 /* w1 := inv(R0') * (v1 - R'* w2) */
318 scf_rt_prod(scf, &w[0], -1.0, &w[n0]);
319 scf_r0_solve(scf, 1, &w[0]);
320 /* compute (x, x~) := P * (w1, w2); x~ is not needed */
321 for (i = 1; i <= n; i++)
322 {
323 #if 1 /* FIXME: currently P = I */
324 xassert(pp_inv[i] == i);
325 #endif
326 x[i] = w[pp_inv[i]];
327 }
328 return;
329 }
330
331 /***********************************************************************
332 * scf_add_r_row - add new row to matrix R
333 *
334 * This routine adds new (nn+1)-th row to matrix R, whose elements are
335 * specified in locations w[1,...,n0]. */
336
scf_add_r_row(SCF * scf,const double w[])337 void scf_add_r_row(SCF *scf, const double w[/*1+n0*/])
338 { int n0 = scf->n0;
339 int nn = scf->nn;
340 SVA *sva = scf->sva;
341 int *sv_ind = sva->ind;
342 double *sv_val = sva->val;
343 int rr_ref = scf->rr_ref;
344 int *rr_ptr = &sva->ptr[rr_ref-1];
345 int *rr_len = &sva->len[rr_ref-1];
346 int j, len, ptr;
347 xassert(0 <= nn && nn < scf->nn_max);
348 /* determine length of new row */
349 len = 0;
350 for (j = 1; j <= n0; j++)
351 { if (w[j] != 0.0)
352 len++;
353 }
354 /* reserve locations for new row in static part of SVA */
355 if (len > 0)
356 { if (sva->r_ptr - sva->m_ptr < len)
357 { sva_more_space(sva, len);
358 sv_ind = sva->ind;
359 sv_val = sva->val;
360 }
361 sva_reserve_cap(sva, rr_ref + nn, len);
362 }
363 /* store new row in sparse format */
364 ptr = rr_ptr[nn+1];
365 for (j = 1; j <= n0; j++)
366 { if (w[j] != 0.0)
367 { sv_ind[ptr] = j;
368 sv_val[ptr] = w[j];
369 ptr++;
370 }
371 }
372 xassert(ptr - rr_ptr[nn+1] == len);
373 rr_len[nn+1] = len;
374 #ifdef GLP_DEBUG
375 sva_check_area(sva);
376 #endif
377 return;
378 }
379
380 /***********************************************************************
381 * scf_add_s_col - add new column to matrix S
382 *
383 * This routine adds new (nn+1)-th column to matrix S, whose elements
384 * are specified in locations v[1,...,n0]. */
385
scf_add_s_col(SCF * scf,const double v[])386 void scf_add_s_col(SCF *scf, const double v[/*1+n0*/])
387 { int n0 = scf->n0;
388 int nn = scf->nn;
389 SVA *sva = scf->sva;
390 int *sv_ind = sva->ind;
391 double *sv_val = sva->val;
392 int ss_ref = scf->ss_ref;
393 int *ss_ptr = &sva->ptr[ss_ref-1];
394 int *ss_len = &sva->len[ss_ref-1];
395 int i, len, ptr;
396 xassert(0 <= nn && nn < scf->nn_max);
397 /* determine length of new column */
398 len = 0;
399 for (i = 1; i <= n0; i++)
400 { if (v[i] != 0.0)
401 len++;
402 }
403 /* reserve locations for new column in static part of SVA */
404 if (len > 0)
405 { if (sva->r_ptr - sva->m_ptr < len)
406 { sva_more_space(sva, len);
407 sv_ind = sva->ind;
408 sv_val = sva->val;
409 }
410 sva_reserve_cap(sva, ss_ref + nn, len);
411 }
412 /* store new column in sparse format */
413 ptr = ss_ptr[nn+1];
414 for (i = 1; i <= n0; i++)
415 { if (v[i] != 0.0)
416 { sv_ind[ptr] = i;
417 sv_val[ptr] = v[i];
418 ptr++;
419 }
420 }
421 xassert(ptr - ss_ptr[nn+1] == len);
422 ss_len[nn+1] = len;
423 #ifdef GLP_DEBUG
424 sva_check_area(sva);
425 #endif
426 return;
427 }
428
429 /***********************************************************************
430 * scf_update_aug - update factorization of augmented matrix
431 *
432 * Given factorization of the current augmented matrix:
433 *
434 * ( A0 A1 ) ( R0 ) ( S0 S )
435 * ( ) = ( ) ( ),
436 * ( A2 A3 ) ( R I ) ( C )
437 *
438 * this routine computes factorization of the new augmented matrix:
439 *
440 * ( A0 | A1 b )
441 * ( ---+------ ) ( A0 A1^ ) ( R0 ) ( S0 S^ )
442 * ( A2 | A3 f ) = ( ) = ( ) ( ),
443 * ( | ) ( A2^ A3^ ) ( R^ I ) ( C^ )
444 * ( d' | g' h )
445 *
446 * where b and d are specified n0-vectors, f and g are specified
447 * nn-vectors, and h is a specified scalar. (Note that corresponding
448 * arrays are clobbered on exit.)
449 *
450 * The parameter upd specifies how to update factorization of the Schur
451 * complement C:
452 *
453 * 1 Bartels-Golub updating.
454 *
455 * 2 Givens rotations updating.
456 *
457 * The working arrays w1, w2, and w3 are used in the same way as in the
458 * routine scf_s0_solve.
459 *
460 * RETURNS
461 *
462 * 0 Factorization has been successfully updated.
463 *
464 * 1 Updating limit has been reached.
465 *
466 * 2 Updating IFU-factorization of matrix C failed.
467 *
468 * For details see the program documentation. */
469
scf_update_aug(SCF * scf,double b[],double d[],double f[],double g[],double h,int upd,double w1[],double w2[],double w3[])470 int scf_update_aug(SCF *scf, double b[/*1+n0*/], double d[/*1+n0*/],
471 double f[/*1+nn*/], double g[/*1+nn*/], double h, int upd,
472 double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/])
473 { int n0 = scf->n0;
474 int k, ret;
475 double *v, *w, *x, *y, z;
476 if (scf->nn == scf->nn_max)
477 { /* updating limit has been reached */
478 return 1;
479 }
480 /* v := inv(R0) * b */
481 scf_r0_solve(scf, 0, (v = b));
482 /* w := inv(S0') * d */
483 scf_s0_solve(scf, 1, (w = d), w1, w2, w3);
484 /* x := f - R * v */
485 scf_r_prod(scf, (x = f), -1.0, v);
486 /* y := g - S'* w */
487 scf_st_prod(scf, (y = g), -1.0, w);
488 /* z := h - v'* w */
489 z = h;
490 for (k = 1; k <= n0; k++)
491 z -= v[k] * w[k];
492 /* new R := R with row w added */
493 scf_add_r_row(scf, w);
494 /* new S := S with column v added */
495 scf_add_s_col(scf, v);
496 /* update IFU-factorization of C */
497 switch (upd)
498 { case 1:
499 ret = ifu_bg_update(&scf->ifu, x, y, z);
500 break;
501 case 2:
502 ret = ifu_gr_update(&scf->ifu, x, y, z);
503 break;
504 default:
505 xassert(upd != upd);
506 }
507 if (ret != 0)
508 { /* updating IFU-factorization failed */
509 return 2;
510 }
511 /* increase number of additional rows and columns */
512 scf->nn++;
513 /* expand P and Q */
514 k = n0 + scf->nn;
515 scf->pp_ind[k] = scf->pp_inv[k] = k;
516 scf->qq_ind[k] = scf->qq_inv[k] = k;
517 /* factorization has been successfully updated */
518 return 0;
519 }
520
521 /* eof */
522