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
|
/* sgf.h (sparse Gaussian factorizer) */
/***********************************************************************
* This code is part of GLPK (GNU Linear Programming Kit).
*
* Copyright (C) 2012-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 SGF_H
#define SGF_H
#include "luf.h"
typedef struct SGF SGF;
struct SGF
{ /* sparse Gaussian factorizer workspace */
LUF *luf;
/* LU-factorization being computed */
/*--------------------------------------------------------------*/
/* to efficiently choose pivot elements according to Markowitz
* strategy, the search technique proposed by Iain Duff is used;
* it is based on using two families of sets {R[0], ..., R[n]}
* and {C[0], ..., C[n]}, where R[k] and C[k], 0 <= k <= n, are,
* respectively, sets of rows and columns of the active submatrix
* of matrix V having k non-zeros (i.e. whose length is k); each
* set R[k] and C[k] is implemented as a doubly linked list */
int *rs_head; /* int rs_head[1+n]; */
/* rs_head[k], 0 <= k <= n, is the number of first row, which
* has k non-zeros in the active submatrix */
int *rs_prev; /* int rs_prev[1+n]; */
/* rs_prev[0] is not used;
* rs_prev[i], 1 <= i <= n, is the number of previous row, which
* has the same number of non-zeros as i-th row;
* rs_prev[i] < 0 means that i-th row is inactive */
int *rs_next; /* int rs_next[1+n]; */
/* rs_next[0] is not used;
* rs_next[i], 1 <= i <= n, is the number of next row, which has
* the same number of non-zeros as i-th row;
* rs_next[i] < 0 means that i-th row is inactive */
int *cs_head; /* int cs_head[1+n]; */
/* cs_head[k], 0 <= k <= n, is the number of first column, which
* has k non-zeros in the active submatrix */
int *cs_prev; /* int cs_prev[1+n]; */
/* cs_prev[0] is not used;
* cs_prev[j], 1 <= j <= n, is the number of previous column,
* which has the same number of non-zeros as j-th column;
* cs_prev[j] < 0 means that j-th column is inactive */
int *cs_next; /* int cs_next[1+n]; */
/* cs_next[0] is not used;
* cs_next[j], 1 <= j <= n, is the number of next column, which
* has the same number of non-zeros as j-th column;
* cs_next[j] < 0 means that j-th column is inactive */
/* NOTE: cs_prev[j] = cs_next[j] = j means that j-th column was
* temporarily removed from corresponding set C[k] by the
* pivoting routine according to Uwe Suhl's heuristic */
/*--------------------------------------------------------------*/
/* working arrays */
double *vr_max; /* int vr_max[1+n]; */
/* vr_max[0] is not used;
* vr_max[i], 1 <= i <= n, is used only if i-th row of matrix V
* is active (i.e. belongs to the active submatrix), and is the
* largest magnitude of elements in that row; if vr_max[i] < 0,
* the largest magnitude is unknown yet */
char *flag; /* char flag[1+n]; */
/* boolean working array */
double *work; /* double work[1+n]; */
/* floating-point working array */
/*--------------------------------------------------------------*/
/* control parameters */
int updat;
/* if this flag is set, the matrix V is assumed to be updatable;
* in this case factorized (non-active) part of V is stored in
* the left part of SVA rather than in its right part */
double piv_tol;
/* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j]
* of the active submatrix fits to be pivot if it satisfies to
* the stability criterion |v[i,j]| >= piv_tol * max |v[i,*]|,
* i.e. if it is not very small in the magnitude among other
* elements in the same row; decreasing this parameter gives
* better sparsity at the expense of numerical accuracy and vice
* versa */
int piv_lim;
/* maximal allowable number of pivot candidates to be considered;
* if piv_lim pivot candidates have been considered, the pivoting
* routine terminates the search with the best candidate found */
int suhl;
/* if this flag is set, the pivoting routine applies a heuristic
* proposed by Uwe Suhl: if a column of the active submatrix has
* no eligible pivot candidates (i.e. all its elements do not
* satisfy to the stability criterion), the routine excludes it
* from futher consideration until it becomes column singleton;
* in many cases this allows reducing the time needed to choose
* the pivot */
double eps_tol;
/* epsilon tolerance; each element of the active submatrix, whose
* magnitude is less than eps_tol, is replaced by exact zero */
#if 0 /* FIXME */
double den_lim;
/* density limit; if the density of the active submatrix reaches
* this limit, the factorization routine switches from sparse to
* dense mode */
#endif
};
#define sgf_activate_row(i) \
do \
{ int len = vr_len[i]; \
rs_prev[i] = 0; \
rs_next[i] = rs_head[len]; \
if (rs_next[i] != 0) \
rs_prev[rs_next[i]] = i; \
rs_head[len] = i; \
} while (0)
/* include i-th row of matrix V in active set R[len] */
#define sgf_deactivate_row(i) \
do \
{ if (rs_prev[i] == 0) \
rs_head[vr_len[i]] = rs_next[i]; \
else \
rs_next[rs_prev[i]] = rs_next[i]; \
if (rs_next[i] == 0) \
; \
else \
rs_prev[rs_next[i]] = rs_prev[i]; \
rs_prev[i] = rs_next[i] = -1; \
} while (0)
/* remove i-th row of matrix V from active set R[len] */
#define sgf_activate_col(j) \
do \
{ int len = vc_len[j]; \
cs_prev[j] = 0; \
cs_next[j] = cs_head[len]; \
if (cs_next[j] != 0) \
cs_prev[cs_next[j]] = j; \
cs_head[len] = j; \
} while (0)
/* include j-th column of matrix V in active set C[len] */
#define sgf_deactivate_col(j) \
do \
{ if (cs_prev[j] == 0) \
cs_head[vc_len[j]] = cs_next[j]; \
else \
cs_next[cs_prev[j]] = cs_next[j]; \
if (cs_next[j] == 0) \
; \
else \
cs_prev[cs_next[j]] = cs_prev[j]; \
cs_prev[j] = cs_next[j] = -1; \
} while (0)
/* remove j-th column of matrix V from active set C[len] */
#define sgf_reduce_nuc _glp_sgf_reduce_nuc
int sgf_reduce_nuc(LUF *luf, int *k1, int *k2, int cnt[/*1+n*/],
int list[/*1+n*/]);
/* initial reordering to minimize nucleus size */
#define sgf_singl_phase _glp_sgf_singl_phase
int sgf_singl_phase(LUF *luf, int k1, int k2, int updat,
int ind[/*1+n*/], double val[/*1+n*/]);
/* compute LU-factorization (singleton phase) */
#define sgf_choose_pivot _glp_sgf_choose_pivot
int sgf_choose_pivot(SGF *sgf, int *p, int *q);
/* choose pivot element v[p,q] */
#define sgf_eliminate _glp_sgf_eliminate
int sgf_eliminate(SGF *sgf, int p, int q);
/* perform gaussian elimination */
#define sgf_dense_lu _glp_sgf_dense_lu
int sgf_dense_lu(int n, double a[], int r[], int c[], double eps);
/* compute dense LU-factorization with full pivoting */
#define sgf_dense_phase _glp_sgf_dense_phase
int sgf_dense_phase(LUF *luf, int k, int updat);
/* compute LU-factorization (dense phase) */
#define sgf_factorize _glp_sgf_factorize
int sgf_factorize(SGF *sgf, int singl);
/* compute LU-factorization (main routine) */
#endif
/* eof */
|
