COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpfhv.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: 26.8 KB
Line 
1/* glpfhv.c (LP basis factorization, FHV eta file 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 "glpfhv.h"
26#include "glpenv.h"
27#define xfault xerror
28
29/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
30
31#define M_MAX 100000000 /* = 100*10^6 */
32/* maximal order of the basis matrix */
33
34/***********************************************************************
35*  NAME
36*
37*  fhv_create_it - create LP basis factorization
38*
39*  SYNOPSIS
40*
41*  #include "glpfhv.h"
42*  FHV *fhv_create_it(void);
43*
44*  DESCRIPTION
45*
46*  The routine fhv_create_it creates a program object, which represents
47*  a factorization of LP basis.
48*
49*  RETURNS
50*
51*  The routine fhv_create_it returns a pointer to the object created. */
52
53FHV *fhv_create_it(void)
54{     FHV *fhv;
55      fhv = xmalloc(sizeof(FHV));
56      fhv->m_max = fhv->m = 0;
57      fhv->valid = 0;
58      fhv->luf = luf_create_it();
59      fhv->hh_max = 50;
60      fhv->hh_nfs = 0;
61      fhv->hh_ind = fhv->hh_ptr = fhv->hh_len = NULL;
62      fhv->p0_row = fhv->p0_col = NULL;
63      fhv->cc_ind = NULL;
64      fhv->cc_val = NULL;
65      fhv->upd_tol = 1e-6;
66      fhv->nnz_h = 0;
67      return fhv;
68}
69
70/***********************************************************************
71*  NAME
72*
73*  fhv_factorize - compute LP basis factorization
74*
75*  SYNOPSIS
76*
77*  #include "glpfhv.h"
78*  int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
79*     int ind[], double val[]), void *info);
80*
81*  DESCRIPTION
82*
83*  The routine fhv_factorize computes the factorization of the basis
84*  matrix B specified by the routine col.
85*
86*  The parameter fhv specified the basis factorization data structure
87*  created by the routine fhv_create_it.
88*
89*  The parameter m specifies the order of B, m > 0.
90*
91*  The formal routine col specifies the matrix B to be factorized. To
92*  obtain j-th column of A the routine fhv_factorize calls the routine
93*  col with the parameter j (1 <= j <= n). In response the routine col
94*  should store row indices and numerical values of non-zero elements
95*  of j-th column of B to locations ind[1,...,len] and val[1,...,len],
96*  respectively, where len is the number of non-zeros in j-th column
97*  returned on exit. Neither zero nor duplicate elements are allowed.
98*
99*  The parameter info is a transit pointer passed to the routine col.
100*
101*  RETURNS
102*
103*  0  The factorization has been successfully computed.
104*
105*  FHV_ESING
106*     The specified matrix is singular within the working precision.
107*
108*  FHV_ECOND
109*     The specified matrix is ill-conditioned.
110*
111*  For more details see comments to the routine luf_factorize.
112*
113*  ALGORITHM
114*
115*  The routine fhv_factorize calls the routine luf_factorize (see the
116*  module GLPLUF), which actually computes LU-factorization of the basis
117*  matrix B in the form
118*
119*     [B] = (F, V, P, Q),
120*
121*  where F and V are such matrices that
122*
123*     B = F * V,
124*
125*  and P and Q are such permutation matrices that the matrix
126*
127*     L = P * F * inv(P)
128*
129*  is lower triangular with unity diagonal, and the matrix
130*
131*     U = P * V * Q
132*
133*  is upper triangular.
134*
135*  In order to build the complete representation of the factorization
136*  (see formula (1) in the file glpfhv.h) the routine fhv_factorize just
137*  additionally sets H = I and P0 = P. */
138
139int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
140      int ind[], double val[]), void *info)
141{     int ret;
142      if (m < 1)
143         xfault("fhv_factorize: m = %d; invalid parameter\n", m);
144      if (m > M_MAX)
145         xfault("fhv_factorize: m = %d; matrix too big\n", m);
146      fhv->m = m;
147      /* invalidate the factorization */
148      fhv->valid = 0;
149      /* allocate/reallocate arrays, if necessary */
150      if (fhv->hh_ind == NULL)
151         fhv->hh_ind = xcalloc(1+fhv->hh_max, sizeof(int));
152      if (fhv->hh_ptr == NULL)
153         fhv->hh_ptr = xcalloc(1+fhv->hh_max, sizeof(int));
154      if (fhv->hh_len == NULL)
155         fhv->hh_len = xcalloc(1+fhv->hh_max, sizeof(int));
156      if (fhv->m_max < m)
157      {  if (fhv->p0_row != NULL) xfree(fhv->p0_row);
158         if (fhv->p0_col != NULL) xfree(fhv->p0_col);
159         if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);
160         if (fhv->cc_val != NULL) xfree(fhv->cc_val);
161         fhv->m_max = m + 100;
162         fhv->p0_row = xcalloc(1+fhv->m_max, sizeof(int));
163         fhv->p0_col = xcalloc(1+fhv->m_max, sizeof(int));
164         fhv->cc_ind = xcalloc(1+fhv->m_max, sizeof(int));
165         fhv->cc_val = xcalloc(1+fhv->m_max, sizeof(double));
166      }
167      /* try to factorize the basis matrix */
168      switch (luf_factorize(fhv->luf, m, col, info))
169      {  case 0:
170            break;
171         case LUF_ESING:
172            ret = FHV_ESING;
173            goto done;
174         case LUF_ECOND:
175            ret = FHV_ECOND;
176            goto done;
177         default:
178            xassert(fhv != fhv);
179      }
180      /* the basis matrix has been successfully factorized */
181      fhv->valid = 1;
182      /* H := I */
183      fhv->hh_nfs = 0;
184      /* P0 := P */
185      memcpy(&fhv->p0_row[1], &fhv->luf->pp_row[1], sizeof(int) * m);
186      memcpy(&fhv->p0_col[1], &fhv->luf->pp_col[1], sizeof(int) * m);
187      /* currently H has no factors */
188      fhv->nnz_h = 0;
189      ret = 0;
190done: /* return to the calling program */
191      return ret;
192}
193
194/***********************************************************************
195*  NAME
196*
197*  fhv_h_solve - solve system H*x = b or H'*x = b
198*
199*  SYNOPSIS
200*
201*  #include "glpfhv.h"
202*  void fhv_h_solve(FHV *fhv, int tr, double x[]);
203*
204*  DESCRIPTION
205*
206*  The routine fhv_h_solve solves either the system H*x = b (if the
207*  flag tr is zero) or the system H'*x = b (if the flag tr is non-zero),
208*  where the matrix H is a component of the factorization specified by
209*  the parameter fhv, H' is a matrix transposed to H.
210*
211*  On entry the array x should contain elements of the right-hand side
212*  vector b in locations x[1], ..., x[m], where m is the order of the
213*  matrix H. On exit this array will contain elements of the solution
214*  vector x in the same locations. */
215
216void fhv_h_solve(FHV *fhv, int tr, double x[])
217{     int nfs = fhv->hh_nfs;
218      int *hh_ind = fhv->hh_ind;
219      int *hh_ptr = fhv->hh_ptr;
220      int *hh_len = fhv->hh_len;
221      int *sv_ind = fhv->luf->sv_ind;
222      double *sv_val = fhv->luf->sv_val;
223      int i, k, beg, end, ptr;
224      double temp;
225      if (!fhv->valid)
226         xfault("fhv_h_solve: the factorization is not valid\n");
227      if (!tr)
228      {  /* solve the system H*x = b */
229         for (k = 1; k <= nfs; k++)
230         {  i = hh_ind[k];
231            temp = x[i];
232            beg = hh_ptr[k];
233            end = beg + hh_len[k] - 1;
234            for (ptr = beg; ptr <= end; ptr++)
235               temp -= sv_val[ptr] * x[sv_ind[ptr]];
236            x[i] = temp;
237         }
238      }
239      else
240      {  /* solve the system H'*x = b */
241         for (k = nfs; k >= 1; k--)
242         {  i = hh_ind[k];
243            temp = x[i];
244            if (temp == 0.0) continue;
245            beg = hh_ptr[k];
246            end = beg + hh_len[k] - 1;
247            for (ptr = beg; ptr <= end; ptr++)
248               x[sv_ind[ptr]] -= sv_val[ptr] * temp;
249         }
250      }
251      return;
252}
253
254/***********************************************************************
255*  NAME
256*
257*  fhv_ftran - perform forward transformation (solve system B*x = b)
258*
259*  SYNOPSIS
260*
261*  #include "glpfhv.h"
262*  void fhv_ftran(FHV *fhv, double x[]);
263*
264*  DESCRIPTION
265*
266*  The routine fhv_ftran performs forward transformation, i.e. solves
267*  the system B*x = b, where B is the basis matrix, x is the vector of
268*  unknowns to be computed, b is the vector of right-hand sides.
269*
270*  On entry elements of the vector b should be stored in dense format
271*  in locations x[1], ..., x[m], where m is the number of rows. On exit
272*  the routine stores elements of the vector x in the same locations. */
273
274void fhv_ftran(FHV *fhv, double x[])
275{     int *pp_row = fhv->luf->pp_row;
276      int *pp_col = fhv->luf->pp_col;
277      int *p0_row = fhv->p0_row;
278      int *p0_col = fhv->p0_col;
279      if (!fhv->valid)
280         xfault("fhv_ftran: the factorization is not valid\n");
281      /* B = F*H*V, therefore inv(B) = inv(V)*inv(H)*inv(F) */
282      fhv->luf->pp_row = p0_row;
283      fhv->luf->pp_col = p0_col;
284      luf_f_solve(fhv->luf, 0, x);
285      fhv->luf->pp_row = pp_row;
286      fhv->luf->pp_col = pp_col;
287      fhv_h_solve(fhv, 0, x);
288      luf_v_solve(fhv->luf, 0, x);
289      return;
290}
291
292/***********************************************************************
293*  NAME
294*
295*  fhv_btran - perform backward transformation (solve system B'*x = b)
296*
297*  SYNOPSIS
298*
299*  #include "glpfhv.h"
300*  void fhv_btran(FHV *fhv, double x[]);
301*
302*  DESCRIPTION
303*
304*  The routine fhv_btran performs backward transformation, i.e. solves
305*  the system B'*x = b, where B' is a matrix transposed to the basis
306*  matrix B, x is the vector of unknowns to be computed, b is the vector
307*  of right-hand sides.
308*
309*  On entry elements of the vector b should be stored in dense format
310*  in locations x[1], ..., x[m], where m is the number of rows. On exit
311*  the routine stores elements of the vector x in the same locations. */
312
313void fhv_btran(FHV *fhv, double x[])
314{     int *pp_row = fhv->luf->pp_row;
315      int *pp_col = fhv->luf->pp_col;
316      int *p0_row = fhv->p0_row;
317      int *p0_col = fhv->p0_col;
318      if (!fhv->valid)
319         xfault("fhv_btran: the factorization is not valid\n");
320      /* B = F*H*V, therefore inv(B') = inv(F')*inv(H')*inv(V') */
321      luf_v_solve(fhv->luf, 1, x);
322      fhv_h_solve(fhv, 1, x);
323      fhv->luf->pp_row = p0_row;
324      fhv->luf->pp_col = p0_col;
325      luf_f_solve(fhv->luf, 1, x);
326      fhv->luf->pp_row = pp_row;
327      fhv->luf->pp_col = pp_col;
328      return;
329}
330
331/***********************************************************************
332*  NAME
333*
334*  fhv_update_it - update LP basis factorization
335*
336*  SYNOPSIS
337*
338*  #include "glpfhv.h"
339*  int fhv_update_it(FHV *fhv, int j, int len, const int ind[],
340*     const double val[]);
341*
342*  DESCRIPTION
343*
344*  The routine fhv_update_it updates the factorization of the basis
345*  matrix B after replacing its j-th column by a new vector.
346*
347*  The parameter j specifies the number of column of B, which has been
348*  replaced, 1 <= j <= m, where m is the order of B.
349*
350*  Row indices and numerical values of non-zero elements of the new
351*  column of B should be placed in locations ind[1], ..., ind[len] and
352*  val[1], ..., val[len], resp., where len is the number of non-zeros
353*  in the column. Neither zero nor duplicate elements are allowed.
354*
355*  RETURNS
356*
357*  0  The factorization has been successfully updated.
358*
359*  FHV_ESING
360*     The adjacent basis matrix is structurally singular, since after
361*     changing j-th column of matrix V by the new column (see algorithm
362*     below) the case k1 > k2 occured.
363*
364*  FHV_ECHECK
365*     The factorization is inaccurate, since after transforming k2-th
366*     row of matrix U = P*V*Q, its diagonal element u[k2,k2] is zero or
367*     close to zero,
368*
369*  FHV_ELIMIT
370*     Maximal number of H factors has been reached.
371*
372*  FHV_EROOM
373*     Overflow of the sparse vector area.
374*
375*  In case of non-zero return code the factorization becomes invalid.
376*  It should not be used until it has been recomputed with the routine
377*  fhv_factorize.
378*
379*  ALGORITHM
380*
381*  The routine fhv_update_it is based on the transformation proposed by
382*  Forrest and Tomlin.
383*
384*  Let j-th column of the basis matrix B have been replaced by new
385*  column B[j]. In order to keep the equality B = F*H*V j-th column of
386*  matrix V should be replaced by the column inv(F*H)*B[j].
387*
388*  From the standpoint of matrix U = P*V*Q, replacement of j-th column
389*  of matrix V is equivalent to replacement of k1-th column of matrix U,
390*  where k1 is determined by permutation matrix Q. Thus, matrix U loses
391*  its upper triangular form and becomes the following:
392*
393*         1   k1       k2   m
394*     1   x x * x x x x x x x
395*         . x * x x x x x x x
396*     k1  . . * x x x x x x x
397*         . . * x x x x x x x
398*         . . * . x x x x x x
399*         . . * . . x x x x x
400*         . . * . . . x x x x
401*     k2  . . * . . . . x x x
402*         . . . . . . . . x x
403*     m   . . . . . . . . . x
404*
405*  where row index k2 corresponds to the lowest non-zero element of
406*  k1-th column.
407*
408*  The routine moves rows and columns k1+1, k1+2, ..., k2 of matrix U
409*  by one position to the left and upwards and moves k1-th row and k1-th
410*  column to position k2. As the result of such symmetric permutations
411*  matrix U becomes the following:
412*
413*         1   k1       k2   m
414*     1   x x x x x x x * x x
415*         . x x x x x x * x x
416*     k1  . . x x x x x * x x
417*         . . . x x x x * x x
418*         . . . . x x x * x x
419*         . . . . . x x * x x
420*         . . . . . . x * x x
421*     k2  . . x x x x x * x x
422*         . . . . . . . . x x
423*     m   . . . . . . . . . x
424*
425*  Then the routine performs gaussian elimination to eliminate elements
426*  u[k2,k1], u[k2,k1+1], ..., u[k2,k2-1] using diagonal elements
427*  u[k1,k1], u[k1+1,k1+1], ..., u[k2-1,k2-1] as pivots in the same way
428*  as described in comments to the routine luf_factorize (see the module
429*  GLPLUF). Note that actually all operations are performed on matrix V,
430*  not on matrix U. During the elimination process the routine permutes
431*  neither rows nor columns, so only k2-th row of matrix U is changed.
432*
433*  To keep the main equality B = F*H*V, each time when the routine
434*  applies elementary gaussian transformation to the transformed row of
435*  matrix V (which corresponds to k2-th row of matrix U), it also adds
436*  a new element (gaussian multiplier) to the current row-like factor
437*  of matrix H, which corresponds to the transformed row of matrix V. */
438
439int fhv_update_it(FHV *fhv, int j, int len, const int ind[],
440      const double val[])
441{     int m = fhv->m;
442      LUF *luf = fhv->luf;
443      int *vr_ptr = luf->vr_ptr;
444      int *vr_len = luf->vr_len;
445      int *vr_cap = luf->vr_cap;
446      double *vr_piv = luf->vr_piv;
447      int *vc_ptr = luf->vc_ptr;
448      int *vc_len = luf->vc_len;
449      int *vc_cap = luf->vc_cap;
450      int *pp_row = luf->pp_row;
451      int *pp_col = luf->pp_col;
452      int *qq_row = luf->qq_row;
453      int *qq_col = luf->qq_col;
454      int *sv_ind = luf->sv_ind;
455      double *sv_val = luf->sv_val;
456      double *work = luf->work;
457      double eps_tol = luf->eps_tol;
458      int *hh_ind = fhv->hh_ind;
459      int *hh_ptr = fhv->hh_ptr;
460      int *hh_len = fhv->hh_len;
461      int *p0_row = fhv->p0_row;
462      int *p0_col = fhv->p0_col;
463      int *cc_ind = fhv->cc_ind;
464      double *cc_val = fhv->cc_val;
465      double upd_tol = fhv->upd_tol;
466      int i, i_beg, i_end, i_ptr, j_beg, j_end, j_ptr, k, k1, k2, p, q,
467         p_beg, p_end, p_ptr, ptr, ret;
468      double f, temp;
469      if (!fhv->valid)
470         xfault("fhv_update_it: the factorization is not valid\n");
471      if (!(1 <= j && j <= m))
472         xfault("fhv_update_it: j = %d; column number out of range\n",
473            j);
474      /* check if the new factor of matrix H can be created */
475      if (fhv->hh_nfs == fhv->hh_max)
476      {  /* maximal number of updates has been reached */
477         fhv->valid = 0;
478         ret = FHV_ELIMIT;
479         goto done;
480      }
481      /* convert new j-th column of B to dense format */
482      for (i = 1; i <= m; i++)
483         cc_val[i] = 0.0;
484      for (k = 1; k <= len; k++)
485      {  i = ind[k];
486         if (!(1 <= i && i <= m))
487            xfault("fhv_update_it: ind[%d] = %d; row number out of rang"
488               "e\n", k, i);
489         if (cc_val[i] != 0.0)
490            xfault("fhv_update_it: ind[%d] = %d; duplicate row index no"
491               "t allowed\n", k, i);
492         if (val[k] == 0.0)
493            xfault("fhv_update_it: val[%d] = %g; zero element not allow"
494               "ed\n", k, val[k]);
495         cc_val[i] = val[k];
496      }
497      /* new j-th column of V := inv(F * H) * (new B[j]) */
498      fhv->luf->pp_row = p0_row;
499      fhv->luf->pp_col = p0_col;
500      luf_f_solve(fhv->luf, 0, cc_val);
501      fhv->luf->pp_row = pp_row;
502      fhv->luf->pp_col = pp_col;
503      fhv_h_solve(fhv, 0, cc_val);
504      /* convert new j-th column of V to sparse format */
505      len = 0;
506      for (i = 1; i <= m; i++)
507      {  temp = cc_val[i];
508         if (temp == 0.0 || fabs(temp) < eps_tol) continue;
509         len++, cc_ind[len] = i, cc_val[len] = temp;
510      }
511      /* clear old content of j-th column of matrix V */
512      j_beg = vc_ptr[j];
513      j_end = j_beg + vc_len[j] - 1;
514      for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
515      {  /* get row index of v[i,j] */
516         i = sv_ind[j_ptr];
517         /* find v[i,j] in the i-th row */
518         i_beg = vr_ptr[i];
519         i_end = i_beg + vr_len[i] - 1;
520         for (i_ptr = i_beg; sv_ind[i_ptr] != j; i_ptr++) /* nop */;
521         xassert(i_ptr <= i_end);
522         /* remove v[i,j] from the i-th row */
523         sv_ind[i_ptr] = sv_ind[i_end];
524         sv_val[i_ptr] = sv_val[i_end];
525         vr_len[i]--;
526      }
527      /* now j-th column of matrix V is empty */
528      luf->nnz_v -= vc_len[j];
529      vc_len[j] = 0;
530      /* add new elements of j-th column of matrix V to corresponding
531         row lists; determine indices k1 and k2 */
532      k1 = qq_row[j], k2 = 0;
533      for (ptr = 1; ptr <= len; ptr++)
534      {  /* get row index of v[i,j] */
535         i = cc_ind[ptr];
536         /* at least one unused location is needed in i-th row */
537         if (vr_len[i] + 1 > vr_cap[i])
538         {  if (luf_enlarge_row(luf, i, vr_len[i] + 10))
539            {  /* overflow of the sparse vector area */
540               fhv->valid = 0;
541               luf->new_sva = luf->sv_size + luf->sv_size;
542               xassert(luf->new_sva > luf->sv_size);
543               ret = FHV_EROOM;
544               goto done;
545            }
546         }
547         /* add v[i,j] to i-th row */
548         i_ptr = vr_ptr[i] + vr_len[i];
549         sv_ind[i_ptr] = j;
550         sv_val[i_ptr] = cc_val[ptr];
551         vr_len[i]++;
552         /* adjust index k2 */
553         if (k2 < pp_col[i]) k2 = pp_col[i];
554      }
555      /* capacity of j-th column (which is currently empty) should be
556         not less than len locations */
557      if (vc_cap[j] < len)
558      {  if (luf_enlarge_col(luf, j, len))
559         {  /* overflow of the sparse vector area */
560            fhv->valid = 0;
561            luf->new_sva = luf->sv_size + luf->sv_size;
562            xassert(luf->new_sva > luf->sv_size);
563            ret = FHV_EROOM;
564            goto done;
565         }
566      }
567      /* add new elements of matrix V to j-th column list */
568      j_ptr = vc_ptr[j];
569      memmove(&sv_ind[j_ptr], &cc_ind[1], len * sizeof(int));
570      memmove(&sv_val[j_ptr], &cc_val[1], len * sizeof(double));
571      vc_len[j] = len;
572      luf->nnz_v += len;
573      /* if k1 > k2, diagonal element u[k2,k2] of matrix U is zero and
574         therefore the adjacent basis matrix is structurally singular */
575      if (k1 > k2)
576      {  fhv->valid = 0;
577         ret = FHV_ESING;
578         goto done;
579      }
580      /* perform implicit symmetric permutations of rows and columns of
581         matrix U */
582      i = pp_row[k1], j = qq_col[k1];
583      for (k = k1; k < k2; k++)
584      {  pp_row[k] = pp_row[k+1], pp_col[pp_row[k]] = k;
585         qq_col[k] = qq_col[k+1], qq_row[qq_col[k]] = k;
586      }
587      pp_row[k2] = i, pp_col[i] = k2;
588      qq_col[k2] = j, qq_row[j] = k2;
589      /* now i-th row of the matrix V is k2-th row of matrix U; since
590         no pivoting is used, only this row will be transformed */
591      /* copy elements of i-th row of matrix V to the working array and
592         remove these elements from matrix V */
593      for (j = 1; j <= m; j++) work[j] = 0.0;
594      i_beg = vr_ptr[i];
595      i_end = i_beg + vr_len[i] - 1;
596      for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
597      {  /* get column index of v[i,j] */
598         j = sv_ind[i_ptr];
599         /* store v[i,j] to the working array */
600         work[j] = sv_val[i_ptr];
601         /* find v[i,j] in the j-th column */
602         j_beg = vc_ptr[j];
603         j_end = j_beg + vc_len[j] - 1;
604         for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++) /* nop */;
605         xassert(j_ptr <= j_end);
606         /* remove v[i,j] from the j-th column */
607         sv_ind[j_ptr] = sv_ind[j_end];
608         sv_val[j_ptr] = sv_val[j_end];
609         vc_len[j]--;
610      }
611      /* now i-th row of matrix V is empty */
612      luf->nnz_v -= vr_len[i];
613      vr_len[i] = 0;
614      /* create the next row-like factor of the matrix H; this factor
615         corresponds to i-th (transformed) row */
616      fhv->hh_nfs++;
617      hh_ind[fhv->hh_nfs] = i;
618      /* hh_ptr[] will be set later */
619      hh_len[fhv->hh_nfs] = 0;
620      /* up to (k2 - k1) free locations are needed to add new elements
621         to the non-trivial row of the row-like factor */
622      if (luf->sv_end - luf->sv_beg < k2 - k1)
623      {  luf_defrag_sva(luf);
624         if (luf->sv_end - luf->sv_beg < k2 - k1)
625         {  /* overflow of the sparse vector area */
626            fhv->valid = luf->valid = 0;
627            luf->new_sva = luf->sv_size + luf->sv_size;
628            xassert(luf->new_sva > luf->sv_size);
629            ret = FHV_EROOM;
630            goto done;
631         }
632      }
633      /* eliminate subdiagonal elements of matrix U */
634      for (k = k1; k < k2; k++)
635      {  /* v[p,q] = u[k,k] */
636         p = pp_row[k], q = qq_col[k];
637         /* this is the crucial point, where even tiny non-zeros should
638            not be dropped */
639         if (work[q] == 0.0) continue;
640         /* compute gaussian multiplier f = v[i,q] / v[p,q] */
641         f = work[q] / vr_piv[p];
642         /* perform gaussian transformation:
643            (i-th row) := (i-th row) - f * (p-th row)
644            in order to eliminate v[i,q] = u[k2,k] */
645         p_beg = vr_ptr[p];
646         p_end = p_beg + vr_len[p] - 1;
647         for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
648            work[sv_ind[p_ptr]] -= f * sv_val[p_ptr];
649         /* store new element (gaussian multiplier that corresponds to
650            p-th row) in the current row-like factor */
651         luf->sv_end--;
652         sv_ind[luf->sv_end] = p;
653         sv_val[luf->sv_end] = f;
654         hh_len[fhv->hh_nfs]++;
655      }
656      /* set pointer to the current row-like factor of the matrix H
657         (if no elements were added to this factor, it is unity matrix
658         and therefore can be discarded) */
659      if (hh_len[fhv->hh_nfs] == 0)
660         fhv->hh_nfs--;
661      else
662      {  hh_ptr[fhv->hh_nfs] = luf->sv_end;
663         fhv->nnz_h += hh_len[fhv->hh_nfs];
664      }
665      /* store new pivot which corresponds to u[k2,k2] */
666      vr_piv[i] = work[qq_col[k2]];
667      /* new elements of i-th row of matrix V (which are non-diagonal
668         elements u[k2,k2+1], ..., u[k2,m] of matrix U = P*V*Q) now are
669         contained in the working array; add them to matrix V */
670      len = 0;
671      for (k = k2+1; k <= m; k++)
672      {  /* get column index and value of v[i,j] = u[k2,k] */
673         j = qq_col[k];
674         temp = work[j];
675         /* if v[i,j] is close to zero, skip it */
676         if (fabs(temp) < eps_tol) continue;
677         /* at least one unused location is needed in j-th column */
678         if (vc_len[j] + 1 > vc_cap[j])
679         {  if (luf_enlarge_col(luf, j, vc_len[j] + 10))
680            {  /* overflow of the sparse vector area */
681               fhv->valid = 0;
682               luf->new_sva = luf->sv_size + luf->sv_size;
683               xassert(luf->new_sva > luf->sv_size);
684               ret = FHV_EROOM;
685               goto done;
686            }
687         }
688         /* add v[i,j] to j-th column */
689         j_ptr = vc_ptr[j] + vc_len[j];
690         sv_ind[j_ptr] = i;
691         sv_val[j_ptr] = temp;
692         vc_len[j]++;
693         /* also store v[i,j] to the auxiliary array */
694         len++, cc_ind[len] = j, cc_val[len] = temp;
695      }
696      /* capacity of i-th row (which is currently empty) should be not
697         less than len locations */
698      if (vr_cap[i] < len)
699      {  if (luf_enlarge_row(luf, i, len))
700         {  /* overflow of the sparse vector area */
701            fhv->valid = 0;
702            luf->new_sva = luf->sv_size + luf->sv_size;
703            xassert(luf->new_sva > luf->sv_size);
704            ret = FHV_EROOM;
705            goto done;
706         }
707      }
708      /* add new elements to i-th row list */
709      i_ptr = vr_ptr[i];
710      memmove(&sv_ind[i_ptr], &cc_ind[1], len * sizeof(int));
711      memmove(&sv_val[i_ptr], &cc_val[1], len * sizeof(double));
712      vr_len[i] = len;
713      luf->nnz_v += len;
714      /* updating is finished; check that diagonal element u[k2,k2] is
715         not very small in absolute value among other elements in k2-th
716         row and k2-th column of matrix U = P*V*Q */
717      /* temp = max(|u[k2,*]|, |u[*,k2]|) */
718      temp = 0.0;
719      /* walk through k2-th row of U which is i-th row of V */
720      i = pp_row[k2];
721      i_beg = vr_ptr[i];
722      i_end = i_beg + vr_len[i] - 1;
723      for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
724         if (temp < fabs(sv_val[i_ptr])) temp = fabs(sv_val[i_ptr]);
725      /* walk through k2-th column of U which is j-th column of V */
726      j = qq_col[k2];
727      j_beg = vc_ptr[j];
728      j_end = j_beg + vc_len[j] - 1;
729      for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
730         if (temp < fabs(sv_val[j_ptr])) temp = fabs(sv_val[j_ptr]);
731      /* check that u[k2,k2] is not very small */
732      if (fabs(vr_piv[i]) < upd_tol * temp)
733      {  /* the factorization seems to be inaccurate and therefore must
734            be recomputed */
735         fhv->valid = 0;
736         ret = FHV_ECHECK;
737         goto done;
738      }
739      /* the factorization has been successfully updated */
740      ret = 0;
741done: /* return to the calling program */
742      return ret;
743}
744
745/***********************************************************************
746*  NAME
747*
748*  fhv_delete_it - delete LP basis factorization
749*
750*  SYNOPSIS
751*
752*  #include "glpfhv.h"
753*  void fhv_delete_it(FHV *fhv);
754*
755*  DESCRIPTION
756*
757*  The routine fhv_delete_it deletes LP basis factorization specified
758*  by the parameter fhv and frees all memory allocated to this program
759*  object. */
760
761void fhv_delete_it(FHV *fhv)
762{     luf_delete_it(fhv->luf);
763      if (fhv->hh_ind != NULL) xfree(fhv->hh_ind);
764      if (fhv->hh_ptr != NULL) xfree(fhv->hh_ptr);
765      if (fhv->hh_len != NULL) xfree(fhv->hh_len);
766      if (fhv->p0_row != NULL) xfree(fhv->p0_row);
767      if (fhv->p0_col != NULL) xfree(fhv->p0_col);
768      if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);
769      if (fhv->cc_val != NULL) xfree(fhv->cc_val);
770      xfree(fhv);
771      return;
772}
773
774/* eof */
Note: See TracBrowser for help on using the repository browser.