rev |
line source |
alpar@9
|
1 /* glpmat.h (linear algebra routines) */
|
alpar@9
|
2
|
alpar@9
|
3 /***********************************************************************
|
alpar@9
|
4 * This code is part of GLPK (GNU Linear Programming Kit).
|
alpar@9
|
5 *
|
alpar@9
|
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
|
alpar@9
|
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
|
alpar@9
|
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
|
alpar@9
|
9 * E-mail: <mao@gnu.org>.
|
alpar@9
|
10 *
|
alpar@9
|
11 * GLPK is free software: you can redistribute it and/or modify it
|
alpar@9
|
12 * under the terms of the GNU General Public License as published by
|
alpar@9
|
13 * the Free Software Foundation, either version 3 of the License, or
|
alpar@9
|
14 * (at your option) any later version.
|
alpar@9
|
15 *
|
alpar@9
|
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
|
alpar@9
|
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
alpar@9
|
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
|
alpar@9
|
19 * License for more details.
|
alpar@9
|
20 *
|
alpar@9
|
21 * You should have received a copy of the GNU General Public License
|
alpar@9
|
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
|
alpar@9
|
23 ***********************************************************************/
|
alpar@9
|
24
|
alpar@9
|
25 #ifndef GLPMAT_H
|
alpar@9
|
26 #define GLPMAT_H
|
alpar@9
|
27
|
alpar@9
|
28 /***********************************************************************
|
alpar@9
|
29 * FULL-VECTOR STORAGE
|
alpar@9
|
30 *
|
alpar@9
|
31 * For a sparse vector x having n elements, ne of which are non-zero,
|
alpar@9
|
32 * the full-vector storage format uses two arrays x_ind and x_vec, which
|
alpar@9
|
33 * are set up as follows:
|
alpar@9
|
34 *
|
alpar@9
|
35 * x_ind is an integer array of length [1+ne]. Location x_ind[0] is
|
alpar@9
|
36 * not used, and locations x_ind[1], ..., x_ind[ne] contain indices of
|
alpar@9
|
37 * non-zero elements in vector x.
|
alpar@9
|
38 *
|
alpar@9
|
39 * x_vec is a floating-point array of length [1+n]. Location x_vec[0]
|
alpar@9
|
40 * is not used, and locations x_vec[1], ..., x_vec[n] contain numeric
|
alpar@9
|
41 * values of ALL elements in vector x, including its zero elements.
|
alpar@9
|
42 *
|
alpar@9
|
43 * Let, for example, the following sparse vector x be given:
|
alpar@9
|
44 *
|
alpar@9
|
45 * (0, 1, 0, 0, 2, 3, 0, 4)
|
alpar@9
|
46 *
|
alpar@9
|
47 * Then the arrays are:
|
alpar@9
|
48 *
|
alpar@9
|
49 * x_ind = { X; 2, 5, 6, 8 }
|
alpar@9
|
50 *
|
alpar@9
|
51 * x_vec = { X; 0, 1, 0, 0, 2, 3, 0, 4 }
|
alpar@9
|
52 *
|
alpar@9
|
53 * COMPRESSED-VECTOR STORAGE
|
alpar@9
|
54 *
|
alpar@9
|
55 * For a sparse vector x having n elements, ne of which are non-zero,
|
alpar@9
|
56 * the compressed-vector storage format uses two arrays x_ind and x_vec,
|
alpar@9
|
57 * which are set up as follows:
|
alpar@9
|
58 *
|
alpar@9
|
59 * x_ind is an integer array of length [1+ne]. Location x_ind[0] is
|
alpar@9
|
60 * not used, and locations x_ind[1], ..., x_ind[ne] contain indices of
|
alpar@9
|
61 * non-zero elements in vector x.
|
alpar@9
|
62 *
|
alpar@9
|
63 * x_vec is a floating-point array of length [1+ne]. Location x_vec[0]
|
alpar@9
|
64 * is not used, and locations x_vec[1], ..., x_vec[ne] contain numeric
|
alpar@9
|
65 * values of corresponding non-zero elements in vector x.
|
alpar@9
|
66 *
|
alpar@9
|
67 * Let, for example, the following sparse vector x be given:
|
alpar@9
|
68 *
|
alpar@9
|
69 * (0, 1, 0, 0, 2, 3, 0, 4)
|
alpar@9
|
70 *
|
alpar@9
|
71 * Then the arrays are:
|
alpar@9
|
72 *
|
alpar@9
|
73 * x_ind = { X; 2, 5, 6, 8 }
|
alpar@9
|
74 *
|
alpar@9
|
75 * x_vec = { X; 1, 2, 3, 4 }
|
alpar@9
|
76 *
|
alpar@9
|
77 * STORAGE-BY-ROWS
|
alpar@9
|
78 *
|
alpar@9
|
79 * For a sparse matrix A, which has m rows, n columns, and ne non-zero
|
alpar@9
|
80 * elements the storage-by-rows format uses three arrays A_ptr, A_ind,
|
alpar@9
|
81 * and A_val, which are set up as follows:
|
alpar@9
|
82 *
|
alpar@9
|
83 * A_ptr is an integer array of length [1+m+1] also called "row pointer
|
alpar@9
|
84 * array". It contains the relative starting positions of each row of A
|
alpar@9
|
85 * in the arrays A_ind and A_val, i.e. element A_ptr[i], 1 <= i <= m,
|
alpar@9
|
86 * indicates where row i begins in the arrays A_ind and A_val. If all
|
alpar@9
|
87 * elements in row i are zero, then A_ptr[i] = A_ptr[i+1]. Location
|
alpar@9
|
88 * A_ptr[0] is not used, location A_ptr[1] must contain 1, and location
|
alpar@9
|
89 * A_ptr[m+1] must contain ne+1 that indicates the position after the
|
alpar@9
|
90 * last element in the arrays A_ind and A_val.
|
alpar@9
|
91 *
|
alpar@9
|
92 * A_ind is an integer array of length [1+ne]. Location A_ind[0] is not
|
alpar@9
|
93 * used, and locations A_ind[1], ..., A_ind[ne] contain column indices
|
alpar@9
|
94 * of (non-zero) elements in matrix A.
|
alpar@9
|
95 *
|
alpar@9
|
96 * A_val is a floating-point array of length [1+ne]. Location A_val[0]
|
alpar@9
|
97 * is not used, and locations A_val[1], ..., A_val[ne] contain numeric
|
alpar@9
|
98 * values of non-zero elements in matrix A.
|
alpar@9
|
99 *
|
alpar@9
|
100 * Non-zero elements of matrix A are stored contiguously, and the rows
|
alpar@9
|
101 * of matrix A are stored consecutively from 1 to m in the arrays A_ind
|
alpar@9
|
102 * and A_val. The elements in each row of A may be stored in any order
|
alpar@9
|
103 * in A_ind and A_val. Note that elements with duplicate column indices
|
alpar@9
|
104 * are not allowed.
|
alpar@9
|
105 *
|
alpar@9
|
106 * Let, for example, the following sparse matrix A be given:
|
alpar@9
|
107 *
|
alpar@9
|
108 * | 11 . 13 . . . |
|
alpar@9
|
109 * | 21 22 . 24 . . |
|
alpar@9
|
110 * | . 32 33 . . . |
|
alpar@9
|
111 * | . . 43 44 . 46 |
|
alpar@9
|
112 * | . . . . . . |
|
alpar@9
|
113 * | 61 62 . . . 66 |
|
alpar@9
|
114 *
|
alpar@9
|
115 * Then the arrays are:
|
alpar@9
|
116 *
|
alpar@9
|
117 * A_ptr = { X; 1, 3, 6, 8, 11, 11; 14 }
|
alpar@9
|
118 *
|
alpar@9
|
119 * A_ind = { X; 1, 3; 4, 2, 1; 2, 3; 4, 3, 6; 1, 2, 6 }
|
alpar@9
|
120 *
|
alpar@9
|
121 * A_val = { X; 11, 13; 24, 22, 21; 32, 33; 44, 43, 46; 61, 62, 66 }
|
alpar@9
|
122 *
|
alpar@9
|
123 * PERMUTATION MATRICES
|
alpar@9
|
124 *
|
alpar@9
|
125 * Let P be a permutation matrix of the order n. It is represented as
|
alpar@9
|
126 * an integer array P_per of length [1+n+n] as follows: if p[i,j] = 1,
|
alpar@9
|
127 * then P_per[i] = j and P_per[n+j] = i. Location P_per[0] is not used.
|
alpar@9
|
128 *
|
alpar@9
|
129 * Let A' = P*A. If i-th row of A corresponds to i'-th row of A', then
|
alpar@9
|
130 * P_per[i'] = i and P_per[n+i] = i'.
|
alpar@9
|
131 *
|
alpar@9
|
132 * References:
|
alpar@9
|
133 *
|
alpar@9
|
134 * 1. Gustavson F.G. Some basic techniques for solving sparse systems of
|
alpar@9
|
135 * linear equations. In Rose and Willoughby (1972), pp. 41-52.
|
alpar@9
|
136 *
|
alpar@9
|
137 * 2. Basic Linear Algebra Subprograms Technical (BLAST) Forum Standard.
|
alpar@9
|
138 * University of Tennessee (2001). */
|
alpar@9
|
139
|
alpar@9
|
140 #define check_fvs _glp_mat_check_fvs
|
alpar@9
|
141 int check_fvs(int n, int nnz, int ind[], double vec[]);
|
alpar@9
|
142 /* check sparse vector in full-vector storage format */
|
alpar@9
|
143
|
alpar@9
|
144 #define check_pattern _glp_mat_check_pattern
|
alpar@9
|
145 int check_pattern(int m, int n, int A_ptr[], int A_ind[]);
|
alpar@9
|
146 /* check pattern of sparse matrix */
|
alpar@9
|
147
|
alpar@9
|
148 #define transpose _glp_mat_transpose
|
alpar@9
|
149 void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[],
|
alpar@9
|
150 int AT_ptr[], int AT_ind[], double AT_val[]);
|
alpar@9
|
151 /* transpose sparse matrix */
|
alpar@9
|
152
|
alpar@9
|
153 #define adat_symbolic _glp_mat_adat_symbolic
|
alpar@9
|
154 int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[],
|
alpar@9
|
155 int S_ptr[]);
|
alpar@9
|
156 /* compute S = P*A*D*A'*P' (symbolic phase) */
|
alpar@9
|
157
|
alpar@9
|
158 #define adat_numeric _glp_mat_adat_numeric
|
alpar@9
|
159 void adat_numeric(int m, int n, int P_per[],
|
alpar@9
|
160 int A_ptr[], int A_ind[], double A_val[], double D_diag[],
|
alpar@9
|
161 int S_ptr[], int S_ind[], double S_val[], double S_diag[]);
|
alpar@9
|
162 /* compute S = P*A*D*A'*P' (numeric phase) */
|
alpar@9
|
163
|
alpar@9
|
164 #define min_degree _glp_mat_min_degree
|
alpar@9
|
165 void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]);
|
alpar@9
|
166 /* minimum degree ordering */
|
alpar@9
|
167
|
alpar@9
|
168 #define amd_order1 _glp_mat_amd_order1
|
alpar@9
|
169 void amd_order1(int n, int A_ptr[], int A_ind[], int P_per[]);
|
alpar@9
|
170 /* approximate minimum degree ordering (AMD) */
|
alpar@9
|
171
|
alpar@9
|
172 #define symamd_ord _glp_mat_symamd_ord
|
alpar@9
|
173 void symamd_ord(int n, int A_ptr[], int A_ind[], int P_per[]);
|
alpar@9
|
174 /* approximate minimum degree ordering (SYMAMD) */
|
alpar@9
|
175
|
alpar@9
|
176 #define chol_symbolic _glp_mat_chol_symbolic
|
alpar@9
|
177 int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]);
|
alpar@9
|
178 /* compute Cholesky factorization (symbolic phase) */
|
alpar@9
|
179
|
alpar@9
|
180 #define chol_numeric _glp_mat_chol_numeric
|
alpar@9
|
181 int chol_numeric(int n,
|
alpar@9
|
182 int A_ptr[], int A_ind[], double A_val[], double A_diag[],
|
alpar@9
|
183 int U_ptr[], int U_ind[], double U_val[], double U_diag[]);
|
alpar@9
|
184 /* compute Cholesky factorization (numeric phase) */
|
alpar@9
|
185
|
alpar@9
|
186 #define u_solve _glp_mat_u_solve
|
alpar@9
|
187 void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
|
alpar@9
|
188 double U_diag[], double x[]);
|
alpar@9
|
189 /* solve upper triangular system U*x = b */
|
alpar@9
|
190
|
alpar@9
|
191 #define ut_solve _glp_mat_ut_solve
|
alpar@9
|
192 void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
|
alpar@9
|
193 double U_diag[], double x[]);
|
alpar@9
|
194 /* solve lower triangular system U'*x = b */
|
alpar@9
|
195
|
alpar@9
|
196 #endif
|
alpar@9
|
197
|
alpar@9
|
198 /* eof */
|