/* glpfhv.c (LP basis factorization, FHV eta file version) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, * Moscow Aviation Institute, Moscow, Russia. All rights reserved. * E-mail: . * * GLPK is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * GLPK is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * License for more details. * * You should have received a copy of the GNU General Public License * along with GLPK. If not, see . ***********************************************************************/ #include "glpfhv.h" #include "glpenv.h" #define xfault xerror /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ #define M_MAX 100000000 /* = 100*10^6 */ /* maximal order of the basis matrix */ /*********************************************************************** * NAME * * fhv_create_it - create LP basis factorization * * SYNOPSIS * * #include "glpfhv.h" * FHV *fhv_create_it(void); * * DESCRIPTION * * The routine fhv_create_it creates a program object, which represents * a factorization of LP basis. * * RETURNS * * The routine fhv_create_it returns a pointer to the object created. */ FHV *fhv_create_it(void) { FHV *fhv; fhv = xmalloc(sizeof(FHV)); fhv->m_max = fhv->m = 0; fhv->valid = 0; fhv->luf = luf_create_it(); fhv->hh_max = 50; fhv->hh_nfs = 0; fhv->hh_ind = fhv->hh_ptr = fhv->hh_len = NULL; fhv->p0_row = fhv->p0_col = NULL; fhv->cc_ind = NULL; fhv->cc_val = NULL; fhv->upd_tol = 1e-6; fhv->nnz_h = 0; return fhv; } /*********************************************************************** * NAME * * fhv_factorize - compute LP basis factorization * * SYNOPSIS * * #include "glpfhv.h" * int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j, * int ind[], double val[]), void *info); * * DESCRIPTION * * The routine fhv_factorize computes the factorization of the basis * matrix B specified by the routine col. * * The parameter fhv specified the basis factorization data structure * created by the routine fhv_create_it. * * The parameter m specifies the order of B, m > 0. * * The formal routine col specifies the matrix B to be factorized. To * obtain j-th column of A the routine fhv_factorize calls the routine * col with the parameter j (1 <= j <= n). In response the routine col * should store row indices and numerical values of non-zero elements * of j-th column of B to locations ind[1,...,len] and val[1,...,len], * respectively, where len is the number of non-zeros in j-th column * returned on exit. Neither zero nor duplicate elements are allowed. * * The parameter info is a transit pointer passed to the routine col. * * RETURNS * * 0 The factorization has been successfully computed. * * FHV_ESING * The specified matrix is singular within the working precision. * * FHV_ECOND * The specified matrix is ill-conditioned. * * For more details see comments to the routine luf_factorize. * * ALGORITHM * * The routine fhv_factorize calls the routine luf_factorize (see the * module GLPLUF), which actually computes LU-factorization of the basis * matrix B in the form * * [B] = (F, V, P, Q), * * where F and V are such matrices that * * B = F * V, * * and P and Q are such permutation matrices that the matrix * * L = P * F * inv(P) * * is lower triangular with unity diagonal, and the matrix * * U = P * V * Q * * is upper triangular. * * In order to build the complete representation of the factorization * (see formula (1) in the file glpfhv.h) the routine fhv_factorize just * additionally sets H = I and P0 = P. */ int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j, int ind[], double val[]), void *info) { int ret; if (m < 1) xfault("fhv_factorize: m = %d; invalid parameter\n", m); if (m > M_MAX) xfault("fhv_factorize: m = %d; matrix too big\n", m); fhv->m = m; /* invalidate the factorization */ fhv->valid = 0; /* allocate/reallocate arrays, if necessary */ if (fhv->hh_ind == NULL) fhv->hh_ind = xcalloc(1+fhv->hh_max, sizeof(int)); if (fhv->hh_ptr == NULL) fhv->hh_ptr = xcalloc(1+fhv->hh_max, sizeof(int)); if (fhv->hh_len == NULL) fhv->hh_len = xcalloc(1+fhv->hh_max, sizeof(int)); if (fhv->m_max < m) { if (fhv->p0_row != NULL) xfree(fhv->p0_row); if (fhv->p0_col != NULL) xfree(fhv->p0_col); if (fhv->cc_ind != NULL) xfree(fhv->cc_ind); if (fhv->cc_val != NULL) xfree(fhv->cc_val); fhv->m_max = m + 100; fhv->p0_row = xcalloc(1+fhv->m_max, sizeof(int)); fhv->p0_col = xcalloc(1+fhv->m_max, sizeof(int)); fhv->cc_ind = xcalloc(1+fhv->m_max, sizeof(int)); fhv->cc_val = xcalloc(1+fhv->m_max, sizeof(double)); } /* try to factorize the basis matrix */ switch (luf_factorize(fhv->luf, m, col, info)) { case 0: break; case LUF_ESING: ret = FHV_ESING; goto done; case LUF_ECOND: ret = FHV_ECOND; goto done; default: xassert(fhv != fhv); } /* the basis matrix has been successfully factorized */ fhv->valid = 1; /* H := I */ fhv->hh_nfs = 0; /* P0 := P */ memcpy(&fhv->p0_row[1], &fhv->luf->pp_row[1], sizeof(int) * m); memcpy(&fhv->p0_col[1], &fhv->luf->pp_col[1], sizeof(int) * m); /* currently H has no factors */ fhv->nnz_h = 0; ret = 0; done: /* return to the calling program */ return ret; } /*********************************************************************** * NAME * * fhv_h_solve - solve system H*x = b or H'*x = b * * SYNOPSIS * * #include "glpfhv.h" * void fhv_h_solve(FHV *fhv, int tr, double x[]); * * DESCRIPTION * * The routine fhv_h_solve solves either the system H*x = b (if the * flag tr is zero) or the system H'*x = b (if the flag tr is non-zero), * where the matrix H is a component of the factorization specified by * the parameter fhv, H' is a matrix transposed to H. * * On entry the array x should contain elements of the right-hand side * vector b in locations x[1], ..., x[m], where m is the order of the * matrix H. On exit this array will contain elements of the solution * vector x in the same locations. */ void fhv_h_solve(FHV *fhv, int tr, double x[]) { int nfs = fhv->hh_nfs; int *hh_ind = fhv->hh_ind; int *hh_ptr = fhv->hh_ptr; int *hh_len = fhv->hh_len; int *sv_ind = fhv->luf->sv_ind; double *sv_val = fhv->luf->sv_val; int i, k, beg, end, ptr; double temp; if (!fhv->valid) xfault("fhv_h_solve: the factorization is not valid\n"); if (!tr) { /* solve the system H*x = b */ for (k = 1; k <= nfs; k++) { i = hh_ind[k]; temp = x[i]; beg = hh_ptr[k]; end = beg + hh_len[k] - 1; for (ptr = beg; ptr <= end; ptr++) temp -= sv_val[ptr] * x[sv_ind[ptr]]; x[i] = temp; } } else { /* solve the system H'*x = b */ for (k = nfs; k >= 1; k--) { i = hh_ind[k]; temp = x[i]; if (temp == 0.0) continue; beg = hh_ptr[k]; end = beg + hh_len[k] - 1; for (ptr = beg; ptr <= end; ptr++) x[sv_ind[ptr]] -= sv_val[ptr] * temp; } } return; } /*********************************************************************** * NAME * * fhv_ftran - perform forward transformation (solve system B*x = b) * * SYNOPSIS * * #include "glpfhv.h" * void fhv_ftran(FHV *fhv, double x[]); * * DESCRIPTION * * The routine fhv_ftran performs forward transformation, i.e. solves * the system B*x = b, where B is the basis matrix, x is the vector of * unknowns to be computed, b is the vector of right-hand sides. * * On entry elements of the vector b should be stored in dense format * in locations x[1], ..., x[m], where m is the number of rows. On exit * the routine stores elements of the vector x in the same locations. */ void fhv_ftran(FHV *fhv, double x[]) { int *pp_row = fhv->luf->pp_row; int *pp_col = fhv->luf->pp_col; int *p0_row = fhv->p0_row; int *p0_col = fhv->p0_col; if (!fhv->valid) xfault("fhv_ftran: the factorization is not valid\n"); /* B = F*H*V, therefore inv(B) = inv(V)*inv(H)*inv(F) */ fhv->luf->pp_row = p0_row; fhv->luf->pp_col = p0_col; luf_f_solve(fhv->luf, 0, x); fhv->luf->pp_row = pp_row; fhv->luf->pp_col = pp_col; fhv_h_solve(fhv, 0, x); luf_v_solve(fhv->luf, 0, x); return; } /*********************************************************************** * NAME * * fhv_btran - perform backward transformation (solve system B'*x = b) * * SYNOPSIS * * #include "glpfhv.h" * void fhv_btran(FHV *fhv, double x[]); * * DESCRIPTION * * The routine fhv_btran performs backward transformation, i.e. solves * the system B'*x = b, where B' is a matrix transposed to the basis * matrix B, x is the vector of unknowns to be computed, b is the vector * of right-hand sides. * * On entry elements of the vector b should be stored in dense format * in locations x[1], ..., x[m], where m is the number of rows. On exit * the routine stores elements of the vector x in the same locations. */ void fhv_btran(FHV *fhv, double x[]) { int *pp_row = fhv->luf->pp_row; int *pp_col = fhv->luf->pp_col; int *p0_row = fhv->p0_row; int *p0_col = fhv->p0_col; if (!fhv->valid) xfault("fhv_btran: the factorization is not valid\n"); /* B = F*H*V, therefore inv(B') = inv(F')*inv(H')*inv(V') */ luf_v_solve(fhv->luf, 1, x); fhv_h_solve(fhv, 1, x); fhv->luf->pp_row = p0_row; fhv->luf->pp_col = p0_col; luf_f_solve(fhv->luf, 1, x); fhv->luf->pp_row = pp_row; fhv->luf->pp_col = pp_col; return; } /*********************************************************************** * NAME * * fhv_update_it - update LP basis factorization * * SYNOPSIS * * #include "glpfhv.h" * int fhv_update_it(FHV *fhv, int j, int len, const int ind[], * const double val[]); * * DESCRIPTION * * The routine fhv_update_it updates the factorization of the basis * matrix B after replacing its j-th column by a new vector. * * The parameter j specifies the number of column of B, which has been * replaced, 1 <= j <= m, where m is the order of B. * * Row indices and numerical values of non-zero elements of the new * column of B should be placed in locations ind[1], ..., ind[len] and * val[1], ..., val[len], resp., where len is the number of non-zeros * in the column. Neither zero nor duplicate elements are allowed. * * RETURNS * * 0 The factorization has been successfully updated. * * FHV_ESING * The adjacent basis matrix is structurally singular, since after * changing j-th column of matrix V by the new column (see algorithm * below) the case k1 > k2 occured. * * FHV_ECHECK * The factorization is inaccurate, since after transforming k2-th * row of matrix U = P*V*Q, its diagonal element u[k2,k2] is zero or * close to zero, * * FHV_ELIMIT * Maximal number of H factors has been reached. * * FHV_EROOM * Overflow of the sparse vector area. * * In case of non-zero return code the factorization becomes invalid. * It should not be used until it has been recomputed with the routine * fhv_factorize. * * ALGORITHM * * The routine fhv_update_it is based on the transformation proposed by * Forrest and Tomlin. * * Let j-th column of the basis matrix B have been replaced by new * column B[j]. In order to keep the equality B = F*H*V j-th column of * matrix V should be replaced by the column inv(F*H)*B[j]. * * From the standpoint of matrix U = P*V*Q, replacement of j-th column * of matrix V is equivalent to replacement of k1-th column of matrix U, * where k1 is determined by permutation matrix Q. Thus, matrix U loses * its upper triangular form and becomes the following: * * 1 k1 k2 m * 1 x x * x x x x x x x * . x * x x x x x x x * k1 . . * x x x x x x x * . . * x x x x x x x * . . * . x x x x x x * . . * . . x x x x x * . . * . . . x x x x * k2 . . * . . . . x x x * . . . . . . . . x x * m . . . . . . . . . x * * where row index k2 corresponds to the lowest non-zero element of * k1-th column. * * The routine moves rows and columns k1+1, k1+2, ..., k2 of matrix U * by one position to the left and upwards and moves k1-th row and k1-th * column to position k2. As the result of such symmetric permutations * matrix U becomes the following: * * 1 k1 k2 m * 1 x x x x x x x * x x * . x x x x x x * x x * k1 . . x x x x x * x x * . . . x x x x * x x * . . . . x x x * x x * . . . . . x x * x x * . . . . . . x * x x * k2 . . x x x x x * x x * . . . . . . . . x x * m . . . . . . . . . x * * Then the routine performs gaussian elimination to eliminate elements * u[k2,k1], u[k2,k1+1], ..., u[k2,k2-1] using diagonal elements * u[k1,k1], u[k1+1,k1+1], ..., u[k2-1,k2-1] as pivots in the same way * as described in comments to the routine luf_factorize (see the module * GLPLUF). Note that actually all operations are performed on matrix V, * not on matrix U. During the elimination process the routine permutes * neither rows nor columns, so only k2-th row of matrix U is changed. * * To keep the main equality B = F*H*V, each time when the routine * applies elementary gaussian transformation to the transformed row of * matrix V (which corresponds to k2-th row of matrix U), it also adds * a new element (gaussian multiplier) to the current row-like factor * of matrix H, which corresponds to the transformed row of matrix V. */ int fhv_update_it(FHV *fhv, int j, int len, const int ind[], const double val[]) { int m = fhv->m; LUF *luf = fhv->luf; int *vr_ptr = luf->vr_ptr; int *vr_len = luf->vr_len; int *vr_cap = luf->vr_cap; double *vr_piv = luf->vr_piv; int *vc_ptr = luf->vc_ptr; int *vc_len = luf->vc_len; int *vc_cap = luf->vc_cap; int *pp_row = luf->pp_row; int *pp_col = luf->pp_col; int *qq_row = luf->qq_row; int *qq_col = luf->qq_col; int *sv_ind = luf->sv_ind; double *sv_val = luf->sv_val; double *work = luf->work; double eps_tol = luf->eps_tol; int *hh_ind = fhv->hh_ind; int *hh_ptr = fhv->hh_ptr; int *hh_len = fhv->hh_len; int *p0_row = fhv->p0_row; int *p0_col = fhv->p0_col; int *cc_ind = fhv->cc_ind; double *cc_val = fhv->cc_val; double upd_tol = fhv->upd_tol; int i, i_beg, i_end, i_ptr, j_beg, j_end, j_ptr, k, k1, k2, p, q, p_beg, p_end, p_ptr, ptr, ret; double f, temp; if (!fhv->valid) xfault("fhv_update_it: the factorization is not valid\n"); if (!(1 <= j && j <= m)) xfault("fhv_update_it: j = %d; column number out of range\n", j); /* check if the new factor of matrix H can be created */ if (fhv->hh_nfs == fhv->hh_max) { /* maximal number of updates has been reached */ fhv->valid = 0; ret = FHV_ELIMIT; goto done; } /* convert new j-th column of B to dense format */ for (i = 1; i <= m; i++) cc_val[i] = 0.0; for (k = 1; k <= len; k++) { i = ind[k]; if (!(1 <= i && i <= m)) xfault("fhv_update_it: ind[%d] = %d; row number out of rang" "e\n", k, i); if (cc_val[i] != 0.0) xfault("fhv_update_it: ind[%d] = %d; duplicate row index no" "t allowed\n", k, i); if (val[k] == 0.0) xfault("fhv_update_it: val[%d] = %g; zero element not allow" "ed\n", k, val[k]); cc_val[i] = val[k]; } /* new j-th column of V := inv(F * H) * (new B[j]) */ fhv->luf->pp_row = p0_row; fhv->luf->pp_col = p0_col; luf_f_solve(fhv->luf, 0, cc_val); fhv->luf->pp_row = pp_row; fhv->luf->pp_col = pp_col; fhv_h_solve(fhv, 0, cc_val); /* convert new j-th column of V to sparse format */ len = 0; for (i = 1; i <= m; i++) { temp = cc_val[i]; if (temp == 0.0 || fabs(temp) < eps_tol) continue; len++, cc_ind[len] = i, cc_val[len] = temp; } /* clear old content of j-th column of matrix V */ j_beg = vc_ptr[j]; j_end = j_beg + vc_len[j] - 1; for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) { /* get row index of v[i,j] */ i = sv_ind[j_ptr]; /* find v[i,j] in the i-th row */ i_beg = vr_ptr[i]; i_end = i_beg + vr_len[i] - 1; for (i_ptr = i_beg; sv_ind[i_ptr] != j; i_ptr++) /* nop */; xassert(i_ptr <= i_end); /* remove v[i,j] from the i-th row */ sv_ind[i_ptr] = sv_ind[i_end]; sv_val[i_ptr] = sv_val[i_end]; vr_len[i]--; } /* now j-th column of matrix V is empty */ luf->nnz_v -= vc_len[j]; vc_len[j] = 0; /* add new elements of j-th column of matrix V to corresponding row lists; determine indices k1 and k2 */ k1 = qq_row[j], k2 = 0; for (ptr = 1; ptr <= len; ptr++) { /* get row index of v[i,j] */ i = cc_ind[ptr]; /* at least one unused location is needed in i-th row */ if (vr_len[i] + 1 > vr_cap[i]) { if (luf_enlarge_row(luf, i, vr_len[i] + 10)) { /* overflow of the sparse vector area */ fhv->valid = 0; luf->new_sva = luf->sv_size + luf->sv_size; xassert(luf->new_sva > luf->sv_size); ret = FHV_EROOM; goto done; } } /* add v[i,j] to i-th row */ i_ptr = vr_ptr[i] + vr_len[i]; sv_ind[i_ptr] = j; sv_val[i_ptr] = cc_val[ptr]; vr_len[i]++; /* adjust index k2 */ if (k2 < pp_col[i]) k2 = pp_col[i]; } /* capacity of j-th column (which is currently empty) should be not less than len locations */ if (vc_cap[j] < len) { if (luf_enlarge_col(luf, j, len)) { /* overflow of the sparse vector area */ fhv->valid = 0; luf->new_sva = luf->sv_size + luf->sv_size; xassert(luf->new_sva > luf->sv_size); ret = FHV_EROOM; goto done; } } /* add new elements of matrix V to j-th column list */ j_ptr = vc_ptr[j]; memmove(&sv_ind[j_ptr], &cc_ind[1], len * sizeof(int)); memmove(&sv_val[j_ptr], &cc_val[1], len * sizeof(double)); vc_len[j] = len; luf->nnz_v += len; /* if k1 > k2, diagonal element u[k2,k2] of matrix U is zero and therefore the adjacent basis matrix is structurally singular */ if (k1 > k2) { fhv->valid = 0; ret = FHV_ESING; goto done; } /* perform implicit symmetric permutations of rows and columns of matrix U */ i = pp_row[k1], j = qq_col[k1]; for (k = k1; k < k2; k++) { pp_row[k] = pp_row[k+1], pp_col[pp_row[k]] = k; qq_col[k] = qq_col[k+1], qq_row[qq_col[k]] = k; } pp_row[k2] = i, pp_col[i] = k2; qq_col[k2] = j, qq_row[j] = k2; /* now i-th row of the matrix V is k2-th row of matrix U; since no pivoting is used, only this row will be transformed */ /* copy elements of i-th row of matrix V to the working array and remove these elements from matrix V */ for (j = 1; j <= m; j++) work[j] = 0.0; i_beg = vr_ptr[i]; i_end = i_beg + vr_len[i] - 1; for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) { /* get column index of v[i,j] */ j = sv_ind[i_ptr]; /* store v[i,j] to the working array */ work[j] = sv_val[i_ptr]; /* find v[i,j] in the j-th column */ j_beg = vc_ptr[j]; j_end = j_beg + vc_len[j] - 1; for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++) /* nop */; xassert(j_ptr <= j_end); /* remove v[i,j] from the j-th column */ sv_ind[j_ptr] = sv_ind[j_end]; sv_val[j_ptr] = sv_val[j_end]; vc_len[j]--; } /* now i-th row of matrix V is empty */ luf->nnz_v -= vr_len[i]; vr_len[i] = 0; /* create the next row-like factor of the matrix H; this factor corresponds to i-th (transformed) row */ fhv->hh_nfs++; hh_ind[fhv->hh_nfs] = i; /* hh_ptr[] will be set later */ hh_len[fhv->hh_nfs] = 0; /* up to (k2 - k1) free locations are needed to add new elements to the non-trivial row of the row-like factor */ if (luf->sv_end - luf->sv_beg < k2 - k1) { luf_defrag_sva(luf); if (luf->sv_end - luf->sv_beg < k2 - k1) { /* overflow of the sparse vector area */ fhv->valid = luf->valid = 0; luf->new_sva = luf->sv_size + luf->sv_size; xassert(luf->new_sva > luf->sv_size); ret = FHV_EROOM; goto done; } } /* eliminate subdiagonal elements of matrix U */ for (k = k1; k < k2; k++) { /* v[p,q] = u[k,k] */ p = pp_row[k], q = qq_col[k]; /* this is the crucial point, where even tiny non-zeros should not be dropped */ if (work[q] == 0.0) continue; /* compute gaussian multiplier f = v[i,q] / v[p,q] */ f = work[q] / vr_piv[p]; /* perform gaussian transformation: (i-th row) := (i-th row) - f * (p-th row) in order to eliminate v[i,q] = u[k2,k] */ p_beg = vr_ptr[p]; p_end = p_beg + vr_len[p] - 1; for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) work[sv_ind[p_ptr]] -= f * sv_val[p_ptr]; /* store new element (gaussian multiplier that corresponds to p-th row) in the current row-like factor */ luf->sv_end--; sv_ind[luf->sv_end] = p; sv_val[luf->sv_end] = f; hh_len[fhv->hh_nfs]++; } /* set pointer to the current row-like factor of the matrix H (if no elements were added to this factor, it is unity matrix and therefore can be discarded) */ if (hh_len[fhv->hh_nfs] == 0) fhv->hh_nfs--; else { hh_ptr[fhv->hh_nfs] = luf->sv_end; fhv->nnz_h += hh_len[fhv->hh_nfs]; } /* store new pivot which corresponds to u[k2,k2] */ vr_piv[i] = work[qq_col[k2]]; /* new elements of i-th row of matrix V (which are non-diagonal elements u[k2,k2+1], ..., u[k2,m] of matrix U = P*V*Q) now are contained in the working array; add them to matrix V */ len = 0; for (k = k2+1; k <= m; k++) { /* get column index and value of v[i,j] = u[k2,k] */ j = qq_col[k]; temp = work[j]; /* if v[i,j] is close to zero, skip it */ if (fabs(temp) < eps_tol) continue; /* at least one unused location is needed in j-th column */ if (vc_len[j] + 1 > vc_cap[j]) { if (luf_enlarge_col(luf, j, vc_len[j] + 10)) { /* overflow of the sparse vector area */ fhv->valid = 0; luf->new_sva = luf->sv_size + luf->sv_size; xassert(luf->new_sva > luf->sv_size); ret = FHV_EROOM; goto done; } } /* add v[i,j] to j-th column */ j_ptr = vc_ptr[j] + vc_len[j]; sv_ind[j_ptr] = i; sv_val[j_ptr] = temp; vc_len[j]++; /* also store v[i,j] to the auxiliary array */ len++, cc_ind[len] = j, cc_val[len] = temp; } /* capacity of i-th row (which is currently empty) should be not less than len locations */ if (vr_cap[i] < len) { if (luf_enlarge_row(luf, i, len)) { /* overflow of the sparse vector area */ fhv->valid = 0; luf->new_sva = luf->sv_size + luf->sv_size; xassert(luf->new_sva > luf->sv_size); ret = FHV_EROOM; goto done; } } /* add new elements to i-th row list */ i_ptr = vr_ptr[i]; memmove(&sv_ind[i_ptr], &cc_ind[1], len * sizeof(int)); memmove(&sv_val[i_ptr], &cc_val[1], len * sizeof(double)); vr_len[i] = len; luf->nnz_v += len; /* updating is finished; check that diagonal element u[k2,k2] is not very small in absolute value among other elements in k2-th row and k2-th column of matrix U = P*V*Q */ /* temp = max(|u[k2,*]|, |u[*,k2]|) */ temp = 0.0; /* walk through k2-th row of U which is i-th row of V */ i = pp_row[k2]; i_beg = vr_ptr[i]; i_end = i_beg + vr_len[i] - 1; for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) if (temp < fabs(sv_val[i_ptr])) temp = fabs(sv_val[i_ptr]); /* walk through k2-th column of U which is j-th column of V */ j = qq_col[k2]; j_beg = vc_ptr[j]; j_end = j_beg + vc_len[j] - 1; for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) if (temp < fabs(sv_val[j_ptr])) temp = fabs(sv_val[j_ptr]); /* check that u[k2,k2] is not very small */ if (fabs(vr_piv[i]) < upd_tol * temp) { /* the factorization seems to be inaccurate and therefore must be recomputed */ fhv->valid = 0; ret = FHV_ECHECK; goto done; } /* the factorization has been successfully updated */ ret = 0; done: /* return to the calling program */ return ret; } /*********************************************************************** * NAME * * fhv_delete_it - delete LP basis factorization * * SYNOPSIS * * #include "glpfhv.h" * void fhv_delete_it(FHV *fhv); * * DESCRIPTION * * The routine fhv_delete_it deletes LP basis factorization specified * by the parameter fhv and frees all memory allocated to this program * object. */ void fhv_delete_it(FHV *fhv) { luf_delete_it(fhv->luf); if (fhv->hh_ind != NULL) xfree(fhv->hh_ind); if (fhv->hh_ptr != NULL) xfree(fhv->hh_ptr); if (fhv->hh_len != NULL) xfree(fhv->hh_len); if (fhv->p0_row != NULL) xfree(fhv->p0_row); if (fhv->p0_col != NULL) xfree(fhv->p0_col); if (fhv->cc_ind != NULL) xfree(fhv->cc_ind); if (fhv->cc_val != NULL) xfree(fhv->cc_val); xfree(fhv); return; } /* eof */