lemon-project-template-glpk

annotate deps/glpk/src/glplux.h @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glplux.h (LU-factorization, bignum arithmetic) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #ifndef GLPLUX_H
alpar@9 26 #define GLPLUX_H
alpar@9 27
alpar@9 28 #include "glpdmp.h"
alpar@9 29 #include "glpgmp.h"
alpar@9 30
alpar@9 31 /*----------------------------------------------------------------------
alpar@9 32 // The structure LUX defines LU-factorization of a square matrix A,
alpar@9 33 // which is the following quartet:
alpar@9 34 //
alpar@9 35 // [A] = (F, V, P, Q), (1)
alpar@9 36 //
alpar@9 37 // where F and V are such matrices that
alpar@9 38 //
alpar@9 39 // A = F * V, (2)
alpar@9 40 //
alpar@9 41 // and P and Q are such permutation matrices that the matrix
alpar@9 42 //
alpar@9 43 // L = P * F * inv(P) (3)
alpar@9 44 //
alpar@9 45 // is lower triangular with unity diagonal, and the matrix
alpar@9 46 //
alpar@9 47 // U = P * V * Q (4)
alpar@9 48 //
alpar@9 49 // is upper triangular. All the matrices have the order n.
alpar@9 50 //
alpar@9 51 // The matrices F and V are stored in row/column-wise sparse format as
alpar@9 52 // row and column linked lists of non-zero elements. Unity elements on
alpar@9 53 // the main diagonal of the matrix F are not stored. Pivot elements of
alpar@9 54 // the matrix V (that correspond to diagonal elements of the matrix U)
alpar@9 55 // are also missing from the row and column lists and stored separately
alpar@9 56 // in an ordinary array.
alpar@9 57 //
alpar@9 58 // The permutation matrices P and Q are stored as ordinary arrays using
alpar@9 59 // both row- and column-like formats.
alpar@9 60 //
alpar@9 61 // The matrices L and U being completely defined by the matrices F, V,
alpar@9 62 // P, and Q are not stored explicitly.
alpar@9 63 //
alpar@9 64 // It is easy to show that the factorization (1)-(3) is some version of
alpar@9 65 // LU-factorization. Indeed, from (3) and (4) it follows that:
alpar@9 66 //
alpar@9 67 // F = inv(P) * L * P,
alpar@9 68 //
alpar@9 69 // V = inv(P) * U * inv(Q),
alpar@9 70 //
alpar@9 71 // and substitution into (2) gives:
alpar@9 72 //
alpar@9 73 // A = F * V = inv(P) * L * U * inv(Q).
alpar@9 74 //
alpar@9 75 // For more details see the program documentation. */
alpar@9 76
alpar@9 77 typedef struct LUX LUX;
alpar@9 78 typedef struct LUXELM LUXELM;
alpar@9 79 typedef struct LUXWKA LUXWKA;
alpar@9 80
alpar@9 81 struct LUX
alpar@9 82 { /* LU-factorization of a square matrix */
alpar@9 83 int n;
alpar@9 84 /* the order of matrices A, F, V, P, Q */
alpar@9 85 DMP *pool;
alpar@9 86 /* memory pool for elements of matrices F and V */
alpar@9 87 LUXELM **F_row; /* LUXELM *F_row[1+n]; */
alpar@9 88 /* F_row[0] is not used;
alpar@9 89 F_row[i], 1 <= i <= n, is a pointer to the list of elements in
alpar@9 90 i-th row of matrix F (diagonal elements are not stored) */
alpar@9 91 LUXELM **F_col; /* LUXELM *F_col[1+n]; */
alpar@9 92 /* F_col[0] is not used;
alpar@9 93 F_col[j], 1 <= j <= n, is a pointer to the list of elements in
alpar@9 94 j-th column of matrix F (diagonal elements are not stored) */
alpar@9 95 mpq_t *V_piv; /* mpq_t V_piv[1+n]; */
alpar@9 96 /* V_piv[0] is not used;
alpar@9 97 V_piv[p], 1 <= p <= n, is a pivot element v[p,q] corresponding
alpar@9 98 to a diagonal element u[k,k] of matrix U = P*V*Q (used on k-th
alpar@9 99 elimination step, k = 1, 2, ..., n) */
alpar@9 100 LUXELM **V_row; /* LUXELM *V_row[1+n]; */
alpar@9 101 /* V_row[0] is not used;
alpar@9 102 V_row[i], 1 <= i <= n, is a pointer to the list of elements in
alpar@9 103 i-th row of matrix V (except pivot elements) */
alpar@9 104 LUXELM **V_col; /* LUXELM *V_col[1+n]; */
alpar@9 105 /* V_col[0] is not used;
alpar@9 106 V_col[j], 1 <= j <= n, is a pointer to the list of elements in
alpar@9 107 j-th column of matrix V (except pivot elements) */
alpar@9 108 int *P_row; /* int P_row[1+n]; */
alpar@9 109 /* P_row[0] is not used;
alpar@9 110 P_row[i] = j means that p[i,j] = 1, where p[i,j] is an element
alpar@9 111 of permutation matrix P */
alpar@9 112 int *P_col; /* int P_col[1+n]; */
alpar@9 113 /* P_col[0] is not used;
alpar@9 114 P_col[j] = i means that p[i,j] = 1, where p[i,j] is an element
alpar@9 115 of permutation matrix P */
alpar@9 116 /* if i-th row or column of matrix F is i'-th row or column of
alpar@9 117 matrix L = P*F*inv(P), or if i-th row of matrix V is i'-th row
alpar@9 118 of matrix U = P*V*Q, then P_row[i'] = i and P_col[i] = i' */
alpar@9 119 int *Q_row; /* int Q_row[1+n]; */
alpar@9 120 /* Q_row[0] is not used;
alpar@9 121 Q_row[i] = j means that q[i,j] = 1, where q[i,j] is an element
alpar@9 122 of permutation matrix Q */
alpar@9 123 int *Q_col; /* int Q_col[1+n]; */
alpar@9 124 /* Q_col[0] is not used;
alpar@9 125 Q_col[j] = i means that q[i,j] = 1, where q[i,j] is an element
alpar@9 126 of permutation matrix Q */
alpar@9 127 /* if j-th column of matrix V is j'-th column of matrix U = P*V*Q,
alpar@9 128 then Q_row[j] = j' and Q_col[j'] = j */
alpar@9 129 int rank;
alpar@9 130 /* the (exact) rank of matrices A and V */
alpar@9 131 };
alpar@9 132
alpar@9 133 struct LUXELM
alpar@9 134 { /* element of matrix F or V */
alpar@9 135 int i;
alpar@9 136 /* row index, 1 <= i <= m */
alpar@9 137 int j;
alpar@9 138 /* column index, 1 <= j <= n */
alpar@9 139 mpq_t val;
alpar@9 140 /* numeric (non-zero) element value */
alpar@9 141 LUXELM *r_prev;
alpar@9 142 /* pointer to previous element in the same row */
alpar@9 143 LUXELM *r_next;
alpar@9 144 /* pointer to next element in the same row */
alpar@9 145 LUXELM *c_prev;
alpar@9 146 /* pointer to previous element in the same column */
alpar@9 147 LUXELM *c_next;
alpar@9 148 /* pointer to next element in the same column */
alpar@9 149 };
alpar@9 150
alpar@9 151 struct LUXWKA
alpar@9 152 { /* working area (used only during factorization) */
alpar@9 153 /* in order to efficiently implement Markowitz strategy and Duff
alpar@9 154 search technique there are two families {R[0], R[1], ..., R[n]}
alpar@9 155 and {C[0], C[1], ..., C[n]}; member R[k] is a set of active
alpar@9 156 rows of matrix V having k non-zeros, and member C[k] is a set
alpar@9 157 of active columns of matrix V having k non-zeros (in the active
alpar@9 158 submatrix); each set R[k] and C[k] is implemented as a separate
alpar@9 159 doubly linked list */
alpar@9 160 int *R_len; /* int R_len[1+n]; */
alpar@9 161 /* R_len[0] is not used;
alpar@9 162 R_len[i], 1 <= i <= n, is the number of non-zero elements in
alpar@9 163 i-th row of matrix V (that is the length of i-th row) */
alpar@9 164 int *R_head; /* int R_head[1+n]; */
alpar@9 165 /* R_head[k], 0 <= k <= n, is the number of a first row, which is
alpar@9 166 active and whose length is k */
alpar@9 167 int *R_prev; /* int R_prev[1+n]; */
alpar@9 168 /* R_prev[0] is not used;
alpar@9 169 R_prev[i], 1 <= i <= n, is the number of a previous row, which
alpar@9 170 is active and has the same length as i-th row */
alpar@9 171 int *R_next; /* int R_next[1+n]; */
alpar@9 172 /* R_prev[0] is not used;
alpar@9 173 R_prev[i], 1 <= i <= n, is the number of a next row, which is
alpar@9 174 active and has the same length as i-th row */
alpar@9 175 int *C_len; /* int C_len[1+n]; */
alpar@9 176 /* C_len[0] is not used;
alpar@9 177 C_len[j], 1 <= j <= n, is the number of non-zero elements in
alpar@9 178 j-th column of the active submatrix of matrix V (that is the
alpar@9 179 length of j-th column in the active submatrix) */
alpar@9 180 int *C_head; /* int C_head[1+n]; */
alpar@9 181 /* C_head[k], 0 <= k <= n, is the number of a first column, which
alpar@9 182 is active and whose length is k */
alpar@9 183 int *C_prev; /* int C_prev[1+n]; */
alpar@9 184 /* C_prev[0] is not used;
alpar@9 185 C_prev[j], 1 <= j <= n, is the number of a previous column,
alpar@9 186 which is active and has the same length as j-th column */
alpar@9 187 int *C_next; /* int C_next[1+n]; */
alpar@9 188 /* C_next[0] is not used;
alpar@9 189 C_next[j], 1 <= j <= n, is the number of a next column, which
alpar@9 190 is active and has the same length as j-th column */
alpar@9 191 };
alpar@9 192
alpar@9 193 #define lux_create _glp_lux_create
alpar@9 194 #define lux_decomp _glp_lux_decomp
alpar@9 195 #define lux_f_solve _glp_lux_f_solve
alpar@9 196 #define lux_v_solve _glp_lux_v_solve
alpar@9 197 #define lux_solve _glp_lux_solve
alpar@9 198 #define lux_delete _glp_lux_delete
alpar@9 199
alpar@9 200 LUX *lux_create(int n);
alpar@9 201 /* create LU-factorization */
alpar@9 202
alpar@9 203 int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
alpar@9 204 mpq_t val[]), void *info);
alpar@9 205 /* compute LU-factorization */
alpar@9 206
alpar@9 207 void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
alpar@9 208 /* solve system F*x = b or F'*x = b */
alpar@9 209
alpar@9 210 void lux_v_solve(LUX *lux, int tr, mpq_t x[]);
alpar@9 211 /* solve system V*x = b or V'*x = b */
alpar@9 212
alpar@9 213 void lux_solve(LUX *lux, int tr, mpq_t x[]);
alpar@9 214 /* solve system A*x = b or A'*x = b */
alpar@9 215
alpar@9 216 void lux_delete(LUX *lux);
alpar@9 217 /* delete LU-factorization */
alpar@9 218
alpar@9 219 #endif
alpar@9 220
alpar@9 221 /* eof */