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