diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/draft/lux.h')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/draft/lux.h | 220 |
1 files changed, 220 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/draft/lux.h b/test/monniaux/glpk-4.65/src/draft/lux.h new file mode 100644 index 00000000..8767bb8e --- /dev/null +++ b/test/monniaux/glpk-4.65/src/draft/lux.h @@ -0,0 +1,220 @@ +/* lux.h (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/>. +***********************************************************************/ + +#ifndef LUX_H +#define LUX_H + +#include "dmp.h" +#include "mygmp.h" + +/*********************************************************************** +* The structure LUX defines LU-factorization of a square matrix A, +* which is the following quartet: +* +* [A] = (F, V, P, Q), (1) +* +* where F and V are such matrices that +* +* A = F * V, (2) +* +* and P and Q are such permutation matrices that the matrix +* +* L = P * F * inv(P) (3) +* +* is lower triangular with unity diagonal, and the matrix +* +* U = P * V * Q (4) +* +* is upper triangular. All the matrices have the order n. +* +* The matrices F and V are stored in row/column-wise sparse format as +* row and column linked lists of non-zero elements. Unity elements on +* the main diagonal of the matrix F are not stored. Pivot elements of +* the matrix V (that correspond to diagonal elements of the matrix U) +* are also missing from the row and column lists and stored separately +* in an ordinary array. +* +* The permutation matrices P and Q are stored as ordinary arrays using +* both row- and column-like formats. +* +* The matrices L and U being completely defined by the matrices F, V, +* P, and Q are not stored explicitly. +* +* It is easy to show that the factorization (1)-(3) is some version of +* LU-factorization. Indeed, from (3) and (4) it follows that: +* +* F = inv(P) * L * P, +* +* V = inv(P) * U * inv(Q), +* +* and substitution into (2) gives: +* +* A = F * V = inv(P) * L * U * inv(Q). +* +* For more details see the program documentation. */ + +typedef struct LUX LUX; +typedef struct LUXELM LUXELM; +typedef struct LUXWKA LUXWKA; + +struct LUX +{ /* LU-factorization of a square matrix */ + int n; + /* the order of matrices A, F, V, P, Q */ + DMP *pool; + /* memory pool for elements of matrices F and V */ + LUXELM **F_row; /* LUXELM *F_row[1+n]; */ + /* F_row[0] is not used; + F_row[i], 1 <= i <= n, is a pointer to the list of elements in + i-th row of matrix F (diagonal elements are not stored) */ + LUXELM **F_col; /* LUXELM *F_col[1+n]; */ + /* F_col[0] is not used; + F_col[j], 1 <= j <= n, is a pointer to the list of elements in + j-th column of matrix F (diagonal elements are not stored) */ + mpq_t *V_piv; /* mpq_t V_piv[1+n]; */ + /* V_piv[0] is not used; + V_piv[p], 1 <= p <= n, is a pivot element v[p,q] corresponding + to a diagonal element u[k,k] of matrix U = P*V*Q (used on k-th + elimination step, k = 1, 2, ..., n) */ + LUXELM **V_row; /* LUXELM *V_row[1+n]; */ + /* V_row[0] is not used; + V_row[i], 1 <= i <= n, is a pointer to the list of elements in + i-th row of matrix V (except pivot elements) */ + LUXELM **V_col; /* LUXELM *V_col[1+n]; */ + /* V_col[0] is not used; + V_col[j], 1 <= j <= n, is a pointer to the list of elements in + j-th column of matrix V (except pivot elements) */ + int *P_row; /* int P_row[1+n]; */ + /* P_row[0] is not used; + P_row[i] = j means that p[i,j] = 1, where p[i,j] is an element + of permutation matrix P */ + int *P_col; /* int P_col[1+n]; */ + /* P_col[0] is not used; + P_col[j] = i means that p[i,j] = 1, where p[i,j] is an element + of permutation matrix P */ + /* if i-th row or column of matrix F is i'-th row or column of + matrix L = P*F*inv(P), or if i-th row of matrix V is i'-th row + of matrix U = P*V*Q, then P_row[i'] = i and P_col[i] = i' */ + int *Q_row; /* int Q_row[1+n]; */ + /* Q_row[0] is not used; + Q_row[i] = j means that q[i,j] = 1, where q[i,j] is an element + of permutation matrix Q */ + int *Q_col; /* int Q_col[1+n]; */ + /* Q_col[0] is not used; + Q_col[j] = i means that q[i,j] = 1, where q[i,j] is an element + of permutation matrix Q */ + /* if j-th column of matrix V is j'-th column of matrix U = P*V*Q, + then Q_row[j] = j' and Q_col[j'] = j */ + int rank; + /* the (exact) rank of matrices A and V */ +}; + +struct LUXELM +{ /* element of matrix F or V */ + int i; + /* row index, 1 <= i <= m */ + int j; + /* column index, 1 <= j <= n */ + mpq_t val; + /* numeric (non-zero) element value */ + LUXELM *r_prev; + /* pointer to previous element in the same row */ + LUXELM *r_next; + /* pointer to next element in the same row */ + LUXELM *c_prev; + /* pointer to previous element in the same column */ + LUXELM *c_next; + /* pointer to next element in the same column */ +}; + +struct LUXWKA +{ /* working area (used only during factorization) */ + /* in order to efficiently implement Markowitz strategy and Duff + search technique there are two families {R[0], R[1], ..., R[n]} + and {C[0], C[1], ..., C[n]}; member R[k] is a set of active + rows of matrix V having k non-zeros, and member C[k] is a set + of active columns of matrix V having k non-zeros (in the active + submatrix); each set R[k] and C[k] is implemented as a separate + doubly linked list */ + int *R_len; /* int R_len[1+n]; */ + /* R_len[0] is not used; + R_len[i], 1 <= i <= n, is the number of non-zero elements in + i-th row of matrix V (that is the length of i-th row) */ + int *R_head; /* int R_head[1+n]; */ + /* R_head[k], 0 <= k <= n, is the number of a first row, which is + active and whose length is k */ + int *R_prev; /* int R_prev[1+n]; */ + /* R_prev[0] is not used; + R_prev[i], 1 <= i <= n, is the number of a previous row, which + is active and has the same length as i-th row */ + int *R_next; /* int R_next[1+n]; */ + /* R_prev[0] is not used; + R_prev[i], 1 <= i <= n, is the number of a next row, which is + active and has the same length as i-th row */ + int *C_len; /* int C_len[1+n]; */ + /* C_len[0] is not used; + C_len[j], 1 <= j <= n, is the number of non-zero elements in + j-th column of the active submatrix of matrix V (that is the + length of j-th column in the active submatrix) */ + int *C_head; /* int C_head[1+n]; */ + /* C_head[k], 0 <= k <= n, is the number of a first column, which + is active and whose length is k */ + int *C_prev; /* int C_prev[1+n]; */ + /* C_prev[0] is not used; + C_prev[j], 1 <= j <= n, is the number of a previous column, + which is active and has the same length as j-th column */ + int *C_next; /* int C_next[1+n]; */ + /* C_next[0] is not used; + C_next[j], 1 <= j <= n, is the number of a next column, which + is active and has the same length as j-th column */ +}; + +#define lux_create _glp_lux_create +LUX *lux_create(int n); +/* create LU-factorization */ + +#define lux_decomp _glp_lux_decomp +int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], + mpq_t val[]), void *info); +/* compute LU-factorization */ + +#define lux_f_solve _glp_lux_f_solve +void lux_f_solve(LUX *lux, int tr, mpq_t x[]); +/* solve system F*x = b or F'*x = b */ + +#define lux_v_solve _glp_lux_v_solve +void lux_v_solve(LUX *lux, int tr, mpq_t x[]); +/* solve system V*x = b or V'*x = b */ + +#define lux_solve _glp_lux_solve +void lux_solve(LUX *lux, int tr, mpq_t x[]); +/* solve system A*x = b or A'*x = b */ + +#define lux_delete _glp_lux_delete +void lux_delete(LUX *lux); +/* delete LU-factorization */ + +#endif + +/* eof */ |