1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpmat.h Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,198 @@
1.4 +/* glpmat.h (linear algebra routines) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#ifndef GLPMAT_H
1.29 +#define GLPMAT_H
1.30 +
1.31 +/***********************************************************************
1.32 +* FULL-VECTOR STORAGE
1.33 +*
1.34 +* For a sparse vector x having n elements, ne of which are non-zero,
1.35 +* the full-vector storage format uses two arrays x_ind and x_vec, which
1.36 +* are set up as follows:
1.37 +*
1.38 +* x_ind is an integer array of length [1+ne]. Location x_ind[0] is
1.39 +* not used, and locations x_ind[1], ..., x_ind[ne] contain indices of
1.40 +* non-zero elements in vector x.
1.41 +*
1.42 +* x_vec is a floating-point array of length [1+n]. Location x_vec[0]
1.43 +* is not used, and locations x_vec[1], ..., x_vec[n] contain numeric
1.44 +* values of ALL elements in vector x, including its zero elements.
1.45 +*
1.46 +* Let, for example, the following sparse vector x be given:
1.47 +*
1.48 +* (0, 1, 0, 0, 2, 3, 0, 4)
1.49 +*
1.50 +* Then the arrays are:
1.51 +*
1.52 +* x_ind = { X; 2, 5, 6, 8 }
1.53 +*
1.54 +* x_vec = { X; 0, 1, 0, 0, 2, 3, 0, 4 }
1.55 +*
1.56 +* COMPRESSED-VECTOR STORAGE
1.57 +*
1.58 +* For a sparse vector x having n elements, ne of which are non-zero,
1.59 +* the compressed-vector storage format uses two arrays x_ind and x_vec,
1.60 +* which are set up as follows:
1.61 +*
1.62 +* x_ind is an integer array of length [1+ne]. Location x_ind[0] is
1.63 +* not used, and locations x_ind[1], ..., x_ind[ne] contain indices of
1.64 +* non-zero elements in vector x.
1.65 +*
1.66 +* x_vec is a floating-point array of length [1+ne]. Location x_vec[0]
1.67 +* is not used, and locations x_vec[1], ..., x_vec[ne] contain numeric
1.68 +* values of corresponding non-zero elements in vector x.
1.69 +*
1.70 +* Let, for example, the following sparse vector x be given:
1.71 +*
1.72 +* (0, 1, 0, 0, 2, 3, 0, 4)
1.73 +*
1.74 +* Then the arrays are:
1.75 +*
1.76 +* x_ind = { X; 2, 5, 6, 8 }
1.77 +*
1.78 +* x_vec = { X; 1, 2, 3, 4 }
1.79 +*
1.80 +* STORAGE-BY-ROWS
1.81 +*
1.82 +* For a sparse matrix A, which has m rows, n columns, and ne non-zero
1.83 +* elements the storage-by-rows format uses three arrays A_ptr, A_ind,
1.84 +* and A_val, which are set up as follows:
1.85 +*
1.86 +* A_ptr is an integer array of length [1+m+1] also called "row pointer
1.87 +* array". It contains the relative starting positions of each row of A
1.88 +* in the arrays A_ind and A_val, i.e. element A_ptr[i], 1 <= i <= m,
1.89 +* indicates where row i begins in the arrays A_ind and A_val. If all
1.90 +* elements in row i are zero, then A_ptr[i] = A_ptr[i+1]. Location
1.91 +* A_ptr[0] is not used, location A_ptr[1] must contain 1, and location
1.92 +* A_ptr[m+1] must contain ne+1 that indicates the position after the
1.93 +* last element in the arrays A_ind and A_val.
1.94 +*
1.95 +* A_ind is an integer array of length [1+ne]. Location A_ind[0] is not
1.96 +* used, and locations A_ind[1], ..., A_ind[ne] contain column indices
1.97 +* of (non-zero) elements in matrix A.
1.98 +*
1.99 +* A_val is a floating-point array of length [1+ne]. Location A_val[0]
1.100 +* is not used, and locations A_val[1], ..., A_val[ne] contain numeric
1.101 +* values of non-zero elements in matrix A.
1.102 +*
1.103 +* Non-zero elements of matrix A are stored contiguously, and the rows
1.104 +* of matrix A are stored consecutively from 1 to m in the arrays A_ind
1.105 +* and A_val. The elements in each row of A may be stored in any order
1.106 +* in A_ind and A_val. Note that elements with duplicate column indices
1.107 +* are not allowed.
1.108 +*
1.109 +* Let, for example, the following sparse matrix A be given:
1.110 +*
1.111 +* | 11 . 13 . . . |
1.112 +* | 21 22 . 24 . . |
1.113 +* | . 32 33 . . . |
1.114 +* | . . 43 44 . 46 |
1.115 +* | . . . . . . |
1.116 +* | 61 62 . . . 66 |
1.117 +*
1.118 +* Then the arrays are:
1.119 +*
1.120 +* A_ptr = { X; 1, 3, 6, 8, 11, 11; 14 }
1.121 +*
1.122 +* A_ind = { X; 1, 3; 4, 2, 1; 2, 3; 4, 3, 6; 1, 2, 6 }
1.123 +*
1.124 +* A_val = { X; 11, 13; 24, 22, 21; 32, 33; 44, 43, 46; 61, 62, 66 }
1.125 +*
1.126 +* PERMUTATION MATRICES
1.127 +*
1.128 +* Let P be a permutation matrix of the order n. It is represented as
1.129 +* an integer array P_per of length [1+n+n] as follows: if p[i,j] = 1,
1.130 +* then P_per[i] = j and P_per[n+j] = i. Location P_per[0] is not used.
1.131 +*
1.132 +* Let A' = P*A. If i-th row of A corresponds to i'-th row of A', then
1.133 +* P_per[i'] = i and P_per[n+i] = i'.
1.134 +*
1.135 +* References:
1.136 +*
1.137 +* 1. Gustavson F.G. Some basic techniques for solving sparse systems of
1.138 +* linear equations. In Rose and Willoughby (1972), pp. 41-52.
1.139 +*
1.140 +* 2. Basic Linear Algebra Subprograms Technical (BLAST) Forum Standard.
1.141 +* University of Tennessee (2001). */
1.142 +
1.143 +#define check_fvs _glp_mat_check_fvs
1.144 +int check_fvs(int n, int nnz, int ind[], double vec[]);
1.145 +/* check sparse vector in full-vector storage format */
1.146 +
1.147 +#define check_pattern _glp_mat_check_pattern
1.148 +int check_pattern(int m, int n, int A_ptr[], int A_ind[]);
1.149 +/* check pattern of sparse matrix */
1.150 +
1.151 +#define transpose _glp_mat_transpose
1.152 +void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[],
1.153 + int AT_ptr[], int AT_ind[], double AT_val[]);
1.154 +/* transpose sparse matrix */
1.155 +
1.156 +#define adat_symbolic _glp_mat_adat_symbolic
1.157 +int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[],
1.158 + int S_ptr[]);
1.159 +/* compute S = P*A*D*A'*P' (symbolic phase) */
1.160 +
1.161 +#define adat_numeric _glp_mat_adat_numeric
1.162 +void adat_numeric(int m, int n, int P_per[],
1.163 + int A_ptr[], int A_ind[], double A_val[], double D_diag[],
1.164 + int S_ptr[], int S_ind[], double S_val[], double S_diag[]);
1.165 +/* compute S = P*A*D*A'*P' (numeric phase) */
1.166 +
1.167 +#define min_degree _glp_mat_min_degree
1.168 +void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]);
1.169 +/* minimum degree ordering */
1.170 +
1.171 +#define amd_order1 _glp_mat_amd_order1
1.172 +void amd_order1(int n, int A_ptr[], int A_ind[], int P_per[]);
1.173 +/* approximate minimum degree ordering (AMD) */
1.174 +
1.175 +#define symamd_ord _glp_mat_symamd_ord
1.176 +void symamd_ord(int n, int A_ptr[], int A_ind[], int P_per[]);
1.177 +/* approximate minimum degree ordering (SYMAMD) */
1.178 +
1.179 +#define chol_symbolic _glp_mat_chol_symbolic
1.180 +int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]);
1.181 +/* compute Cholesky factorization (symbolic phase) */
1.182 +
1.183 +#define chol_numeric _glp_mat_chol_numeric
1.184 +int chol_numeric(int n,
1.185 + int A_ptr[], int A_ind[], double A_val[], double A_diag[],
1.186 + int U_ptr[], int U_ind[], double U_val[], double U_diag[]);
1.187 +/* compute Cholesky factorization (numeric phase) */
1.188 +
1.189 +#define u_solve _glp_mat_u_solve
1.190 +void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
1.191 + double U_diag[], double x[]);
1.192 +/* solve upper triangular system U*x = b */
1.193 +
1.194 +#define ut_solve _glp_mat_ut_solve
1.195 +void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
1.196 + double U_diag[], double x[]);
1.197 +/* solve lower triangular system U'*x = b */
1.198 +
1.199 +#endif
1.200 +
1.201 +/* eof */