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
|
/* btf.h (sparse block triangular LU-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 BTF_H
#define BTF_H
#include "sva.h"
/***********************************************************************
* The structure BTF describes BT-factorization, which is sparse block
* triangular LU-factorization.
*
* The BT-factorization has the following format:
*
* A = P * A~ * Q, (1)
*
* where A is a given (unsymmetric) square matrix, A~ is an upper block
* triangular matrix (see below), P and Q are permutation matrices. All
* the matrices have the same order n.
*
* The matrix A~, which is a permuted version of the original matrix A,
* has the following structure:
*
* A~[1,1] A~[1,2] ... A~[1,num-1] A~[1,num]
*
* A~[2,2] ... A~[2,num-1] A~[2,num]
*
* . . . . . . . . . (2)
*
* A~[num-1,num-1] A~[num-1,num]
*
* A~[num,num]
*
* where A~[i,j] is a submatrix called a "block," num is the number of
* blocks. Each diagonal block A~[k,k] is a non-singular square matrix,
* and each subdiagonal block A~[i,j], i > j, is a zero submatrix, thus
* A~ is an upper block triangular matrix.
*
* Permutation matrices P and Q are stored in ordinary arrays in both
* row- and column-like formats.
*
* The original matrix A is stored in both row- and column-wise sparse
* formats in the associated sparse vector area (SVA). Should note that
* elements of all diagonal blocks A~[k,k] in matrix A are set to zero
* (i.e. removed), so only elements of non-diagonal blocks are stored.
*
* Each diagonal block A~[k,k], 1 <= k <= num, is stored in the form of
* LU-factorization (see the module LUF). */
typedef struct BTF BTF;
struct BTF
{ /* sparse block triangular LU-factorization */
int n;
/* order of matrices A, A~, P, Q */
SVA *sva;
/* associated sparse vector area used to store rows and columns
* of matrix A as well as sparse vectors for LU-factorizations of
* all diagonal blocks A~[k,k] */
/*--------------------------------------------------------------*/
/* matrix P */
int *pp_ind; /* int pp_ind[1+n]; */
/* pp_ind[i] = j means that P[i,j] = 1 */
int *pp_inv; /* int pp_inv[1+n]; */
/* pp_inv[j] = i means that P[i,j] = 1 */
/* if i-th row of matrix A is i'-th row of matrix A~, then
* pp_ind[i] = i' and pp_inv[i'] = i */
/*--------------------------------------------------------------*/
/* matrix Q */
int *qq_ind; /* int qq_ind[1+n]; */
/* qq_ind[i] = j means that Q[i,j] = 1 */
int *qq_inv; /* int qq_inv[1+n]; */
/* qq_inv[j] = i means that Q[i,j] = 1 */
/* if j-th column of matrix A is j'-th column of matrix A~, then
* qq_ind[j'] = j and qq_inv[j] = j' */
/*--------------------------------------------------------------*/
/* block triangular structure of matrix A~ */
int num;
/* number of diagonal blocks, 1 <= num <= n */
int *beg; /* int beg[1+num+1]; */
/* beg[0] is not used;
* beg[k], 1 <= k <= num, is index of first row/column of k-th
* block of matrix A~;
* beg[num+1] is always n+1;
* note that order (size) of k-th diagonal block can be computed
* as beg[k+1] - beg[k] */
/*--------------------------------------------------------------*/
/* original matrix A in row-wise format */
/* NOTE: elements of all diagonal blocks A~[k,k] are removed */
int ar_ref;
/* reference number of sparse vector in SVA, which is the first
* row of matrix A */
#if 0 + 0
int *ar_ptr = &sva->ptr[ar_ref-1];
/* ar_ptr[0] is not used;
* ar_ptr[i], 1 <= i <= n, is pointer to i-th row in SVA */
int *ar_len = &sva->ptr[ar_ref-1];
/* ar_len[0] is not used;
* ar_len[i], 1 <= i <= n, is length of i-th row */
#endif
/*--------------------------------------------------------------*/
/* original matrix A in column-wise format */
/* NOTE: elements of all diagonal blocks A~[k,k] are removed */
int ac_ref;
/* reference number of sparse vector in SVA, which is the first
* column of matrix A */
#if 0 + 0
int *ac_ptr = &sva->ptr[ac_ref-1];
/* ac_ptr[0] is not used;
* ac_ptr[j], 1 <= j <= n, is pointer to j-th column in SVA */
int *ac_len = &sva->ptr[ac_ref-1];
/* ac_len[0] is not used;
* ac_len[j], 1 <= j <= n, is length of j-th column */
#endif
/*--------------------------------------------------------------*/
/* LU-factorizations of diagonal blocks A~[k,k] */
/* to decrease overhead expenses similar arrays for all LUFs are
* packed into a single array; for example, elements fr_ptr[1],
* ..., fr_ptr[n1], where n1 = beg[2] - beg[1], are related to
* LUF for first diagonal block A~[1,1], elements fr_ptr[n1+1],
* ..., fr_ptr[n1+n2], where n2 = beg[3] - beg[2], are related to
* LUF for second diagonal block A~[2,2], etc.; in other words,
* elements related to LUF for k-th diagonal block A~[k,k] have
* indices beg[k], beg[k]+1, ..., beg[k+1]-1 */
/* for details about LUF see description of the LUF module */
int fr_ref;
/* reference number of sparse vector in SVA, which is the first
row of matrix F for first diagonal block A~[1,1] */
int fc_ref;
/* reference number of sparse vector in SVA, which is the first
column of matrix F for first diagonal block A~[1,1] */
int vr_ref;
/* reference number of sparse vector in SVA, which is the first
row of matrix V for first diagonal block A~[1,1] */
double *vr_piv; /* double vr_piv[1+n]; */
/* vr_piv[0] is not used;
vr_piv[1,...,n] are pivot elements for all diagonal blocks */
int vc_ref;
/* reference number of sparse vector in SVA, which is the first
column of matrix V for first diagonal block A~[1,1] */
int *p1_ind; /* int p1_ind[1+n]; */
int *p1_inv; /* int p1_inv[1+n]; */
int *q1_ind; /* int q1_ind[1+n]; */
int *q1_inv; /* int q1_inv[1+n]; */
/* permutation matrices P and Q for all diagonal blocks */
};
#define btf_store_a_cols _glp_btf_store_a_cols
int btf_store_a_cols(BTF *btf, int (*col)(void *info, int j, int ind[],
double val[]), void *info, int ind[], double val[]);
/* store pattern of matrix A in column-wise format */
#define btf_make_blocks _glp_btf_make_blocks
int btf_make_blocks(BTF *btf);
/* permutations to block triangular form */
#define btf_check_blocks _glp_btf_check_blocks
void btf_check_blocks(BTF *btf);
/* check structure of matrix A~ */
#define btf_build_a_rows _glp_btf_build_a_rows
void btf_build_a_rows(BTF *btf, int len[/*1+n*/]);
/* build matrix A in row-wise format */
#define btf_a_solve _glp_btf_a_solve
void btf_a_solve(BTF *btf, double b[/*1+n*/], double x[/*1+n*/],
double w1[/*1+n*/], double w2[/*1+n*/]);
/* solve system A * x = b */
#define btf_at_solve _glp_btf_at_solve
void btf_at_solve(BTF *btf, double b[/*1+n*/], double x[/*1+n*/],
double w1[/*1+n*/], double w2[/*1+n*/]);
/* solve system A'* x = b */
#define btf_at_solve1 _glp_btf_at_solve1
void btf_at_solve1(BTF *btf, double e[/*1+n*/], double y[/*1+n*/],
double w1[/*1+n*/], double w2[/*1+n*/]);
/* solve system A'* y = e' to cause growth in y */
#define btf_estimate_norm _glp_btf_estimate_norm
double btf_estimate_norm(BTF *btf, double w1[/*1+n*/], double
w2[/*1+n*/], double w3[/*1+n*/], double w4[/*1+n*/]);
/* estimate 1-norm of inv(A) */
#endif
/* eof */
|
