1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glplux.h	Sun Nov 06 20:59:10 2011 +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, 2011 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 */