1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
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 */
|
