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
|
/* fhv.h (sparse updatable FHV-factorization) */
/***********************************************************************
* 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 FHV_H
#define FHV_H
#include "luf.h"
/***********************************************************************
* The structure FHV describes sparse updatable FHV-factorization.
*
* The FHV-factorization has the following format:
*
* A = F * H * V, (1)
*
* F = P0 * L * P0', (2)
*
* H = H[1] * H[2] * ... * H[nfs], (3)
*
* V = P * U * Q, (4)
*
* where: A is a given (unsymmetric) square matrix; F, H, V are matrix
* factors actually computed; L is a lower triangular matrix with unity
* diagonal; U is an upper tringular matrix; H[k], k = 1, 2, ..., nfs,
* is a row-like factor, which differs from unity matrix only in one
* row called a non-trivial row; P0, P, Q are permutation matrices; and
* P0' is a matrix transposed to P0.
*
* Matrices F, V, P, Q are stored in the underlying LUF object.
*
* Non-trivial rows of factors H[k] are stored as sparse vectors in the
* right (static) part of the sparse vector area (SVA). Note that unity
* diagonal elements of non-trivial rows are not stored.
*
* Matrix P0 is stored in the same way as matrix P.
*
* Matrices L and U are completely defined by matrices F, V, P, and Q,
* and therefore not stored explicitly. */
typedef struct FHV FHV;
struct FHV
{ /* FHV-factorization */
LUF *luf;
/* LU-factorization (contains matrices F, V, P, Q) */
/*--------------------------------------------------------------*/
/* matrix H in the form of eta file */
int nfs_max;
/* maximal number of row-like factors (this limits the number of
* updates of the factorization) */
int nfs;
/* current number of row-like factors, 0 <= nfs <= nfs_max */
int *hh_ind; /* int hh_ind[1+nfs_max]; */
/* hh_ind[0] is not used;
* hh_ind[k], 1 <= k <= nfs, is number of non-trivial row of
* factor H[k] */
int hh_ref;
/* reference number of sparse vector in SVA, which is non-trivial
* row of factor H[1] */
#if 0 + 0
int *hh_ptr = &sva->ptr[hh_ref-1];
/* hh_ptr[0] is not used;
* hh_ptr[k], 1 <= k <= nfs, is pointer to non-trivial row of
* factor H[k] */
int *hh_len = &sva->len[hh_ref-1];
/* hh_len[0] is not used;
* hh_len[k], 1 <= k <= nfs, is number of non-zero elements in
* non-trivial row of factor H[k] */
#endif
/*--------------------------------------------------------------*/
/* matrix P0 */
int *p0_ind; /* int p0_ind[1+n]; */
/* p0_ind[i] = j means that P0[i,j] = 1 */
int *p0_inv; /* int p0_inv[1+n]; */
/* p0_inv[j] = i means that P0[i,j] = 1 */
};
#define fhv_ft_update _glp_fhv_ft_update
int fhv_ft_update(FHV *fhv, int q, int aq_len, const int aq_ind[],
const double aq_val[], int ind[/*1+n*/], double val[/*1+n*/],
double work[/*1+n*/]);
/* update FHV-factorization (Forrest-Tomlin) */
#define fhv_h_solve _glp_fhv_h_solve
void fhv_h_solve(FHV *fhv, double x[/*1+n*/]);
/* solve system H * x = b */
#define fhv_ht_solve _glp_fhv_ht_solve
void fhv_ht_solve(FHV *fhv, double x[/*1+n*/]);
/* solve system H' * x = b */
#endif
/* eof */
|
