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
|
/* glpmat.h (linear algebra routines) */
/***********************************************************************
* 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 GLPMAT_H
#define GLPMAT_H
/***********************************************************************
* FULL-VECTOR STORAGE
*
* For a sparse vector x having n elements, ne of which are non-zero,
* the full-vector storage format uses two arrays x_ind and x_vec, which
* are set up as follows:
*
* x_ind is an integer array of length [1+ne]. Location x_ind[0] is
* not used, and locations x_ind[1], ..., x_ind[ne] contain indices of
* non-zero elements in vector x.
*
* x_vec is a floating-point array of length [1+n]. Location x_vec[0]
* is not used, and locations x_vec[1], ..., x_vec[n] contain numeric
* values of ALL elements in vector x, including its zero elements.
*
* Let, for example, the following sparse vector x be given:
*
* (0, 1, 0, 0, 2, 3, 0, 4)
*
* Then the arrays are:
*
* x_ind = { X; 2, 5, 6, 8 }
*
* x_vec = { X; 0, 1, 0, 0, 2, 3, 0, 4 }
*
* COMPRESSED-VECTOR STORAGE
*
* For a sparse vector x having n elements, ne of which are non-zero,
* the compressed-vector storage format uses two arrays x_ind and x_vec,
* which are set up as follows:
*
* x_ind is an integer array of length [1+ne]. Location x_ind[0] is
* not used, and locations x_ind[1], ..., x_ind[ne] contain indices of
* non-zero elements in vector x.
*
* x_vec is a floating-point array of length [1+ne]. Location x_vec[0]
* is not used, and locations x_vec[1], ..., x_vec[ne] contain numeric
* values of corresponding non-zero elements in vector x.
*
* Let, for example, the following sparse vector x be given:
*
* (0, 1, 0, 0, 2, 3, 0, 4)
*
* Then the arrays are:
*
* x_ind = { X; 2, 5, 6, 8 }
*
* x_vec = { X; 1, 2, 3, 4 }
*
* STORAGE-BY-ROWS
*
* For a sparse matrix A, which has m rows, n columns, and ne non-zero
* elements the storage-by-rows format uses three arrays A_ptr, A_ind,
* and A_val, which are set up as follows:
*
* A_ptr is an integer array of length [1+m+1] also called "row pointer
* array". It contains the relative starting positions of each row of A
* in the arrays A_ind and A_val, i.e. element A_ptr[i], 1 <= i <= m,
* indicates where row i begins in the arrays A_ind and A_val. If all
* elements in row i are zero, then A_ptr[i] = A_ptr[i+1]. Location
* A_ptr[0] is not used, location A_ptr[1] must contain 1, and location
* A_ptr[m+1] must contain ne+1 that indicates the position after the
* last element in the arrays A_ind and A_val.
*
* A_ind is an integer array of length [1+ne]. Location A_ind[0] is not
* used, and locations A_ind[1], ..., A_ind[ne] contain column indices
* of (non-zero) elements in matrix A.
*
* A_val is a floating-point array of length [1+ne]. Location A_val[0]
* is not used, and locations A_val[1], ..., A_val[ne] contain numeric
* values of non-zero elements in matrix A.
*
* Non-zero elements of matrix A are stored contiguously, and the rows
* of matrix A are stored consecutively from 1 to m in the arrays A_ind
* and A_val. The elements in each row of A may be stored in any order
* in A_ind and A_val. Note that elements with duplicate column indices
* are not allowed.
*
* Let, for example, the following sparse matrix A be given:
*
* | 11 . 13 . . . |
* | 21 22 . 24 . . |
* | . 32 33 . . . |
* | . . 43 44 . 46 |
* | . . . . . . |
* | 61 62 . . . 66 |
*
* Then the arrays are:
*
* A_ptr = { X; 1, 3, 6, 8, 11, 11; 14 }
*
* A_ind = { X; 1, 3; 4, 2, 1; 2, 3; 4, 3, 6; 1, 2, 6 }
*
* A_val = { X; 11, 13; 24, 22, 21; 32, 33; 44, 43, 46; 61, 62, 66 }
*
* PERMUTATION MATRICES
*
* Let P be a permutation matrix of the order n. It is represented as
* an integer array P_per of length [1+n+n] as follows: if p[i,j] = 1,
* then P_per[i] = j and P_per[n+j] = i. Location P_per[0] is not used.
*
* Let A' = P*A. If i-th row of A corresponds to i'-th row of A', then
* P_per[i'] = i and P_per[n+i] = i'.
*
* References:
*
* 1. Gustavson F.G. Some basic techniques for solving sparse systems of
* linear equations. In Rose and Willoughby (1972), pp. 41-52.
*
* 2. Basic Linear Algebra Subprograms Technical (BLAST) Forum Standard.
* University of Tennessee (2001). */
#define check_fvs _glp_mat_check_fvs
int check_fvs(int n, int nnz, int ind[], double vec[]);
/* check sparse vector in full-vector storage format */
#define check_pattern _glp_mat_check_pattern
int check_pattern(int m, int n, int A_ptr[], int A_ind[]);
/* check pattern of sparse matrix */
#define transpose _glp_mat_transpose
void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[],
int AT_ptr[], int AT_ind[], double AT_val[]);
/* transpose sparse matrix */
#define adat_symbolic _glp_mat_adat_symbolic
int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[],
int S_ptr[]);
/* compute S = P*A*D*A'*P' (symbolic phase) */
#define adat_numeric _glp_mat_adat_numeric
void adat_numeric(int m, int n, int P_per[],
int A_ptr[], int A_ind[], double A_val[], double D_diag[],
int S_ptr[], int S_ind[], double S_val[], double S_diag[]);
/* compute S = P*A*D*A'*P' (numeric phase) */
#define min_degree _glp_mat_min_degree
void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]);
/* minimum degree ordering */
#define amd_order1 _glp_mat_amd_order1
void amd_order1(int n, int A_ptr[], int A_ind[], int P_per[]);
/* approximate minimum degree ordering (AMD) */
#define symamd_ord _glp_mat_symamd_ord
void symamd_ord(int n, int A_ptr[], int A_ind[], int P_per[]);
/* approximate minimum degree ordering (SYMAMD) */
#define chol_symbolic _glp_mat_chol_symbolic
int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]);
/* compute Cholesky factorization (symbolic phase) */
#define chol_numeric _glp_mat_chol_numeric
int chol_numeric(int n,
int A_ptr[], int A_ind[], double A_val[], double A_diag[],
int U_ptr[], int U_ind[], double U_val[], double U_diag[]);
/* compute Cholesky factorization (numeric phase) */
#define u_solve _glp_mat_u_solve
void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
double U_diag[], double x[]);
/* solve upper triangular system U*x = b */
#define ut_solve _glp_mat_ut_solve
void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
double U_diag[], double x[]);
/* solve lower triangular system U'*x = b */
#endif
/* eof */
|
