1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glplux.h Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,221 @@
1.4 +/* glplux.h (LU-factorization, bignum arithmetic) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#ifndef GLPLUX_H
1.29 +#define GLPLUX_H
1.30 +
1.31 +#include "glpdmp.h"
1.32 +#include "glpgmp.h"
1.33 +
1.34 +/*----------------------------------------------------------------------
1.35 +// The structure LUX defines LU-factorization of a square matrix A,
1.36 +// which is the following quartet:
1.37 +//
1.38 +// [A] = (F, V, P, Q), (1)
1.39 +//
1.40 +// where F and V are such matrices that
1.41 +//
1.42 +// A = F * V, (2)
1.43 +//
1.44 +// and P and Q are such permutation matrices that the matrix
1.45 +//
1.46 +// L = P * F * inv(P) (3)
1.47 +//
1.48 +// is lower triangular with unity diagonal, and the matrix
1.49 +//
1.50 +// U = P * V * Q (4)
1.51 +//
1.52 +// is upper triangular. All the matrices have the order n.
1.53 +//
1.54 +// The matrices F and V are stored in row/column-wise sparse format as
1.55 +// row and column linked lists of non-zero elements. Unity elements on
1.56 +// the main diagonal of the matrix F are not stored. Pivot elements of
1.57 +// the matrix V (that correspond to diagonal elements of the matrix U)
1.58 +// are also missing from the row and column lists and stored separately
1.59 +// in an ordinary array.
1.60 +//
1.61 +// The permutation matrices P and Q are stored as ordinary arrays using
1.62 +// both row- and column-like formats.
1.63 +//
1.64 +// The matrices L and U being completely defined by the matrices F, V,
1.65 +// P, and Q are not stored explicitly.
1.66 +//
1.67 +// It is easy to show that the factorization (1)-(3) is some version of
1.68 +// LU-factorization. Indeed, from (3) and (4) it follows that:
1.69 +//
1.70 +// F = inv(P) * L * P,
1.71 +//
1.72 +// V = inv(P) * U * inv(Q),
1.73 +//
1.74 +// and substitution into (2) gives:
1.75 +//
1.76 +// A = F * V = inv(P) * L * U * inv(Q).
1.77 +//
1.78 +// For more details see the program documentation. */
1.79 +
1.80 +typedef struct LUX LUX;
1.81 +typedef struct LUXELM LUXELM;
1.82 +typedef struct LUXWKA LUXWKA;
1.83 +
1.84 +struct LUX
1.85 +{ /* LU-factorization of a square matrix */
1.86 + int n;
1.87 + /* the order of matrices A, F, V, P, Q */
1.88 + DMP *pool;
1.89 + /* memory pool for elements of matrices F and V */
1.90 + LUXELM **F_row; /* LUXELM *F_row[1+n]; */
1.91 + /* F_row[0] is not used;
1.92 + F_row[i], 1 <= i <= n, is a pointer to the list of elements in
1.93 + i-th row of matrix F (diagonal elements are not stored) */
1.94 + LUXELM **F_col; /* LUXELM *F_col[1+n]; */
1.95 + /* F_col[0] is not used;
1.96 + F_col[j], 1 <= j <= n, is a pointer to the list of elements in
1.97 + j-th column of matrix F (diagonal elements are not stored) */
1.98 + mpq_t *V_piv; /* mpq_t V_piv[1+n]; */
1.99 + /* V_piv[0] is not used;
1.100 + V_piv[p], 1 <= p <= n, is a pivot element v[p,q] corresponding
1.101 + to a diagonal element u[k,k] of matrix U = P*V*Q (used on k-th
1.102 + elimination step, k = 1, 2, ..., n) */
1.103 + LUXELM **V_row; /* LUXELM *V_row[1+n]; */
1.104 + /* V_row[0] is not used;
1.105 + V_row[i], 1 <= i <= n, is a pointer to the list of elements in
1.106 + i-th row of matrix V (except pivot elements) */
1.107 + LUXELM **V_col; /* LUXELM *V_col[1+n]; */
1.108 + /* V_col[0] is not used;
1.109 + V_col[j], 1 <= j <= n, is a pointer to the list of elements in
1.110 + j-th column of matrix V (except pivot elements) */
1.111 + int *P_row; /* int P_row[1+n]; */
1.112 + /* P_row[0] is not used;
1.113 + P_row[i] = j means that p[i,j] = 1, where p[i,j] is an element
1.114 + of permutation matrix P */
1.115 + int *P_col; /* int P_col[1+n]; */
1.116 + /* P_col[0] is not used;
1.117 + P_col[j] = i means that p[i,j] = 1, where p[i,j] is an element
1.118 + of permutation matrix P */
1.119 + /* if i-th row or column of matrix F is i'-th row or column of
1.120 + matrix L = P*F*inv(P), or if i-th row of matrix V is i'-th row
1.121 + of matrix U = P*V*Q, then P_row[i'] = i and P_col[i] = i' */
1.122 + int *Q_row; /* int Q_row[1+n]; */
1.123 + /* Q_row[0] is not used;
1.124 + Q_row[i] = j means that q[i,j] = 1, where q[i,j] is an element
1.125 + of permutation matrix Q */
1.126 + int *Q_col; /* int Q_col[1+n]; */
1.127 + /* Q_col[0] is not used;
1.128 + Q_col[j] = i means that q[i,j] = 1, where q[i,j] is an element
1.129 + of permutation matrix Q */
1.130 + /* if j-th column of matrix V is j'-th column of matrix U = P*V*Q,
1.131 + then Q_row[j] = j' and Q_col[j'] = j */
1.132 + int rank;
1.133 + /* the (exact) rank of matrices A and V */
1.134 +};
1.135 +
1.136 +struct LUXELM
1.137 +{ /* element of matrix F or V */
1.138 + int i;
1.139 + /* row index, 1 <= i <= m */
1.140 + int j;
1.141 + /* column index, 1 <= j <= n */
1.142 + mpq_t val;
1.143 + /* numeric (non-zero) element value */
1.144 + LUXELM *r_prev;
1.145 + /* pointer to previous element in the same row */
1.146 + LUXELM *r_next;
1.147 + /* pointer to next element in the same row */
1.148 + LUXELM *c_prev;
1.149 + /* pointer to previous element in the same column */
1.150 + LUXELM *c_next;
1.151 + /* pointer to next element in the same column */
1.152 +};
1.153 +
1.154 +struct LUXWKA
1.155 +{ /* working area (used only during factorization) */
1.156 + /* in order to efficiently implement Markowitz strategy and Duff
1.157 + search technique there are two families {R[0], R[1], ..., R[n]}
1.158 + and {C[0], C[1], ..., C[n]}; member R[k] is a set of active
1.159 + rows of matrix V having k non-zeros, and member C[k] is a set
1.160 + of active columns of matrix V having k non-zeros (in the active
1.161 + submatrix); each set R[k] and C[k] is implemented as a separate
1.162 + doubly linked list */
1.163 + int *R_len; /* int R_len[1+n]; */
1.164 + /* R_len[0] is not used;
1.165 + R_len[i], 1 <= i <= n, is the number of non-zero elements in
1.166 + i-th row of matrix V (that is the length of i-th row) */
1.167 + int *R_head; /* int R_head[1+n]; */
1.168 + /* R_head[k], 0 <= k <= n, is the number of a first row, which is
1.169 + active and whose length is k */
1.170 + int *R_prev; /* int R_prev[1+n]; */
1.171 + /* R_prev[0] is not used;
1.172 + R_prev[i], 1 <= i <= n, is the number of a previous row, which
1.173 + is active and has the same length as i-th row */
1.174 + int *R_next; /* int R_next[1+n]; */
1.175 + /* R_prev[0] is not used;
1.176 + R_prev[i], 1 <= i <= n, is the number of a next row, which is
1.177 + active and has the same length as i-th row */
1.178 + int *C_len; /* int C_len[1+n]; */
1.179 + /* C_len[0] is not used;
1.180 + C_len[j], 1 <= j <= n, is the number of non-zero elements in
1.181 + j-th column of the active submatrix of matrix V (that is the
1.182 + length of j-th column in the active submatrix) */
1.183 + int *C_head; /* int C_head[1+n]; */
1.184 + /* C_head[k], 0 <= k <= n, is the number of a first column, which
1.185 + is active and whose length is k */
1.186 + int *C_prev; /* int C_prev[1+n]; */
1.187 + /* C_prev[0] is not used;
1.188 + C_prev[j], 1 <= j <= n, is the number of a previous column,
1.189 + which is active and has the same length as j-th column */
1.190 + int *C_next; /* int C_next[1+n]; */
1.191 + /* C_next[0] is not used;
1.192 + C_next[j], 1 <= j <= n, is the number of a next column, which
1.193 + is active and has the same length as j-th column */
1.194 +};
1.195 +
1.196 +#define lux_create _glp_lux_create
1.197 +#define lux_decomp _glp_lux_decomp
1.198 +#define lux_f_solve _glp_lux_f_solve
1.199 +#define lux_v_solve _glp_lux_v_solve
1.200 +#define lux_solve _glp_lux_solve
1.201 +#define lux_delete _glp_lux_delete
1.202 +
1.203 +LUX *lux_create(int n);
1.204 +/* create LU-factorization */
1.205 +
1.206 +int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
1.207 + mpq_t val[]), void *info);
1.208 +/* compute LU-factorization */
1.209 +
1.210 +void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
1.211 +/* solve system F*x = b or F'*x = b */
1.212 +
1.213 +void lux_v_solve(LUX *lux, int tr, mpq_t x[]);
1.214 +/* solve system V*x = b or V'*x = b */
1.215 +
1.216 +void lux_solve(LUX *lux, int tr, mpq_t x[]);
1.217 +/* solve system A*x = b or A'*x = b */
1.218 +
1.219 +void lux_delete(LUX *lux);
1.220 +/* delete LU-factorization */
1.221 +
1.222 +#endif
1.223 +
1.224 +/* eof */