COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glplpf.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 30.7 KB
Line 
1/* glplpf.c (LP basis factorization, Schur complement version) */
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 "glplpf.h"
26#include "glpenv.h"
27#define xfault xerror
28
29#define _GLPLPF_DEBUG 0
30
31/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
32
33#define M_MAX 100000000 /* = 100*10^6 */
34/* maximal order of the basis matrix */
35
36/***********************************************************************
37*  NAME
38*
39*  lpf_create_it - create LP basis factorization
40*
41*  SYNOPSIS
42*
43*  #include "glplpf.h"
44*  LPF *lpf_create_it(void);
45*
46*  DESCRIPTION
47*
48*  The routine lpf_create_it creates a program object, which represents
49*  a factorization of LP basis.
50*
51*  RETURNS
52*
53*  The routine lpf_create_it returns a pointer to the object created. */
54
55LPF *lpf_create_it(void)
56{     LPF *lpf;
57#if _GLPLPF_DEBUG
58      xprintf("lpf_create_it: warning: debug mode enabled\n");
59#endif
60      lpf = xmalloc(sizeof(LPF));
61      lpf->valid = 0;
62      lpf->m0_max = lpf->m0 = 0;
63      lpf->luf = luf_create_it();
64      lpf->m = 0;
65      lpf->B = NULL;
66      lpf->n_max = 50;
67      lpf->n = 0;
68      lpf->R_ptr = lpf->R_len = NULL;
69      lpf->S_ptr = lpf->S_len = NULL;
70      lpf->scf = NULL;
71      lpf->P_row = lpf->P_col = NULL;
72      lpf->Q_row = lpf->Q_col = NULL;
73      lpf->v_size = 1000;
74      lpf->v_ptr = 0;
75      lpf->v_ind = NULL;
76      lpf->v_val = NULL;
77      lpf->work1 = lpf->work2 = NULL;
78      return lpf;
79}
80
81/***********************************************************************
82*  NAME
83*
84*  lpf_factorize - compute LP basis factorization
85*
86*  SYNOPSIS
87*
88*  #include "glplpf.h"
89*  int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
90*     (void *info, int j, int ind[], double val[]), void *info);
91*
92*  DESCRIPTION
93*
94*  The routine lpf_factorize computes the factorization of the basis
95*  matrix B specified by the routine col.
96*
97*  The parameter lpf specified the basis factorization data structure
98*  created with the routine lpf_create_it.
99*
100*  The parameter m specifies the order of B, m > 0.
101*
102*  The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
103*  number of j-th column of B in some original matrix. The array bh is
104*  optional and can be specified as NULL.
105*
106*  The formal routine col specifies the matrix B to be factorized. To
107*  obtain j-th column of A the routine lpf_factorize calls the routine
108*  col with the parameter j (1 <= j <= n). In response the routine col
109*  should store row indices and numerical values of non-zero elements
110*  of j-th column of B to locations ind[1,...,len] and val[1,...,len],
111*  respectively, where len is the number of non-zeros in j-th column
112*  returned on exit. Neither zero nor duplicate elements are allowed.
113*
114*  The parameter info is a transit pointer passed to the routine col.
115*
116*  RETURNS
117*
118*  0  The factorization has been successfully computed.
119*
120*  LPF_ESING
121*     The specified matrix is singular within the working precision.
122*
123*  LPF_ECOND
124*     The specified matrix is ill-conditioned.
125*
126*  For more details see comments to the routine luf_factorize. */
127
128int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
129      (void *info, int j, int ind[], double val[]), void *info)
130{     int k, ret;
131#if _GLPLPF_DEBUG
132      int i, j, len, *ind;
133      double *B, *val;
134#endif
135      xassert(bh == bh);
136      if (m < 1)
137         xfault("lpf_factorize: m = %d; invalid parameter\n", m);
138      if (m > M_MAX)
139         xfault("lpf_factorize: m = %d; matrix too big\n", m);
140      lpf->m0 = lpf->m = m;
141      /* invalidate the factorization */
142      lpf->valid = 0;
143      /* allocate/reallocate arrays, if necessary */
144      if (lpf->R_ptr == NULL)
145         lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int));
146      if (lpf->R_len == NULL)
147         lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int));
148      if (lpf->S_ptr == NULL)
149         lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int));
150      if (lpf->S_len == NULL)
151         lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int));
152      if (lpf->scf == NULL)
153         lpf->scf = scf_create_it(lpf->n_max);
154      if (lpf->v_ind == NULL)
155         lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int));
156      if (lpf->v_val == NULL)
157         lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double));
158      if (lpf->m0_max < m)
159      {  if (lpf->P_row != NULL) xfree(lpf->P_row);
160         if (lpf->P_col != NULL) xfree(lpf->P_col);
161         if (lpf->Q_row != NULL) xfree(lpf->Q_row);
162         if (lpf->Q_col != NULL) xfree(lpf->Q_col);
163         if (lpf->work1 != NULL) xfree(lpf->work1);
164         if (lpf->work2 != NULL) xfree(lpf->work2);
165         lpf->m0_max = m + 100;
166         lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
167         lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
168         lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
169         lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
170         lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
171         lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
172      }
173      /* try to factorize the basis matrix */
174      switch (luf_factorize(lpf->luf, m, col, info))
175      {  case 0:
176            break;
177         case LUF_ESING:
178            ret = LPF_ESING;
179            goto done;
180         case LUF_ECOND:
181            ret = LPF_ECOND;
182            goto done;
183         default:
184            xassert(lpf != lpf);
185      }
186      /* the basis matrix has been successfully factorized */
187      lpf->valid = 1;
188#if _GLPLPF_DEBUG
189      /* store the basis matrix for debugging */
190      if (lpf->B != NULL) xfree(lpf->B);
191      xassert(m <= 32767);
192      lpf->B = B = xcalloc(1+m*m, sizeof(double));
193      ind = xcalloc(1+m, sizeof(int));
194      val = xcalloc(1+m, sizeof(double));
195      for (k = 1; k <= m * m; k++)
196         B[k] = 0.0;
197      for (j = 1; j <= m; j++)
198      {  len = col(info, j, ind, val);
199         xassert(0 <= len && len <= m);
200         for (k = 1; k <= len; k++)
201         {  i = ind[k];
202            xassert(1 <= i && i <= m);
203            xassert(B[(i - 1) * m + j] == 0.0);
204            xassert(val[k] != 0.0);
205            B[(i - 1) * m + j] = val[k];
206         }
207      }
208      xfree(ind);
209      xfree(val);
210#endif
211      /* B = B0, so there are no additional rows/columns */
212      lpf->n = 0;
213      /* reset the Schur complement factorization */
214      scf_reset_it(lpf->scf);
215      /* P := Q := I */
216      for (k = 1; k <= m; k++)
217      {  lpf->P_row[k] = lpf->P_col[k] = k;
218         lpf->Q_row[k] = lpf->Q_col[k] = k;
219      }
220      /* make all SVA locations free */
221      lpf->v_ptr = 1;
222      ret = 0;
223done: /* return to the calling program */
224      return ret;
225}
226
227/***********************************************************************
228*  The routine r_prod computes the product y := y + alpha * R * x,
229*  where x is a n-vector, alpha is a scalar, y is a m0-vector.
230*
231*  Since matrix R is available by columns, the product is computed as
232*  a linear combination:
233*
234*     y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]),
235*
236*  where R[j] is j-th column of R. */
237
238static void r_prod(LPF *lpf, double y[], double a, const double x[])
239{     int n = lpf->n;
240      int *R_ptr = lpf->R_ptr;
241      int *R_len = lpf->R_len;
242      int *v_ind = lpf->v_ind;
243      double *v_val = lpf->v_val;
244      int j, beg, end, ptr;
245      double t;
246      for (j = 1; j <= n; j++)
247      {  if (x[j] == 0.0) continue;
248         /* y := y + alpha * R[j] * x[j] */
249         t = a * x[j];
250         beg = R_ptr[j];
251         end = beg + R_len[j];
252         for (ptr = beg; ptr < end; ptr++)
253            y[v_ind[ptr]] += t * v_val[ptr];
254      }
255      return;
256}
257
258/***********************************************************************
259*  The routine rt_prod computes the product y := y + alpha * R' * x,
260*  where R' is a matrix transposed to R, x is a m0-vector, alpha is a
261*  scalar, y is a n-vector.
262*
263*  Since matrix R is available by columns, the product components are
264*  computed as inner products:
265*
266*     y[j] := y[j] + alpha * (j-th column of R) * x
267*
268*  for j = 1, 2, ..., n. */
269
270static void rt_prod(LPF *lpf, double y[], double a, const double x[])
271{     int n = lpf->n;
272      int *R_ptr = lpf->R_ptr;
273      int *R_len = lpf->R_len;
274      int *v_ind = lpf->v_ind;
275      double *v_val = lpf->v_val;
276      int j, beg, end, ptr;
277      double t;
278      for (j = 1; j <= n; j++)
279      {  /* t := (j-th column of R) * x */
280         t = 0.0;
281         beg = R_ptr[j];
282         end = beg + R_len[j];
283         for (ptr = beg; ptr < end; ptr++)
284            t += v_val[ptr] * x[v_ind[ptr]];
285         /* y[j] := y[j] + alpha * t */
286         y[j] += a * t;
287      }
288      return;
289}
290
291/***********************************************************************
292*  The routine s_prod computes the product y := y + alpha * S * x,
293*  where x is a m0-vector, alpha is a scalar, y is a n-vector.
294*
295*  Since matrix S is available by rows, the product components are
296*  computed as inner products:
297*
298*     y[i] = y[i] + alpha * (i-th row of S) * x
299*
300*  for i = 1, 2, ..., n. */
301
302static void s_prod(LPF *lpf, double y[], double a, const double x[])
303{     int n = lpf->n;
304      int *S_ptr = lpf->S_ptr;
305      int *S_len = lpf->S_len;
306      int *v_ind = lpf->v_ind;
307      double *v_val = lpf->v_val;
308      int i, beg, end, ptr;
309      double t;
310      for (i = 1; i <= n; i++)
311      {  /* t := (i-th row of S) * x */
312         t = 0.0;
313         beg = S_ptr[i];
314         end = beg + S_len[i];
315         for (ptr = beg; ptr < end; ptr++)
316            t += v_val[ptr] * x[v_ind[ptr]];
317         /* y[i] := y[i] + alpha * t */
318         y[i] += a * t;
319      }
320      return;
321}
322
323/***********************************************************************
324*  The routine st_prod computes the product y := y + alpha * S' * x,
325*  where S' is a matrix transposed to S, x is a n-vector, alpha is a
326*  scalar, y is m0-vector.
327*
328*  Since matrix R is available by rows, the product is computed as a
329*  linear combination:
330*
331*     y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]),
332*
333*  where S'[i] is i-th row of S. */
334
335static void st_prod(LPF *lpf, double y[], double a, const double x[])
336{     int n = lpf->n;
337      int *S_ptr = lpf->S_ptr;
338      int *S_len = lpf->S_len;
339      int *v_ind = lpf->v_ind;
340      double *v_val = lpf->v_val;
341      int i, beg, end, ptr;
342      double t;
343      for (i = 1; i <= n; i++)
344      {  if (x[i] == 0.0) continue;
345         /* y := y + alpha * S'[i] * x[i] */
346         t = a * x[i];
347         beg = S_ptr[i];
348         end = beg + S_len[i];
349         for (ptr = beg; ptr < end; ptr++)
350            y[v_ind[ptr]] += t * v_val[ptr];
351      }
352      return;
353}
354
355#if _GLPLPF_DEBUG
356/***********************************************************************
357*  The routine check_error computes the maximal relative error between
358*  left- and right-hand sides for the system B * x = b (if tr is zero)
359*  or B' * x = b (if tr is non-zero), where B' is a matrix transposed
360*  to B. (This routine is intended for debugging only.) */
361
362static void check_error(LPF *lpf, int tr, const double x[],
363      const double b[])
364{     int m = lpf->m;
365      double *B = lpf->B;
366      int i, j;
367      double  d, dmax = 0.0, s, t, tmax;
368      for (i = 1; i <= m; i++)
369      {  s = 0.0;
370         tmax = 1.0;
371         for (j = 1; j <= m; j++)
372         {  if (!tr)
373               t = B[m * (i - 1) + j] * x[j];
374            else
375               t = B[m * (j - 1) + i] * x[j];
376            if (tmax < fabs(t)) tmax = fabs(t);
377            s += t;
378         }
379         d = fabs(s - b[i]) / tmax;
380         if (dmax < d) dmax = d;
381      }
382      if (dmax > 1e-8)
383         xprintf("%s: dmax = %g; relative error too large\n",
384            !tr ? "lpf_ftran" : "lpf_btran", dmax);
385      return;
386}
387#endif
388
389/***********************************************************************
390*  NAME
391*
392*  lpf_ftran - perform forward transformation (solve system B*x = b)
393*
394*  SYNOPSIS
395*
396*  #include "glplpf.h"
397*  void lpf_ftran(LPF *lpf, double x[]);
398*
399*  DESCRIPTION
400*
401*  The routine lpf_ftran performs forward transformation, i.e. solves
402*  the system B*x = b, where B is the basis matrix, x is the vector of
403*  unknowns to be computed, b is the vector of right-hand sides.
404*
405*  On entry elements of the vector b should be stored in dense format
406*  in locations x[1], ..., x[m], where m is the number of rows. On exit
407*  the routine stores elements of the vector x in the same locations.
408*
409*  BACKGROUND
410*
411*  Solution of the system B * x = b can be obtained by solving the
412*  following augmented system:
413*
414*     ( B  F^) ( x )   ( b )
415*     (      ) (   ) = (   )
416*     ( G^ H^) ( y )   ( 0 )
417*
418*  which, using the main equality, can be written as follows:
419*
420*       ( L0 0 ) ( U0 R )   ( x )   ( b )
421*     P (      ) (      ) Q (   ) = (   )
422*       ( S  I ) ( 0  C )   ( y )   ( 0 )
423*
424*  therefore,
425*
426*     ( x )      ( U0 R )-1 ( L0 0 )-1    ( b )
427*     (   ) = Q' (      )   (      )   P' (   )
428*     ( y )      ( 0  C )   ( S  I )      ( 0 )
429*
430*  Thus, computing the solution includes the following steps:
431*
432*  1. Compute
433*
434*     ( f )      ( b )
435*     (   ) = P' (   )
436*     ( g )      ( 0 )
437*
438*  2. Solve the system
439*
440*     ( f1 )   ( L0 0 )-1 ( f )      ( L0 0 ) ( f1 )   ( f )
441*     (    ) = (      )   (   )  =>  (      ) (    ) = (   )
442*     ( g1 )   ( S  I )   ( g )      ( S  I ) ( g1 )   ( g )
443*
444*     from which it follows that:
445*
446*     { L0 * f1      = f      f1 = inv(L0) * f
447*     {                   =>
448*     {  S * f1 + g1 = g      g1 = g - S * f1
449*
450*  3. Solve the system
451*
452*     ( f2 )   ( U0 R )-1 ( f1 )      ( U0 R ) ( f2 )   ( f1 )
453*     (    ) = (      )   (    )  =>  (      ) (    ) = (    )
454*     ( g2 )   ( 0  C )   ( g1 )      ( 0  C ) ( g2 )   ( g1 )
455*
456*     from which it follows that:
457*
458*     { U0 * f2 + R * g2 = f1      f2 = inv(U0) * (f1 - R * g2)
459*     {                        =>
460*     {           C * g2 = g1      g2 = inv(C) * g1
461*
462*  4. Compute
463*
464*     ( x )      ( f2 )
465*     (   ) = Q' (    )
466*     ( y )      ( g2 )                                               */
467
468void lpf_ftran(LPF *lpf, double x[])
469{     int m0 = lpf->m0;
470      int m = lpf->m;
471      int n  = lpf->n;
472      int *P_col = lpf->P_col;
473      int *Q_col = lpf->Q_col;
474      double *fg = lpf->work1;
475      double *f = fg;
476      double *g = fg + m0;
477      int i, ii;
478#if _GLPLPF_DEBUG
479      double *b;
480#endif
481      if (!lpf->valid)
482         xfault("lpf_ftran: the factorization is not valid\n");
483      xassert(0 <= m && m <= m0 + n);
484#if _GLPLPF_DEBUG
485      /* save the right-hand side vector */
486      b = xcalloc(1+m, sizeof(double));
487      for (i = 1; i <= m; i++) b[i] = x[i];
488#endif
489      /* (f g) := inv(P) * (b 0) */
490      for (i = 1; i <= m0 + n; i++)
491         fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0);
492      /* f1 := inv(L0) * f */
493      luf_f_solve(lpf->luf, 0, f);
494      /* g1 := g - S * f1 */
495      s_prod(lpf, g, -1.0, f);
496      /* g2 := inv(C) * g1 */
497      scf_solve_it(lpf->scf, 0, g);
498      /* f2 := inv(U0) * (f1 - R * g2) */
499      r_prod(lpf, f, -1.0, g);
500      luf_v_solve(lpf->luf, 0, f);
501      /* (x y) := inv(Q) * (f2 g2) */
502      for (i = 1; i <= m; i++)
503         x[i] = fg[Q_col[i]];
504#if _GLPLPF_DEBUG
505      /* check relative error in solution */
506      check_error(lpf, 0, x, b);
507      xfree(b);
508#endif
509      return;
510}
511
512/***********************************************************************
513*  NAME
514*
515*  lpf_btran - perform backward transformation (solve system B'*x = b)
516*
517*  SYNOPSIS
518*
519*  #include "glplpf.h"
520*  void lpf_btran(LPF *lpf, double x[]);
521*
522*  DESCRIPTION
523*
524*  The routine lpf_btran performs backward transformation, i.e. solves
525*  the system B'*x = b, where B' is a matrix transposed to the basis
526*  matrix B, x is the vector of unknowns to be computed, b is the vector
527*  of right-hand sides.
528*
529*  On entry elements of the vector b should be stored in dense format
530*  in locations x[1], ..., x[m], where m is the number of rows. On exit
531*  the routine stores elements of the vector x in the same locations.
532*
533*  BACKGROUND
534*
535*  Solution of the system B' * x = b, where B' is a matrix transposed
536*  to B, can be obtained by solving the following augmented system:
537*
538*     ( B  F^)T ( x )   ( b )
539*     (      )  (   ) = (   )
540*     ( G^ H^)  ( y )   ( 0 )
541*
542*  which, using the main equality, can be written as follows:
543*
544*      T ( U0 R )T ( L0 0 )T  T ( x )   ( b )
545*     Q  (      )  (      )  P  (   ) = (   )
546*        ( 0  C )  ( S  I )     ( y )   ( 0 )
547*
548*  or, equivalently, as follows:
549*
550*        ( U'0 0 ) ( L'0 S')    ( x )   ( b )
551*     Q' (       ) (       ) P' (   ) = (   )
552*        ( R'  C') ( 0   I )    ( y )   ( 0 )
553*
554*  therefore,
555*
556*     ( x )     ( L'0 S')-1 ( U'0 0 )-1   ( b )
557*     (   ) = P (       )   (       )   Q (   )
558*     ( y )     ( 0   I )   ( R'  C')     ( 0 )
559*
560*  Thus, computing the solution includes the following steps:
561*
562*  1. Compute
563*
564*     ( f )     ( b )
565*     (   ) = Q (   )
566*     ( g )     ( 0 )
567*
568*  2. Solve the system
569*
570*     ( f1 )   ( U'0 0 )-1 ( f )      ( U'0 0 ) ( f1 )   ( f )
571*     (    ) = (       )   (   )  =>  (       ) (    ) = (   )
572*     ( g1 )   ( R'  C')   ( g )      ( R'  C') ( g1 )   ( g )
573*
574*     from which it follows that:
575*
576*     { U'0 * f1           = f      f1 = inv(U'0) * f
577*     {                         =>
578*     { R'  * f1 + C' * g1 = g      g1 = inv(C') * (g - R' * f1)
579*
580*  3. Solve the system
581*
582*     ( f2 )   ( L'0 S')-1 ( f1 )      ( L'0 S') ( f2 )   ( f1 )
583*     (    ) = (       )   (    )  =>  (       ) (    ) = (    )
584*     ( g2 )   ( 0   I )   ( g1 )      ( 0   I ) ( g2 )   ( g1 )
585*
586*     from which it follows that:
587*
588*     { L'0 * f2 + S' * g2 = f1
589*     {                          =>  f2 = inv(L'0) * ( f1 - S' * g2)
590*     {                 g2 = g1
591*
592*  4. Compute
593*
594*     ( x )     ( f2 )
595*     (   ) = P (    )
596*     ( y )     ( g2 )                                                */
597
598void lpf_btran(LPF *lpf, double x[])
599{     int m0 = lpf->m0;
600      int m = lpf->m;
601      int n = lpf->n;
602      int *P_row = lpf->P_row;
603      int *Q_row = lpf->Q_row;
604      double *fg = lpf->work1;
605      double *f = fg;
606      double *g = fg + m0;
607      int i, ii;
608#if _GLPLPF_DEBUG
609      double *b;
610#endif
611      if (!lpf->valid)
612         xfault("lpf_btran: the factorization is not valid\n");
613      xassert(0 <= m && m <= m0 + n);
614#if _GLPLPF_DEBUG
615      /* save the right-hand side vector */
616      b = xcalloc(1+m, sizeof(double));
617      for (i = 1; i <= m; i++) b[i] = x[i];
618#endif
619      /* (f g) := Q * (b 0) */
620      for (i = 1; i <= m0 + n; i++)
621         fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0);
622      /* f1 := inv(U'0) * f */
623      luf_v_solve(lpf->luf, 1, f);
624      /* g1 := inv(C') * (g - R' * f1) */
625      rt_prod(lpf, g, -1.0, f);
626      scf_solve_it(lpf->scf, 1, g);
627      /* g2 := g1 */
628      g = g;
629      /* f2 := inv(L'0) * (f1 - S' * g2) */
630      st_prod(lpf, f, -1.0, g);
631      luf_f_solve(lpf->luf, 1, f);
632      /* (x y) := P * (f2 g2) */
633      for (i = 1; i <= m; i++)
634         x[i] = fg[P_row[i]];
635#if _GLPLPF_DEBUG
636      /* check relative error in solution */
637      check_error(lpf, 1, x, b);
638      xfree(b);
639#endif
640      return;
641}
642
643/***********************************************************************
644*  The routine enlarge_sva enlarges the Sparse Vector Area to new_size
645*  locations by reallocating the arrays v_ind and v_val. */
646
647static void enlarge_sva(LPF *lpf, int new_size)
648{     int v_size = lpf->v_size;
649      int used = lpf->v_ptr - 1;
650      int *v_ind = lpf->v_ind;
651      double *v_val = lpf->v_val;
652      xassert(v_size < new_size);
653      while (v_size < new_size) v_size += v_size;
654      lpf->v_size = v_size;
655      lpf->v_ind = xcalloc(1+v_size, sizeof(int));
656      lpf->v_val = xcalloc(1+v_size, sizeof(double));
657      xassert(used >= 0);
658      memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int));
659      memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double));
660      xfree(v_ind);
661      xfree(v_val);
662      return;
663}
664
665/***********************************************************************
666*  NAME
667*
668*  lpf_update_it - update LP basis factorization
669*
670*  SYNOPSIS
671*
672*  #include "glplpf.h"
673*  int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
674*     const double val[]);
675*
676*  DESCRIPTION
677*
678*  The routine lpf_update_it updates the factorization of the basis
679*  matrix B after replacing its j-th column by a new vector.
680*
681*  The parameter j specifies the number of column of B, which has been
682*  replaced, 1 <= j <= m, where m is the order of B.
683*
684*  The parameter bh specifies the basis header entry for the new column
685*  of B, which is the number of the new column in some original matrix.
686*  This parameter is optional and can be specified as 0.
687*
688*  Row indices and numerical values of non-zero elements of the new
689*  column of B should be placed in locations ind[1], ..., ind[len] and
690*  val[1], ..., val[len], resp., where len is the number of non-zeros
691*  in the column. Neither zero nor duplicate elements are allowed.
692*
693*  RETURNS
694*
695*  0  The factorization has been successfully updated.
696*
697*  LPF_ESING
698*     New basis B is singular within the working precision.
699*
700*  LPF_ELIMIT
701*     Maximal number of additional rows and columns has been reached.
702*
703*  BACKGROUND
704*
705*  Let j-th column of the current basis matrix B have to be replaced by
706*  a new column a. This replacement is equivalent to removing the old
707*  j-th column by fixing it at zero and introducing the new column as
708*  follows:
709*
710*                   ( B   F^| a )
711*     ( B  F^)      (       |   )
712*     (      ) ---> ( G^  H^| 0 )
713*     ( G^ H^)      (-------+---)
714*                   ( e'j 0 | 0 )
715*
716*  where ej is a unit vector with 1 in j-th position which used to fix
717*  the old j-th column of B (at zero). Then using the main equality we
718*  have:
719*
720*     ( B   F^| a )            ( B0  F | f )
721*     (       |   )   ( P  0 ) (       |   ) ( Q  0 )
722*     ( G^  H^| 0 ) = (      ) ( G   H | g ) (      ) =
723*     (-------+---)   ( 0  1 ) (-------+---) ( 0  1 )
724*     ( e'j 0 | 0 )            ( v'  w'| 0 )
725*
726*       [   ( B0  F )|   ( f ) ]            [   ( B0 F )  |   ( f ) ]
727*       [ P (       )| P (   ) ] ( Q  0 )   [ P (      ) Q| P (   ) ]
728*     = [   ( G   H )|   ( g ) ] (      ) = [   ( G  H )  |   ( g ) ]
729*       [------------+-------- ] ( 0  1 )   [-------------+---------]
730*       [   ( v'  w')|     0   ]            [   ( v' w') Q|    0    ]
731*
732*  where:
733*
734*     ( a )     ( f )      ( f )        ( a )
735*     (   ) = P (   )  =>  (   ) = P' * (   )
736*     ( 0 )     ( g )      ( g )        ( 0 )
737*
738*                                 ( ej )      ( v )    ( v )     ( ej )
739*     ( e'j  0 ) = ( v' w' ) Q => (    ) = Q' (   ) => (   ) = Q (    )
740*                                 ( 0  )      ( w )    ( w )     ( 0  )
741*
742*  On the other hand:
743*
744*              ( B0| F  f )
745*     ( P  0 ) (---+------) ( Q  0 )         ( B0    new F )
746*     (      ) ( G | H  g ) (      ) = new P (             ) new Q
747*     ( 0  1 ) (   |      ) ( 0  1 )         ( new G new H )
748*              ( v'| w' 0 )
749*
750*  where:
751*                               ( G )           ( H  g )
752*     new F = ( F  f ), new G = (   ),  new H = (      ),
753*                               ( v')           ( w' 0 )
754*
755*             ( P  0 )           ( Q  0 )
756*     new P = (      ) , new Q = (      ) .
757*             ( 0  1 )           ( 0  1 )
758*
759*  The factorization structure for the new augmented matrix remains the
760*  same, therefore:
761*
762*           ( B0    new F )         ( L0     0 ) ( U0 new R )
763*     new P (             ) new Q = (          ) (          )
764*           ( new G new H )         ( new S  I ) ( 0  new C )
765*
766*  where:
767*
768*     new F = L0 * new R  =>
769*
770*     new R = inv(L0) * new F = inv(L0) * (F  f) = ( R  inv(L0)*f )
771*
772*     new G = new S * U0  =>
773*
774*                               ( G )             (     S      )
775*     new S = new G * inv(U0) = (   ) * inv(U0) = (            )
776*                               ( v')             ( v'*inv(U0) )
777*
778*     new H = new S * new R + new C  =>
779*
780*     new C = new H - new S * new R =
781*
782*             ( H  g )   (     S      )
783*           = (      ) - (            ) * ( R  inv(L0)*f ) =
784*             ( w' 0 )   ( v'*inv(U0) )
785*
786*             ( H - S*R           g - S*inv(L0)*f      )   ( C  x )
787*           = (                                        ) = (      )
788*             ( w'- v'*inv(U0)*R  -v'*inv(U0)*inv(L0)*f)   ( y' z )
789*
790*  Note that new C is resulted by expanding old C with new column x,
791*  row y', and diagonal element z, where:
792*
793*     x = g - S * inv(L0) * f = g - S * (new column of R)
794*
795*     y = w - R'* inv(U'0)* v = w - R'* (new row of S)
796*
797*     z = - (new row of S) * (new column of R)
798*
799*  Finally, to replace old B by new B we have to permute j-th and last
800*  (just added) columns of the matrix
801*
802*     ( B   F^| a )
803*     (       |   )
804*     ( G^  H^| 0 )
805*     (-------+---)
806*     ( e'j 0 | 0 )
807*
808*  and to keep the main equality do the same for matrix Q. */
809
810int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
811      const double val[])
812{     int m0 = lpf->m0;
813      int m = lpf->m;
814#if _GLPLPF_DEBUG
815      double *B = lpf->B;
816#endif
817      int n = lpf->n;
818      int *R_ptr = lpf->R_ptr;
819      int *R_len = lpf->R_len;
820      int *S_ptr = lpf->S_ptr;
821      int *S_len = lpf->S_len;
822      int *P_row = lpf->P_row;
823      int *P_col = lpf->P_col;
824      int *Q_row = lpf->Q_row;
825      int *Q_col = lpf->Q_col;
826      int v_ptr = lpf->v_ptr;
827      int *v_ind = lpf->v_ind;
828      double *v_val = lpf->v_val;
829      double *a = lpf->work2; /* new column */
830      double *fg = lpf->work1, *f = fg, *g = fg + m0;
831      double *vw = lpf->work2, *v = vw, *w = vw + m0;
832      double *x = g, *y = w, z;
833      int i, ii, k, ret;
834      xassert(bh == bh);
835      if (!lpf->valid)
836         xfault("lpf_update_it: the factorization is not valid\n");
837      if (!(1 <= j && j <= m))
838         xfault("lpf_update_it: j = %d; column number out of range\n",
839            j);
840      xassert(0 <= m && m <= m0 + n);
841      /* check if the basis factorization can be expanded */
842      if (n == lpf->n_max)
843      {  lpf->valid = 0;
844         ret = LPF_ELIMIT;
845         goto done;
846      }
847      /* convert new j-th column of B to dense format */
848      for (i = 1; i <= m; i++)
849         a[i] = 0.0;
850      for (k = 1; k <= len; k++)
851      {  i = ind[k];
852         if (!(1 <= i && i <= m))
853            xfault("lpf_update_it: ind[%d] = %d; row number out of rang"
854               "e\n", k, i);
855         if (a[i] != 0.0)
856            xfault("lpf_update_it: ind[%d] = %d; duplicate row index no"
857               "t allowed\n", k, i);
858         if (val[k] == 0.0)
859            xfault("lpf_update_it: val[%d] = %g; zero element not allow"
860               "ed\n", k, val[k]);
861         a[i] = val[k];
862      }
863#if _GLPLPF_DEBUG
864      /* change column in the basis matrix for debugging */
865      for (i = 1; i <= m; i++)
866         B[(i - 1) * m + j] = a[i];
867#endif
868      /* (f g) := inv(P) * (a 0) */
869      for (i = 1; i <= m0+n; i++)
870         fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0);
871      /* (v w) := Q * (ej 0) */
872      for (i = 1; i <= m0+n; i++) vw[i] = 0.0;
873      vw[Q_col[j]] = 1.0;
874      /* f1 := inv(L0) * f (new column of R) */
875      luf_f_solve(lpf->luf, 0, f);
876      /* v1 := inv(U'0) * v (new row of S) */
877      luf_v_solve(lpf->luf, 1, v);
878      /* we need at most 2 * m0 available locations in the SVA to store
879         new column of matrix R and new row of matrix S */
880      if (lpf->v_size < v_ptr + m0 + m0)
881      {  enlarge_sva(lpf, v_ptr + m0 + m0);
882         v_ind = lpf->v_ind;
883         v_val = lpf->v_val;
884      }
885      /* store new column of R */
886      R_ptr[n+1] = v_ptr;
887      for (i = 1; i <= m0; i++)
888      {  if (f[i] != 0.0)
889            v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++;
890      }
891      R_len[n+1] = v_ptr - lpf->v_ptr;
892      lpf->v_ptr = v_ptr;
893      /* store new row of S */
894      S_ptr[n+1] = v_ptr;
895      for (i = 1; i <= m0; i++)
896      {  if (v[i] != 0.0)
897            v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++;
898      }
899      S_len[n+1] = v_ptr - lpf->v_ptr;
900      lpf->v_ptr = v_ptr;
901      /* x := g - S * f1 (new column of C) */
902      s_prod(lpf, x, -1.0, f);
903      /* y := w - R' * v1 (new row of C) */
904      rt_prod(lpf, y, -1.0, v);
905      /* z := - v1 * f1 (new diagonal element of C) */
906      z = 0.0;
907      for (i = 1; i <= m0; i++) z -= v[i] * f[i];
908      /* update factorization of new matrix C */
909      switch (scf_update_exp(lpf->scf, x, y, z))
910      {  case 0:
911            break;
912         case SCF_ESING:
913            lpf->valid = 0;
914            ret = LPF_ESING;
915            goto done;
916         case SCF_ELIMIT:
917            xassert(lpf != lpf);
918         default:
919            xassert(lpf != lpf);
920      }
921      /* expand matrix P */
922      P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1;
923      /* expand matrix Q */
924      Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1;
925      /* permute j-th and last (just added) column of matrix Q */
926      i = Q_col[j], ii = Q_col[m0+n+1];
927      Q_row[i] = m0+n+1, Q_col[m0+n+1] = i;
928      Q_row[ii] = j, Q_col[j] = ii;
929      /* increase the number of additional rows and columns */
930      lpf->n++;
931      xassert(lpf->n <= lpf->n_max);
932      /* the factorization has been successfully updated */
933      ret = 0;
934done: /* return to the calling program */
935      return ret;
936}
937
938/***********************************************************************
939*  NAME
940*
941*  lpf_delete_it - delete LP basis factorization
942*
943*  SYNOPSIS
944*
945*  #include "glplpf.h"
946*  void lpf_delete_it(LPF *lpf)
947*
948*  DESCRIPTION
949*
950*  The routine lpf_delete_it deletes LP basis factorization specified
951*  by the parameter lpf and frees all memory allocated to this program
952*  object. */
953
954void lpf_delete_it(LPF *lpf)
955{     luf_delete_it(lpf->luf);
956#if _GLPLPF_DEBUG
957      if (lpf->B != NULL) xfree(lpf->B);
958#else
959      xassert(lpf->B == NULL);
960#endif
961      if (lpf->R_ptr != NULL) xfree(lpf->R_ptr);
962      if (lpf->R_len != NULL) xfree(lpf->R_len);
963      if (lpf->S_ptr != NULL) xfree(lpf->S_ptr);
964      if (lpf->S_len != NULL) xfree(lpf->S_len);
965      if (lpf->scf != NULL) scf_delete_it(lpf->scf);
966      if (lpf->P_row != NULL) xfree(lpf->P_row);
967      if (lpf->P_col != NULL) xfree(lpf->P_col);
968      if (lpf->Q_row != NULL) xfree(lpf->Q_row);
969      if (lpf->Q_col != NULL) xfree(lpf->Q_col);
970      if (lpf->v_ind != NULL) xfree(lpf->v_ind);
971      if (lpf->v_val != NULL) xfree(lpf->v_val);
972      if (lpf->work1 != NULL) xfree(lpf->work1);
973      if (lpf->work2 != NULL) xfree(lpf->work2);
974      xfree(lpf);
975      return;
976}
977
978/* eof */
Note: See TracBrowser for help on using the repository browser.