alpar@1: /* glplux.c */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, alpar@1: * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, alpar@1: * Moscow Aviation Institute, Moscow, Russia. All rights reserved. alpar@1: * E-mail: . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #include "glplux.h" alpar@1: #define xfault xerror alpar@1: #define dmp_create_poolx(size) dmp_create_pool() alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // lux_create - create LU-factorization. alpar@1: // alpar@1: // SYNOPSIS alpar@1: // alpar@1: // #include "glplux.h" alpar@1: // LUX *lux_create(int n); alpar@1: // alpar@1: // DESCRIPTION alpar@1: // alpar@1: // The routine lux_create creates LU-factorization data structure for alpar@1: // a matrix of the order n. Initially the factorization corresponds to alpar@1: // the unity matrix (F = V = P = Q = I, so A = I). alpar@1: // alpar@1: // RETURNS alpar@1: // alpar@1: // The routine returns a pointer to the created LU-factorization data alpar@1: // structure, which represents the unity matrix of the order n. */ alpar@1: alpar@1: LUX *lux_create(int n) alpar@1: { LUX *lux; alpar@1: int k; alpar@1: if (n < 1) alpar@1: xfault("lux_create: n = %d; invalid parameter\n", n); alpar@1: lux = xmalloc(sizeof(LUX)); alpar@1: lux->n = n; alpar@1: lux->pool = dmp_create_poolx(sizeof(LUXELM)); alpar@1: lux->F_row = xcalloc(1+n, sizeof(LUXELM *)); alpar@1: lux->F_col = xcalloc(1+n, sizeof(LUXELM *)); alpar@1: lux->V_piv = xcalloc(1+n, sizeof(mpq_t)); alpar@1: lux->V_row = xcalloc(1+n, sizeof(LUXELM *)); alpar@1: lux->V_col = xcalloc(1+n, sizeof(LUXELM *)); alpar@1: lux->P_row = xcalloc(1+n, sizeof(int)); alpar@1: lux->P_col = xcalloc(1+n, sizeof(int)); alpar@1: lux->Q_row = xcalloc(1+n, sizeof(int)); alpar@1: lux->Q_col = xcalloc(1+n, sizeof(int)); alpar@1: for (k = 1; k <= n; k++) alpar@1: { lux->F_row[k] = lux->F_col[k] = NULL; alpar@1: mpq_init(lux->V_piv[k]); alpar@1: mpq_set_si(lux->V_piv[k], 1, 1); alpar@1: lux->V_row[k] = lux->V_col[k] = NULL; alpar@1: lux->P_row[k] = lux->P_col[k] = k; alpar@1: lux->Q_row[k] = lux->Q_col[k] = k; alpar@1: } alpar@1: lux->rank = n; alpar@1: return lux; alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // initialize - initialize LU-factorization data structures. alpar@1: // alpar@1: // This routine initializes data structures for subsequent computing alpar@1: // the LU-factorization of a given matrix A, which is specified by the alpar@1: // formal routine col. On exit V = A and F = P = Q = I, where I is the alpar@1: // unity matrix. */ alpar@1: alpar@1: static void initialize(LUX *lux, int (*col)(void *info, int j, alpar@1: int ind[], mpq_t val[]), void *info, LUXWKA *wka) alpar@1: { int n = lux->n; alpar@1: DMP *pool = lux->pool; alpar@1: LUXELM **F_row = lux->F_row; alpar@1: LUXELM **F_col = lux->F_col; alpar@1: mpq_t *V_piv = lux->V_piv; alpar@1: LUXELM **V_row = lux->V_row; alpar@1: LUXELM **V_col = lux->V_col; alpar@1: int *P_row = lux->P_row; alpar@1: int *P_col = lux->P_col; alpar@1: int *Q_row = lux->Q_row; alpar@1: int *Q_col = lux->Q_col; alpar@1: int *R_len = wka->R_len; alpar@1: int *R_head = wka->R_head; alpar@1: int *R_prev = wka->R_prev; alpar@1: int *R_next = wka->R_next; alpar@1: int *C_len = wka->C_len; alpar@1: int *C_head = wka->C_head; alpar@1: int *C_prev = wka->C_prev; alpar@1: int *C_next = wka->C_next; alpar@1: LUXELM *fij, *vij; alpar@1: int i, j, k, len, *ind; alpar@1: mpq_t *val; alpar@1: /* F := I */ alpar@1: for (i = 1; i <= n; i++) alpar@1: { while (F_row[i] != NULL) alpar@1: { fij = F_row[i], F_row[i] = fij->r_next; alpar@1: mpq_clear(fij->val); alpar@1: dmp_free_atom(pool, fij, sizeof(LUXELM)); alpar@1: } alpar@1: } alpar@1: for (j = 1; j <= n; j++) F_col[j] = NULL; alpar@1: /* V := 0 */ alpar@1: for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1); alpar@1: for (i = 1; i <= n; i++) alpar@1: { while (V_row[i] != NULL) alpar@1: { vij = V_row[i], V_row[i] = vij->r_next; alpar@1: mpq_clear(vij->val); alpar@1: dmp_free_atom(pool, vij, sizeof(LUXELM)); alpar@1: } alpar@1: } alpar@1: for (j = 1; j <= n; j++) V_col[j] = NULL; alpar@1: /* V := A */ alpar@1: ind = xcalloc(1+n, sizeof(int)); alpar@1: val = xcalloc(1+n, sizeof(mpq_t)); alpar@1: for (k = 1; k <= n; k++) mpq_init(val[k]); alpar@1: for (j = 1; j <= n; j++) alpar@1: { /* obtain j-th column of matrix A */ alpar@1: len = col(info, j, ind, val); alpar@1: if (!(0 <= len && len <= n)) alpar@1: xfault("lux_decomp: j = %d: len = %d; invalid column length" alpar@1: "\n", j, len); alpar@1: /* copy elements of j-th column to matrix V */ alpar@1: for (k = 1; k <= len; k++) alpar@1: { /* get row index of a[i,j] */ alpar@1: i = ind[k]; alpar@1: if (!(1 <= i && i <= n)) alpar@1: xfault("lux_decomp: j = %d: i = %d; row index out of ran" alpar@1: "ge\n", j, i); alpar@1: /* check for duplicate indices */ alpar@1: if (V_row[i] != NULL && V_row[i]->j == j) alpar@1: xfault("lux_decomp: j = %d: i = %d; duplicate row indice" alpar@1: "s not allowed\n", j, i); alpar@1: /* check for zero value */ alpar@1: if (mpq_sgn(val[k]) == 0) alpar@1: xfault("lux_decomp: j = %d: i = %d; zero elements not al" alpar@1: "lowed\n", j, i); alpar@1: /* add new element v[i,j] = a[i,j] to V */ alpar@1: vij = dmp_get_atom(pool, sizeof(LUXELM)); alpar@1: vij->i = i, vij->j = j; alpar@1: mpq_init(vij->val); alpar@1: mpq_set(vij->val, val[k]); alpar@1: vij->r_prev = NULL; alpar@1: vij->r_next = V_row[i]; alpar@1: vij->c_prev = NULL; alpar@1: vij->c_next = V_col[j]; alpar@1: if (vij->r_next != NULL) vij->r_next->r_prev = vij; alpar@1: if (vij->c_next != NULL) vij->c_next->c_prev = vij; alpar@1: V_row[i] = V_col[j] = vij; alpar@1: } alpar@1: } alpar@1: xfree(ind); alpar@1: for (k = 1; k <= n; k++) mpq_clear(val[k]); alpar@1: xfree(val); alpar@1: /* P := Q := I */ alpar@1: for (k = 1; k <= n; k++) alpar@1: P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k; alpar@1: /* the rank of A and V is not determined yet */ alpar@1: lux->rank = -1; alpar@1: /* initially the entire matrix V is active */ alpar@1: /* determine its row lengths */ alpar@1: for (i = 1; i <= n; i++) alpar@1: { len = 0; alpar@1: for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++; alpar@1: R_len[i] = len; alpar@1: } alpar@1: /* build linked lists of active rows */ alpar@1: for (len = 0; len <= n; len++) R_head[len] = 0; alpar@1: for (i = 1; i <= n; i++) alpar@1: { len = R_len[i]; alpar@1: R_prev[i] = 0; alpar@1: R_next[i] = R_head[len]; alpar@1: if (R_next[i] != 0) R_prev[R_next[i]] = i; alpar@1: R_head[len] = i; alpar@1: } alpar@1: /* determine its column lengths */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { len = 0; alpar@1: for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++; alpar@1: C_len[j] = len; alpar@1: } alpar@1: /* build linked lists of active columns */ alpar@1: for (len = 0; len <= n; len++) C_head[len] = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { len = C_len[j]; alpar@1: C_prev[j] = 0; alpar@1: C_next[j] = C_head[len]; alpar@1: if (C_next[j] != 0) C_prev[C_next[j]] = j; alpar@1: C_head[len] = j; alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // find_pivot - choose a pivot element. alpar@1: // alpar@1: // This routine chooses a pivot element v[p,q] in the active submatrix alpar@1: // of matrix U = P*V*Q. alpar@1: // alpar@1: // It is assumed that on entry the matrix U has the following partially alpar@1: // triangularized form: alpar@1: // alpar@1: // 1 k n alpar@1: // 1 x x x x x x x x x x alpar@1: // . x x x x x x x x x alpar@1: // . . x x x x x x x x alpar@1: // . . . x x x x x x x alpar@1: // k . . . . * * * * * * alpar@1: // . . . . * * * * * * alpar@1: // . . . . * * * * * * alpar@1: // . . . . * * * * * * alpar@1: // . . . . * * * * * * alpar@1: // n . . . . * * * * * * alpar@1: // alpar@1: // where rows and columns k, k+1, ..., n belong to the active submatrix alpar@1: // (elements of the active submatrix are marked by '*'). alpar@1: // alpar@1: // Since the matrix U = P*V*Q is not stored, the routine works with the alpar@1: // matrix V. It is assumed that the row-wise representation corresponds alpar@1: // to the matrix V, but the column-wise representation corresponds to alpar@1: // the active submatrix of the matrix V, i.e. elements of the matrix V, alpar@1: // which does not belong to the active submatrix, are missing from the alpar@1: // column linked lists. It is also assumed that each active row of the alpar@1: // matrix V is in the set R[len], where len is number of non-zeros in alpar@1: // the row, and each active column of the matrix V is in the set C[len], alpar@1: // where len is number of non-zeros in the column (in the latter case alpar@1: // only elements of the active submatrix are counted; such elements are alpar@1: // marked by '*' on the figure above). alpar@1: // alpar@1: // Due to exact arithmetic any non-zero element of the active submatrix alpar@1: // can be chosen as a pivot. However, to keep sparsity of the matrix V alpar@1: // the routine uses Markowitz strategy, trying to choose such element alpar@1: // v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1), alpar@1: // where nr[p] and nc[q] are the number of non-zero elements, resp., in alpar@1: // p-th row and in q-th column of the active submatrix. alpar@1: // alpar@1: // In order to reduce the search, i.e. not to walk through all elements alpar@1: // of the active submatrix, the routine exploits a technique proposed by alpar@1: // I.Duff. This technique is based on using the sets R[len] and C[len] alpar@1: // of active rows and columns. alpar@1: // alpar@1: // On exit the routine returns a pointer to a pivot v[p,q] chosen, or alpar@1: // NULL, if the active submatrix is empty. */ alpar@1: alpar@1: static LUXELM *find_pivot(LUX *lux, LUXWKA *wka) alpar@1: { int n = lux->n; alpar@1: LUXELM **V_row = lux->V_row; alpar@1: LUXELM **V_col = lux->V_col; alpar@1: int *R_len = wka->R_len; alpar@1: int *R_head = wka->R_head; alpar@1: int *R_next = wka->R_next; alpar@1: int *C_len = wka->C_len; alpar@1: int *C_head = wka->C_head; alpar@1: int *C_next = wka->C_next; alpar@1: LUXELM *piv, *some, *vij; alpar@1: int i, j, len, min_len, ncand, piv_lim = 5; alpar@1: double best, cost; alpar@1: /* nothing is chosen so far */ alpar@1: piv = NULL, best = DBL_MAX, ncand = 0; alpar@1: /* if in the active submatrix there is a column that has the only alpar@1: non-zero (column singleton), choose it as a pivot */ alpar@1: j = C_head[1]; alpar@1: if (j != 0) alpar@1: { xassert(C_len[j] == 1); alpar@1: piv = V_col[j]; alpar@1: xassert(piv != NULL && piv->c_next == NULL); alpar@1: goto done; alpar@1: } alpar@1: /* if in the active submatrix there is a row that has the only alpar@1: non-zero (row singleton), choose it as a pivot */ alpar@1: i = R_head[1]; alpar@1: if (i != 0) alpar@1: { xassert(R_len[i] == 1); alpar@1: piv = V_row[i]; alpar@1: xassert(piv != NULL && piv->r_next == NULL); alpar@1: goto done; alpar@1: } alpar@1: /* there are no singletons in the active submatrix; walk through alpar@1: other non-empty rows and columns */ alpar@1: for (len = 2; len <= n; len++) alpar@1: { /* consider active columns having len non-zeros */ alpar@1: for (j = C_head[len]; j != 0; j = C_next[j]) alpar@1: { /* j-th column has len non-zeros */ alpar@1: /* find an element in the row of minimal length */ alpar@1: some = NULL, min_len = INT_MAX; alpar@1: for (vij = V_col[j]; vij != NULL; vij = vij->c_next) alpar@1: { if (min_len > R_len[vij->i]) alpar@1: some = vij, min_len = R_len[vij->i]; alpar@1: /* if Markowitz cost of this element is not greater than alpar@1: (len-1)**2, it can be chosen right now; this heuristic alpar@1: reduces the search and works well in many cases */ alpar@1: if (min_len <= len) alpar@1: { piv = some; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* j-th column has been scanned */ alpar@1: /* the minimal element found is a next pivot candidate */ alpar@1: xassert(some != NULL); alpar@1: ncand++; alpar@1: /* compute its Markowitz cost */ alpar@1: cost = (double)(min_len - 1) * (double)(len - 1); alpar@1: /* choose between the current candidate and this element */ alpar@1: if (cost < best) piv = some, best = cost; alpar@1: /* if piv_lim candidates have been considered, there is a alpar@1: doubt that a much better candidate exists; therefore it alpar@1: is the time to terminate the search */ alpar@1: if (ncand == piv_lim) goto done; alpar@1: } alpar@1: /* now consider active rows having len non-zeros */ alpar@1: for (i = R_head[len]; i != 0; i = R_next[i]) alpar@1: { /* i-th row has len non-zeros */ alpar@1: /* find an element in the column of minimal length */ alpar@1: some = NULL, min_len = INT_MAX; alpar@1: for (vij = V_row[i]; vij != NULL; vij = vij->r_next) alpar@1: { if (min_len > C_len[vij->j]) alpar@1: some = vij, min_len = C_len[vij->j]; alpar@1: /* if Markowitz cost of this element is not greater than alpar@1: (len-1)**2, it can be chosen right now; this heuristic alpar@1: reduces the search and works well in many cases */ alpar@1: if (min_len <= len) alpar@1: { piv = some; alpar@1: goto done; alpar@1: } alpar@1: } alpar@1: /* i-th row has been scanned */ alpar@1: /* the minimal element found is a next pivot candidate */ alpar@1: xassert(some != NULL); alpar@1: ncand++; alpar@1: /* compute its Markowitz cost */ alpar@1: cost = (double)(len - 1) * (double)(min_len - 1); alpar@1: /* choose between the current candidate and this element */ alpar@1: if (cost < best) piv = some, best = cost; alpar@1: /* if piv_lim candidates have been considered, there is a alpar@1: doubt that a much better candidate exists; therefore it alpar@1: is the time to terminate the search */ alpar@1: if (ncand == piv_lim) goto done; alpar@1: } alpar@1: } alpar@1: done: /* bring the pivot v[p,q] to the factorizing routine */ alpar@1: return piv; alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // eliminate - perform gaussian elimination. alpar@1: // alpar@1: // This routine performs elementary gaussian transformations in order alpar@1: // to eliminate subdiagonal elements in the k-th column of the matrix alpar@1: // U = P*V*Q using the pivot element u[k,k], where k is the number of alpar@1: // the current elimination step. alpar@1: // alpar@1: // The parameter piv specifies the pivot element v[p,q] = u[k,k]. alpar@1: // alpar@1: // Each time when the routine applies the elementary transformation to alpar@1: // a non-pivot row of the matrix V, it stores the corresponding element alpar@1: // to the matrix F in order to keep the main equality A = F*V. alpar@1: // alpar@1: // The routine assumes that on entry the matrices L = P*F*inv(P) and alpar@1: // U = P*V*Q are the following: alpar@1: // alpar@1: // 1 k 1 k n alpar@1: // 1 1 . . . . . . . . . 1 x x x x x x x x x x alpar@1: // x 1 . . . . . . . . . x x x x x x x x x alpar@1: // x x 1 . . . . . . . . . x x x x x x x x alpar@1: // x x x 1 . . . . . . . . . x x x x x x x alpar@1: // k x x x x 1 . . . . . k . . . . * * * * * * alpar@1: // x x x x _ 1 . . . . . . . . # * * * * * alpar@1: // x x x x _ . 1 . . . . . . . # * * * * * alpar@1: // x x x x _ . . 1 . . . . . . # * * * * * alpar@1: // x x x x _ . . . 1 . . . . . # * * * * * alpar@1: // n x x x x _ . . . . 1 n . . . . # * * * * * alpar@1: // alpar@1: // matrix L matrix U alpar@1: // alpar@1: // where rows and columns of the matrix U with numbers k, k+1, ..., n alpar@1: // form the active submatrix (eliminated elements are marked by '#' and alpar@1: // other elements of the active submatrix are marked by '*'). Note that alpar@1: // each eliminated non-zero element u[i,k] of the matrix U gives the alpar@1: // corresponding element l[i,k] of the matrix L (marked by '_'). alpar@1: // alpar@1: // Actually all operations are performed on the matrix V. Should note alpar@1: // that the row-wise representation corresponds to the matrix V, but the alpar@1: // column-wise representation corresponds to the active submatrix of the alpar@1: // matrix V, i.e. elements of the matrix V, which doesn't belong to the alpar@1: // active submatrix, are missing from the column linked lists. alpar@1: // alpar@1: // Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal alpar@1: // elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies alpar@1: // the following elementary gaussian transformations: alpar@1: // alpar@1: // (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V), alpar@1: // alpar@1: // where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier. alpar@1: // alpar@1: // Additionally, in order to keep the main equality A = F*V, each time alpar@1: // when the routine applies the transformation to i-th row of the matrix alpar@1: // V, it also adds f[i,p] as a new element to the matrix F. alpar@1: // alpar@1: // IMPORTANT: On entry the working arrays flag and work should contain alpar@1: // zeros. This status is provided by the routine on exit. */ alpar@1: alpar@1: static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[], alpar@1: mpq_t work[]) alpar@1: { DMP *pool = lux->pool; alpar@1: LUXELM **F_row = lux->F_row; alpar@1: LUXELM **F_col = lux->F_col; alpar@1: mpq_t *V_piv = lux->V_piv; alpar@1: LUXELM **V_row = lux->V_row; alpar@1: LUXELM **V_col = lux->V_col; alpar@1: int *R_len = wka->R_len; alpar@1: int *R_head = wka->R_head; alpar@1: int *R_prev = wka->R_prev; alpar@1: int *R_next = wka->R_next; alpar@1: int *C_len = wka->C_len; alpar@1: int *C_head = wka->C_head; alpar@1: int *C_prev = wka->C_prev; alpar@1: int *C_next = wka->C_next; alpar@1: LUXELM *fip, *vij, *vpj, *viq, *next; alpar@1: mpq_t temp; alpar@1: int i, j, p, q; alpar@1: mpq_init(temp); alpar@1: /* determine row and column indices of the pivot v[p,q] */ alpar@1: xassert(piv != NULL); alpar@1: p = piv->i, q = piv->j; alpar@1: /* remove p-th (pivot) row from the active set; it will never alpar@1: return there */ alpar@1: if (R_prev[p] == 0) alpar@1: R_head[R_len[p]] = R_next[p]; alpar@1: else alpar@1: R_next[R_prev[p]] = R_next[p]; alpar@1: if (R_next[p] == 0) alpar@1: ; alpar@1: else alpar@1: R_prev[R_next[p]] = R_prev[p]; alpar@1: /* remove q-th (pivot) column from the active set; it will never alpar@1: return there */ alpar@1: if (C_prev[q] == 0) alpar@1: C_head[C_len[q]] = C_next[q]; alpar@1: else alpar@1: C_next[C_prev[q]] = C_next[q]; alpar@1: if (C_next[q] == 0) alpar@1: ; alpar@1: else alpar@1: C_prev[C_next[q]] = C_prev[q]; alpar@1: /* store the pivot value in a separate array */ alpar@1: mpq_set(V_piv[p], piv->val); alpar@1: /* remove the pivot from p-th row */ alpar@1: if (piv->r_prev == NULL) alpar@1: V_row[p] = piv->r_next; alpar@1: else alpar@1: piv->r_prev->r_next = piv->r_next; alpar@1: if (piv->r_next == NULL) alpar@1: ; alpar@1: else alpar@1: piv->r_next->r_prev = piv->r_prev; alpar@1: R_len[p]--; alpar@1: /* remove the pivot from q-th column */ alpar@1: if (piv->c_prev == NULL) alpar@1: V_col[q] = piv->c_next; alpar@1: else alpar@1: piv->c_prev->c_next = piv->c_next; alpar@1: if (piv->c_next == NULL) alpar@1: ; alpar@1: else alpar@1: piv->c_next->c_prev = piv->c_prev; alpar@1: C_len[q]--; alpar@1: /* free the space occupied by the pivot */ alpar@1: mpq_clear(piv->val); alpar@1: dmp_free_atom(pool, piv, sizeof(LUXELM)); alpar@1: /* walk through p-th (pivot) row, which already does not contain alpar@1: the pivot v[p,q], and do the following... */ alpar@1: for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) alpar@1: { /* get column index of v[p,j] */ alpar@1: j = vpj->j; alpar@1: /* store v[p,j] in the working array */ alpar@1: flag[j] = 1; alpar@1: mpq_set(work[j], vpj->val); alpar@1: /* remove j-th column from the active set; it will return there alpar@1: later with a new length */ alpar@1: if (C_prev[j] == 0) alpar@1: C_head[C_len[j]] = C_next[j]; alpar@1: else alpar@1: C_next[C_prev[j]] = C_next[j]; alpar@1: if (C_next[j] == 0) alpar@1: ; alpar@1: else alpar@1: C_prev[C_next[j]] = C_prev[j]; alpar@1: /* v[p,j] leaves the active submatrix, so remove it from j-th alpar@1: column; however, v[p,j] is kept in p-th row */ alpar@1: if (vpj->c_prev == NULL) alpar@1: V_col[j] = vpj->c_next; alpar@1: else alpar@1: vpj->c_prev->c_next = vpj->c_next; alpar@1: if (vpj->c_next == NULL) alpar@1: ; alpar@1: else alpar@1: vpj->c_next->c_prev = vpj->c_prev; alpar@1: C_len[j]--; alpar@1: } alpar@1: /* now walk through q-th (pivot) column, which already does not alpar@1: contain the pivot v[p,q], and perform gaussian elimination */ alpar@1: while (V_col[q] != NULL) alpar@1: { /* element v[i,q] has to be eliminated */ alpar@1: viq = V_col[q]; alpar@1: /* get row index of v[i,q] */ alpar@1: i = viq->i; alpar@1: /* remove i-th row from the active set; later it will return alpar@1: there with a new length */ alpar@1: if (R_prev[i] == 0) alpar@1: R_head[R_len[i]] = R_next[i]; alpar@1: else alpar@1: R_next[R_prev[i]] = R_next[i]; alpar@1: if (R_next[i] == 0) alpar@1: ; alpar@1: else alpar@1: R_prev[R_next[i]] = R_prev[i]; alpar@1: /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and alpar@1: store it in the matrix F */ alpar@1: fip = dmp_get_atom(pool, sizeof(LUXELM)); alpar@1: fip->i = i, fip->j = p; alpar@1: mpq_init(fip->val); alpar@1: mpq_div(fip->val, viq->val, V_piv[p]); alpar@1: fip->r_prev = NULL; alpar@1: fip->r_next = F_row[i]; alpar@1: fip->c_prev = NULL; alpar@1: fip->c_next = F_col[p]; alpar@1: if (fip->r_next != NULL) fip->r_next->r_prev = fip; alpar@1: if (fip->c_next != NULL) fip->c_next->c_prev = fip; alpar@1: F_row[i] = F_col[p] = fip; alpar@1: /* v[i,q] has to be eliminated, so remove it from i-th row */ alpar@1: if (viq->r_prev == NULL) alpar@1: V_row[i] = viq->r_next; alpar@1: else alpar@1: viq->r_prev->r_next = viq->r_next; alpar@1: if (viq->r_next == NULL) alpar@1: ; alpar@1: else alpar@1: viq->r_next->r_prev = viq->r_prev; alpar@1: R_len[i]--; alpar@1: /* and also from q-th column */ alpar@1: V_col[q] = viq->c_next; alpar@1: C_len[q]--; alpar@1: /* free the space occupied by v[i,q] */ alpar@1: mpq_clear(viq->val); alpar@1: dmp_free_atom(pool, viq, sizeof(LUXELM)); alpar@1: /* perform gaussian transformation: alpar@1: (i-th row) := (i-th row) - f[i,p] * (p-th row) alpar@1: note that now p-th row, which is in the working array, alpar@1: does not contain the pivot v[p,q], and i-th row does not alpar@1: contain the element v[i,q] to be eliminated */ alpar@1: /* walk through i-th row and transform existing non-zero alpar@1: elements */ alpar@1: for (vij = V_row[i]; vij != NULL; vij = next) alpar@1: { next = vij->r_next; alpar@1: /* get column index of v[i,j] */ alpar@1: j = vij->j; alpar@1: /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */ alpar@1: if (flag[j]) alpar@1: { /* v[p,j] != 0 */ alpar@1: flag[j] = 0; alpar@1: mpq_mul(temp, fip->val, work[j]); alpar@1: mpq_sub(vij->val, vij->val, temp); alpar@1: if (mpq_sgn(vij->val) == 0) alpar@1: { /* new v[i,j] is zero, so remove it from the active alpar@1: submatrix */ alpar@1: /* remove v[i,j] from i-th row */ alpar@1: if (vij->r_prev == NULL) alpar@1: V_row[i] = vij->r_next; alpar@1: else alpar@1: vij->r_prev->r_next = vij->r_next; alpar@1: if (vij->r_next == NULL) alpar@1: ; alpar@1: else alpar@1: vij->r_next->r_prev = vij->r_prev; alpar@1: R_len[i]--; alpar@1: /* remove v[i,j] from j-th column */ alpar@1: if (vij->c_prev == NULL) alpar@1: V_col[j] = vij->c_next; alpar@1: else alpar@1: vij->c_prev->c_next = vij->c_next; alpar@1: if (vij->c_next == NULL) alpar@1: ; alpar@1: else alpar@1: vij->c_next->c_prev = vij->c_prev; alpar@1: C_len[j]--; alpar@1: /* free the space occupied by v[i,j] */ alpar@1: mpq_clear(vij->val); alpar@1: dmp_free_atom(pool, vij, sizeof(LUXELM)); alpar@1: } alpar@1: } alpar@1: } alpar@1: /* now flag is the pattern of the set v[p,*] \ v[i,*] */ alpar@1: /* walk through p-th (pivot) row and create new elements in alpar@1: i-th row, which appear due to fill-in */ alpar@1: for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) alpar@1: { j = vpj->j; alpar@1: if (flag[j]) alpar@1: { /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and alpar@1: add it to i-th row and j-th column */ alpar@1: vij = dmp_get_atom(pool, sizeof(LUXELM)); alpar@1: vij->i = i, vij->j = j; alpar@1: mpq_init(vij->val); alpar@1: mpq_mul(vij->val, fip->val, work[j]); alpar@1: mpq_neg(vij->val, vij->val); alpar@1: vij->r_prev = NULL; alpar@1: vij->r_next = V_row[i]; alpar@1: vij->c_prev = NULL; alpar@1: vij->c_next = V_col[j]; alpar@1: if (vij->r_next != NULL) vij->r_next->r_prev = vij; alpar@1: if (vij->c_next != NULL) vij->c_next->c_prev = vij; alpar@1: V_row[i] = V_col[j] = vij; alpar@1: R_len[i]++, C_len[j]++; alpar@1: } alpar@1: else alpar@1: { /* there is no fill-in, because v[i,j] already exists in alpar@1: i-th row; restore the flag, which was reset before */ alpar@1: flag[j] = 1; alpar@1: } alpar@1: } alpar@1: /* now i-th row has been completely transformed and can return alpar@1: to the active set with a new length */ alpar@1: R_prev[i] = 0; alpar@1: R_next[i] = R_head[R_len[i]]; alpar@1: if (R_next[i] != 0) R_prev[R_next[i]] = i; alpar@1: R_head[R_len[i]] = i; alpar@1: } alpar@1: /* at this point q-th (pivot) column must be empty */ alpar@1: xassert(C_len[q] == 0); alpar@1: /* walk through p-th (pivot) row again and do the following... */ alpar@1: for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) alpar@1: { /* get column index of v[p,j] */ alpar@1: j = vpj->j; alpar@1: /* erase v[p,j] from the working array */ alpar@1: flag[j] = 0; alpar@1: mpq_set_si(work[j], 0, 1); alpar@1: /* now j-th column has been completely transformed, so it can alpar@1: return to the active list with a new length */ alpar@1: C_prev[j] = 0; alpar@1: C_next[j] = C_head[C_len[j]]; alpar@1: if (C_next[j] != 0) C_prev[C_next[j]] = j; alpar@1: C_head[C_len[j]] = j; alpar@1: } alpar@1: mpq_clear(temp); alpar@1: /* return to the factorizing routine */ alpar@1: return; alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // lux_decomp - compute LU-factorization. alpar@1: // alpar@1: // SYNOPSIS alpar@1: // alpar@1: // #include "glplux.h" alpar@1: // int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], alpar@1: // mpq_t val[]), void *info); alpar@1: // alpar@1: // DESCRIPTION alpar@1: // alpar@1: // The routine lux_decomp computes LU-factorization of a given square alpar@1: // matrix A. alpar@1: // alpar@1: // The parameter lux specifies LU-factorization data structure built by alpar@1: // means of the routine lux_create. alpar@1: // alpar@1: // The formal routine col specifies the original matrix A. In order to alpar@1: // obtain j-th column of the matrix A the routine lux_decomp calls the alpar@1: // routine col with the parameter j (1 <= j <= n, where n is the order alpar@1: // of A). In response the routine col should store row indices and alpar@1: // numerical values of non-zero elements of j-th column of A to the alpar@1: // locations ind[1], ..., ind[len] and val[1], ..., val[len], resp., alpar@1: // where len is the number of non-zeros in j-th column, which should be alpar@1: // returned on exit. Neiter zero nor duplicate elements are allowed. alpar@1: // alpar@1: // The parameter info is a transit pointer passed to the formal routine alpar@1: // col; it can be used for various purposes. alpar@1: // alpar@1: // RETURNS alpar@1: // alpar@1: // The routine lux_decomp returns the singularity flag. Zero flag means alpar@1: // that the original matrix A is non-singular while non-zero flag means alpar@1: // that A is (exactly!) singular. alpar@1: // alpar@1: // Note that LU-factorization is valid in both cases, however, in case alpar@1: // of singularity some rows of the matrix V (including pivot elements) alpar@1: // will be empty. alpar@1: // alpar@1: // REPAIRING SINGULAR MATRIX alpar@1: // alpar@1: // If the routine lux_decomp returns non-zero flag, it provides all alpar@1: // necessary information that can be used for "repairing" the matrix A, alpar@1: // where "repairing" means replacing linearly dependent columns of the alpar@1: // matrix A by appropriate columns of the unity matrix. This feature is alpar@1: // needed when the routine lux_decomp is used for reinverting the basis alpar@1: // matrix within the simplex method procedure. alpar@1: // alpar@1: // On exit linearly dependent columns of the matrix U have the numbers alpar@1: // rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A alpar@1: // stored by the routine to the member lux->rank. The correspondence alpar@1: // between columns of A and U is the same as between columns of V and U. alpar@1: // Thus, linearly dependent columns of the matrix A have the numbers alpar@1: // Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array alpar@1: // representing the permutation matrix Q in column-like format. It is alpar@1: // understood that each j-th linearly dependent column of the matrix U alpar@1: // should be replaced by the unity vector, where all elements are zero alpar@1: // except the unity diagonal element u[j,j]. On the other hand j-th row alpar@1: // of the matrix U corresponds to the row of the matrix V (and therefore alpar@1: // of the matrix A) with the number P_row[j], where P_row is an array alpar@1: // representing the permutation matrix P in row-like format. Thus, each alpar@1: // j-th linearly dependent column of the matrix U should be replaced by alpar@1: // a column of the unity matrix with the number P_row[j]. alpar@1: // alpar@1: // The code that repairs the matrix A may look like follows: alpar@1: // alpar@1: // for (j = rank+1; j <= n; j++) alpar@1: // { replace column Q_col[j] of the matrix A by column P_row[j] of alpar@1: // the unity matrix; alpar@1: // } alpar@1: // alpar@1: // where rank, P_row, and Q_col are members of the structure LUX. */ alpar@1: alpar@1: int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], alpar@1: mpq_t val[]), void *info) alpar@1: { int n = lux->n; alpar@1: LUXELM **V_row = lux->V_row; alpar@1: LUXELM **V_col = lux->V_col; alpar@1: int *P_row = lux->P_row; alpar@1: int *P_col = lux->P_col; alpar@1: int *Q_row = lux->Q_row; alpar@1: int *Q_col = lux->Q_col; alpar@1: LUXELM *piv, *vij; alpar@1: LUXWKA *wka; alpar@1: int i, j, k, p, q, t, *flag; alpar@1: mpq_t *work; alpar@1: /* allocate working area */ alpar@1: wka = xmalloc(sizeof(LUXWKA)); alpar@1: wka->R_len = xcalloc(1+n, sizeof(int)); alpar@1: wka->R_head = xcalloc(1+n, sizeof(int)); alpar@1: wka->R_prev = xcalloc(1+n, sizeof(int)); alpar@1: wka->R_next = xcalloc(1+n, sizeof(int)); alpar@1: wka->C_len = xcalloc(1+n, sizeof(int)); alpar@1: wka->C_head = xcalloc(1+n, sizeof(int)); alpar@1: wka->C_prev = xcalloc(1+n, sizeof(int)); alpar@1: wka->C_next = xcalloc(1+n, sizeof(int)); alpar@1: /* initialize LU-factorization data structures */ alpar@1: initialize(lux, col, info, wka); alpar@1: /* allocate working arrays */ alpar@1: flag = xcalloc(1+n, sizeof(int)); alpar@1: work = xcalloc(1+n, sizeof(mpq_t)); alpar@1: for (k = 1; k <= n; k++) alpar@1: { flag[k] = 0; alpar@1: mpq_init(work[k]); alpar@1: } alpar@1: /* main elimination loop */ alpar@1: for (k = 1; k <= n; k++) alpar@1: { /* choose a pivot element v[p,q] */ alpar@1: piv = find_pivot(lux, wka); alpar@1: if (piv == NULL) alpar@1: { /* no pivot can be chosen, because the active submatrix is alpar@1: empty */ alpar@1: break; alpar@1: } alpar@1: /* determine row and column indices of the pivot element */ alpar@1: p = piv->i, q = piv->j; alpar@1: /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th alpar@1: rows and k-th and j'-th columns of the matrix U = P*V*Q to alpar@1: move the element u[i',j'] to the position u[k,k] */ alpar@1: i = P_col[p], j = Q_row[q]; alpar@1: xassert(k <= i && i <= n && k <= j && j <= n); alpar@1: /* permute k-th and i-th rows of the matrix U */ alpar@1: t = P_row[k]; alpar@1: P_row[i] = t, P_col[t] = i; alpar@1: P_row[k] = p, P_col[p] = k; alpar@1: /* permute k-th and j-th columns of the matrix U */ alpar@1: t = Q_col[k]; alpar@1: Q_col[j] = t, Q_row[t] = j; alpar@1: Q_col[k] = q, Q_row[q] = k; alpar@1: /* eliminate subdiagonal elements of k-th column of the matrix alpar@1: U = P*V*Q using the pivot element u[k,k] = v[p,q] */ alpar@1: eliminate(lux, wka, piv, flag, work); alpar@1: } alpar@1: /* determine the rank of A (and V) */ alpar@1: lux->rank = k - 1; alpar@1: /* free working arrays */ alpar@1: xfree(flag); alpar@1: for (k = 1; k <= n; k++) mpq_clear(work[k]); alpar@1: xfree(work); alpar@1: /* build column lists of the matrix V using its row lists */ alpar@1: for (j = 1; j <= n; j++) alpar@1: xassert(V_col[j] == NULL); alpar@1: for (i = 1; i <= n; i++) alpar@1: { for (vij = V_row[i]; vij != NULL; vij = vij->r_next) alpar@1: { j = vij->j; alpar@1: vij->c_prev = NULL; alpar@1: vij->c_next = V_col[j]; alpar@1: if (vij->c_next != NULL) vij->c_next->c_prev = vij; alpar@1: V_col[j] = vij; alpar@1: } alpar@1: } alpar@1: /* free working area */ alpar@1: xfree(wka->R_len); alpar@1: xfree(wka->R_head); alpar@1: xfree(wka->R_prev); alpar@1: xfree(wka->R_next); alpar@1: xfree(wka->C_len); alpar@1: xfree(wka->C_head); alpar@1: xfree(wka->C_prev); alpar@1: xfree(wka->C_next); alpar@1: xfree(wka); alpar@1: /* return to the calling program */ alpar@1: return (lux->rank < n); alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // lux_f_solve - solve system F*x = b or F'*x = b. alpar@1: // alpar@1: // SYNOPSIS alpar@1: // alpar@1: // #include "glplux.h" alpar@1: // void lux_f_solve(LUX *lux, int tr, mpq_t x[]); alpar@1: // alpar@1: // DESCRIPTION alpar@1: // alpar@1: // The routine lux_f_solve solves either the system F*x = b (if the alpar@1: // flag tr is zero) or the system F'*x = b (if the flag tr is non-zero), alpar@1: // where the matrix F is a component of LU-factorization specified by alpar@1: // the parameter lux, F' is a matrix transposed to F. alpar@1: // alpar@1: // On entry the array x should contain elements of the right-hand side alpar@1: // vector b in locations x[1], ..., x[n], where n is the order of the alpar@1: // matrix F. On exit this array will contain elements of the solution alpar@1: // vector x in the same locations. */ alpar@1: alpar@1: void lux_f_solve(LUX *lux, int tr, mpq_t x[]) alpar@1: { int n = lux->n; alpar@1: LUXELM **F_row = lux->F_row; alpar@1: LUXELM **F_col = lux->F_col; alpar@1: int *P_row = lux->P_row; alpar@1: LUXELM *fik, *fkj; alpar@1: int i, j, k; alpar@1: mpq_t temp; alpar@1: mpq_init(temp); alpar@1: if (!tr) alpar@1: { /* solve the system F*x = b */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { k = P_row[j]; alpar@1: if (mpq_sgn(x[k]) != 0) alpar@1: { for (fik = F_col[k]; fik != NULL; fik = fik->c_next) alpar@1: { mpq_mul(temp, fik->val, x[k]); alpar@1: mpq_sub(x[fik->i], x[fik->i], temp); alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* solve the system F'*x = b */ alpar@1: for (i = n; i >= 1; i--) alpar@1: { k = P_row[i]; alpar@1: if (mpq_sgn(x[k]) != 0) alpar@1: { for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next) alpar@1: { mpq_mul(temp, fkj->val, x[k]); alpar@1: mpq_sub(x[fkj->j], x[fkj->j], temp); alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: mpq_clear(temp); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // lux_v_solve - solve system V*x = b or V'*x = b. alpar@1: // alpar@1: // SYNOPSIS alpar@1: // alpar@1: // #include "glplux.h" alpar@1: // void lux_v_solve(LUX *lux, int tr, double x[]); alpar@1: // alpar@1: // DESCRIPTION alpar@1: // alpar@1: // The routine lux_v_solve solves either the system V*x = b (if the alpar@1: // flag tr is zero) or the system V'*x = b (if the flag tr is non-zero), alpar@1: // where the matrix V is a component of LU-factorization specified by alpar@1: // the parameter lux, V' is a matrix transposed to V. alpar@1: // alpar@1: // On entry the array x should contain elements of the right-hand side alpar@1: // vector b in locations x[1], ..., x[n], where n is the order of the alpar@1: // matrix V. On exit this array will contain elements of the solution alpar@1: // vector x in the same locations. */ alpar@1: alpar@1: void lux_v_solve(LUX *lux, int tr, mpq_t x[]) alpar@1: { int n = lux->n; alpar@1: mpq_t *V_piv = lux->V_piv; alpar@1: LUXELM **V_row = lux->V_row; alpar@1: LUXELM **V_col = lux->V_col; alpar@1: int *P_row = lux->P_row; alpar@1: int *Q_col = lux->Q_col; alpar@1: LUXELM *vij; alpar@1: int i, j, k; alpar@1: mpq_t *b, temp; alpar@1: b = xcalloc(1+n, sizeof(mpq_t)); alpar@1: for (k = 1; k <= n; k++) alpar@1: mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1); alpar@1: mpq_init(temp); alpar@1: if (!tr) alpar@1: { /* solve the system V*x = b */ alpar@1: for (k = n; k >= 1; k--) alpar@1: { i = P_row[k], j = Q_col[k]; alpar@1: if (mpq_sgn(b[i]) != 0) alpar@1: { mpq_set(x[j], b[i]); alpar@1: mpq_div(x[j], x[j], V_piv[i]); alpar@1: for (vij = V_col[j]; vij != NULL; vij = vij->c_next) alpar@1: { mpq_mul(temp, vij->val, x[j]); alpar@1: mpq_sub(b[vij->i], b[vij->i], temp); alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: else alpar@1: { /* solve the system V'*x = b */ alpar@1: for (k = 1; k <= n; k++) alpar@1: { i = P_row[k], j = Q_col[k]; alpar@1: if (mpq_sgn(b[j]) != 0) alpar@1: { mpq_set(x[i], b[j]); alpar@1: mpq_div(x[i], x[i], V_piv[i]); alpar@1: for (vij = V_row[i]; vij != NULL; vij = vij->r_next) alpar@1: { mpq_mul(temp, vij->val, x[i]); alpar@1: mpq_sub(b[vij->j], b[vij->j], temp); alpar@1: } alpar@1: } alpar@1: } alpar@1: } alpar@1: for (k = 1; k <= n; k++) mpq_clear(b[k]); alpar@1: mpq_clear(temp); alpar@1: xfree(b); alpar@1: return; alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // lux_solve - solve system A*x = b or A'*x = b. alpar@1: // alpar@1: // SYNOPSIS alpar@1: // alpar@1: // #include "glplux.h" alpar@1: // void lux_solve(LUX *lux, int tr, mpq_t x[]); alpar@1: // alpar@1: // DESCRIPTION alpar@1: // alpar@1: // The routine lux_solve solves either the system A*x = b (if the flag alpar@1: // tr is zero) or the system A'*x = b (if the flag tr is non-zero), alpar@1: // where the parameter lux specifies LU-factorization of the matrix A, alpar@1: // A' is a matrix transposed to A. alpar@1: // alpar@1: // On entry the array x should contain elements of the right-hand side alpar@1: // vector b in locations x[1], ..., x[n], where n is the order of the alpar@1: // matrix A. On exit this array will contain elements of the solution alpar@1: // vector x in the same locations. */ alpar@1: alpar@1: void lux_solve(LUX *lux, int tr, mpq_t x[]) alpar@1: { if (lux->rank < lux->n) alpar@1: xfault("lux_solve: LU-factorization has incomplete rank\n"); alpar@1: if (!tr) alpar@1: { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */ alpar@1: lux_f_solve(lux, 0, x); alpar@1: lux_v_solve(lux, 0, x); alpar@1: } alpar@1: else alpar@1: { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */ alpar@1: lux_v_solve(lux, 1, x); alpar@1: lux_f_solve(lux, 1, x); alpar@1: } alpar@1: return; alpar@1: } alpar@1: alpar@1: /*---------------------------------------------------------------------- alpar@1: // lux_delete - delete LU-factorization. alpar@1: // alpar@1: // SYNOPSIS alpar@1: // alpar@1: // #include "glplux.h" alpar@1: // void lux_delete(LUX *lux); alpar@1: // alpar@1: // DESCRIPTION alpar@1: // alpar@1: // The routine lux_delete deletes LU-factorization data structure, alpar@1: // which the parameter lux points to, freeing all the memory allocated alpar@1: // to this object. */ alpar@1: alpar@1: void lux_delete(LUX *lux) alpar@1: { int n = lux->n; alpar@1: LUXELM *fij, *vij; alpar@1: int i; alpar@1: for (i = 1; i <= n; i++) alpar@1: { for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next) alpar@1: mpq_clear(fij->val); alpar@1: mpq_clear(lux->V_piv[i]); alpar@1: for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next) alpar@1: mpq_clear(vij->val); alpar@1: } alpar@1: dmp_delete_pool(lux->pool); alpar@1: xfree(lux->F_row); alpar@1: xfree(lux->F_col); alpar@1: xfree(lux->V_piv); alpar@1: xfree(lux->V_row); alpar@1: xfree(lux->V_col); alpar@1: xfree(lux->P_row); alpar@1: xfree(lux->P_col); alpar@1: xfree(lux->Q_row); alpar@1: xfree(lux->Q_col); alpar@1: xfree(lux); alpar@1: return; alpar@1: } alpar@1: alpar@1: /* eof */