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
|
/* scf.h (sparse updatable Schur-complement-based factorization) */
/***********************************************************************
* This code is part of GLPK (GNU Linear Programming Kit).
*
* Copyright (C) 2013-2014 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 SCF_H
#define SCF_H
#include "btf.h"
#include "ifu.h"
#include "luf.h"
/***********************************************************************
* The structure SCF describes sparse updatable factorization based on
* Schur complement.
*
* The SCF-factorization has the following format:
*
* ( A A1~ ) ( A0 A1 ) ( R0 ) ( S0 S )
* ( ) = P ( ) Q = P ( ) ( ) Q, (1)
* ( A2~ A3~ ) ( A2 A3 ) ( R I ) ( C )
*
* where:
*
* A is current (unsymmetric) square matrix (not stored);
*
* A1~, A2~, A3~ are some additional matrices (not stored);
*
* A0 is initial (unsymmetric) square matrix (not stored);
*
* A1, A2, A3 are some additional matrices (not stored);
*
* R0 and S0 are matrices that define factorization of the initial
* matrix A0 = R0 * S0 (stored in an invertable form);
*
* R is a matrix defined from R * S0 = A2, so R = A2 * inv(S0) (stored
* in row-wise sparse format);
*
* S is a matrix defined from R0 * S = A1, so S = inv(R0) * A1 (stored
* in column-wise sparse format);
*
* C is Schur complement (to matrix A0) defined from R * S + C = A3,
* so C = A3 - R * S = A3 - A2 * inv(A0) * A1 (stored in an invertable
* form).
*
* P, Q are permutation matrices (stored in both row- and column-like
* formats). */
typedef struct SCF SCF;
struct SCF
{ /* Schur-complement-based factorization */
int n;
/* order of current matrix A */
/*--------------------------------------------------------------*/
/* initial matrix A0 = R0 * S0 of order n0 in invertable form */
int n0;
/* order of matrix A0 */
int type;
/* type of factorization used:
* 1 - LU-factorization (R0 = F0, S0 = V0)
* 2 - BT-factorization (R0 = I, S0 = A0) */
union
{ LUF *luf; /* type = 1 */
BTF *btf; /* type = 2 */
} a0;
/* factorization of matrix A0 */
/*--------------------------------------------------------------*/
/* augmented matrix (A0, A1; A2, A3) of order n0+nn */
int nn_max;
/* maximal number of additional rows and columns in the augmented
* matrix (this limits the number of updates) */
int nn;
/* current number of additional rows and columns in the augmented
* matrix, 0 <= nn <= nn_max */
SVA *sva;
/* associated sparse vector area (SVA) used to store rows of
* matrix R and columns of matrix S */
/*--------------------------------------------------------------*/
/* nn*n0-matrix R in row-wise format */
int rr_ref;
/* reference number of sparse vector in SVA, which is the first
* row of matrix R */
#if 0 + 0
int *rr_ptr = &sva->ptr[rr_ref-1];
/* rr_ptr[0] is not used;
* rr_ptr[i], 1 <= i <= nn, is pointer to i-th row in SVA;
* rr_ptr[nn+1,...,nn_max] are reserved locations */
int *rr_len = &sva->len[rr_ref-1];
/* rr_len[0] is not used;
* rr_len[i], 1 <= i <= nn, is length of i-th row;
* rr_len[nn+1,...,nn_max] are reserved locations */
#endif
/*--------------------------------------------------------------*/
/* n0*nn-matrix S in column-wise format */
int ss_ref;
/* reference number of sparse vector in SVA, which is the first
* column of matrix S */
#if 0 + 0
int *ss_ptr = &sva->ptr[ss_ref-1];
/* ss_ptr[0] is not used;
* ss_ptr[j], 1 <= j <= nn, is pointer to j-th column in SVA;
* ss_ptr[nn+1,...,nn_max] are reserved locations */
int *ss_len = &sva->len[ss_ref-1];
/* ss_len[0] is not used;
* ss_len[j], 1 <= j <= nn, is length of j-th column;
* ss_len[nn+1,...,nn_max] are reserved locations */
#endif
/*--------------------------------------------------------------*/
/* Schur complement C of order nn in invertable form */
IFU ifu;
/* IFU-factorization of matrix C */
/*--------------------------------------------------------------*/
/* permutation matrix P of order n0+nn */
int *pp_ind; /* int pp_ind[1+n0+nn_max]; */
/* pp_ind[i] = j means that P[i,j] = 1 */
int *pp_inv; /* int pp_inv[1+n0+nn_max]; */
/* pp_inv[j] = i means that P[i,j] = 1 */
/*--------------------------------------------------------------*/
/* permutation matrix Q of order n0+nn */
int *qq_ind; /* int qq_ind[1+n0+nn_max]; */
/* qq_ind[i] = j means that Q[i,j] = 1 */
int *qq_inv; /* int qq_inv[1+n0+nn_max]; */
/* qq_inv[j] = i means that Q[i,j] = 1 */
};
#define scf_swap_q_cols(j1, j2) \
do \
{ int i1, i2; \
i1 = qq_inv[j1], i2 = qq_inv[j2]; \
qq_ind[i1] = j2, qq_inv[j2] = i1; \
qq_ind[i2] = j1, qq_inv[j1] = i2; \
} while (0)
/* swap columns j1 and j2 of permutation matrix Q */
#define scf_r0_solve _glp_scf_r0_solve
void scf_r0_solve(SCF *scf, int tr, double x[/*1+n0*/]);
/* solve system R0 * x = b or R0'* x = b */
#define scf_s0_solve _glp_scf_s0_solve
void scf_s0_solve(SCF *scf, int tr, double x[/*1+n0*/],
double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]);
/* solve system S0 * x = b or S0'* x = b */
#define scf_r_prod _glp_scf_r_prod
void scf_r_prod(SCF *scf, double y[/*1+nn*/], double a, const double
x[/*1+n0*/]);
/* compute product y := y + alpha * R * x */
#define scf_rt_prod _glp_scf_rt_prod
void scf_rt_prod(SCF *scf, double y[/*1+n0*/], double a, const double
x[/*1+nn*/]);
/* compute product y := y + alpha * R'* x */
#define scf_s_prod _glp_scf_s_prod
void scf_s_prod(SCF *scf, double y[/*1+n0*/], double a, const double
x[/*1+nn*/]);
/* compute product y := y + alpha * S * x */
#define scf_st_prod _glp_scf_st_prod
void scf_st_prod(SCF *scf, double y[/*1+nn*/], double a, const double
x[/*1+n0*/]);
/* compute product y := y + alpha * S'* x */
#define scf_a_solve _glp_scf_a_solve
void scf_a_solve(SCF *scf, double x[/*1+n*/],
double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/],
double work2[/*1+n*/], double work3[/*1+n*/]);
/* solve system A * x = b */
#define scf_at_solve _glp_scf_at_solve
void scf_at_solve(SCF *scf, double x[/*1+n*/],
double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/],
double work2[/*1+n*/], double work3[/*1+n*/]);
/* solve system A'* x = b */
#define scf_add_r_row _glp_scf_add_r_row
void scf_add_r_row(SCF *scf, const double w[/*1+n0*/]);
/* add new row to matrix R */
#define scf_add_s_col _glp_scf_add_s_col
void scf_add_s_col(SCF *scf, const double v[/*1+n0*/]);
/* add new column to matrix S */
#define scf_update_aug _glp_scf_update_aug
int scf_update_aug(SCF *scf, double b[/*1+n0*/], double d[/*1+n0*/],
double f[/*1+nn*/], double g[/*1+nn*/], double h, int upd,
double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]);
/* update factorization of augmented matrix */
#endif
/* eof */
|
