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