1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpmat.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,924 @@
1.4 +/* glpmat.c */
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 +#include "glpenv.h"
1.29 +#include "glpmat.h"
1.30 +#include "glpqmd.h"
1.31 +#include "amd/amd.h"
1.32 +#include "colamd/colamd.h"
1.33 +
1.34 +/*----------------------------------------------------------------------
1.35 +-- check_fvs - check sparse vector in full-vector storage format.
1.36 +--
1.37 +-- SYNOPSIS
1.38 +--
1.39 +-- #include "glpmat.h"
1.40 +-- int check_fvs(int n, int nnz, int ind[], double vec[]);
1.41 +--
1.42 +-- DESCRIPTION
1.43 +--
1.44 +-- The routine check_fvs checks if a given vector of dimension n in
1.45 +-- full-vector storage format has correct representation.
1.46 +--
1.47 +-- RETURNS
1.48 +--
1.49 +-- The routine returns one of the following codes:
1.50 +--
1.51 +-- 0 - the vector is correct;
1.52 +-- 1 - the number of elements (n) is negative;
1.53 +-- 2 - the number of non-zero elements (nnz) is negative;
1.54 +-- 3 - some element index is out of range;
1.55 +-- 4 - some element index is duplicate;
1.56 +-- 5 - some non-zero element is out of pattern. */
1.57 +
1.58 +int check_fvs(int n, int nnz, int ind[], double vec[])
1.59 +{ int i, t, ret, *flag = NULL;
1.60 + /* check the number of elements */
1.61 + if (n < 0)
1.62 + { ret = 1;
1.63 + goto done;
1.64 + }
1.65 + /* check the number of non-zero elements */
1.66 + if (nnz < 0)
1.67 + { ret = 2;
1.68 + goto done;
1.69 + }
1.70 + /* check vector indices */
1.71 + flag = xcalloc(1+n, sizeof(int));
1.72 + for (i = 1; i <= n; i++) flag[i] = 0;
1.73 + for (t = 1; t <= nnz; t++)
1.74 + { i = ind[t];
1.75 + if (!(1 <= i && i <= n))
1.76 + { ret = 3;
1.77 + goto done;
1.78 + }
1.79 + if (flag[i])
1.80 + { ret = 4;
1.81 + goto done;
1.82 + }
1.83 + flag[i] = 1;
1.84 + }
1.85 + /* check vector elements */
1.86 + for (i = 1; i <= n; i++)
1.87 + { if (!flag[i] && vec[i] != 0.0)
1.88 + { ret = 5;
1.89 + goto done;
1.90 + }
1.91 + }
1.92 + /* the vector is ok */
1.93 + ret = 0;
1.94 +done: if (flag != NULL) xfree(flag);
1.95 + return ret;
1.96 +}
1.97 +
1.98 +/*----------------------------------------------------------------------
1.99 +-- check_pattern - check pattern of sparse matrix.
1.100 +--
1.101 +-- SYNOPSIS
1.102 +--
1.103 +-- #include "glpmat.h"
1.104 +-- int check_pattern(int m, int n, int A_ptr[], int A_ind[]);
1.105 +--
1.106 +-- DESCRIPTION
1.107 +--
1.108 +-- The routine check_pattern checks the pattern of a given mxn matrix
1.109 +-- in storage-by-rows format.
1.110 +--
1.111 +-- RETURNS
1.112 +--
1.113 +-- The routine returns one of the following codes:
1.114 +--
1.115 +-- 0 - the pattern is correct;
1.116 +-- 1 - the number of rows (m) is negative;
1.117 +-- 2 - the number of columns (n) is negative;
1.118 +-- 3 - A_ptr[1] is not 1;
1.119 +-- 4 - some column index is out of range;
1.120 +-- 5 - some column indices are duplicate. */
1.121 +
1.122 +int check_pattern(int m, int n, int A_ptr[], int A_ind[])
1.123 +{ int i, j, ptr, ret, *flag = NULL;
1.124 + /* check the number of rows */
1.125 + if (m < 0)
1.126 + { ret = 1;
1.127 + goto done;
1.128 + }
1.129 + /* check the number of columns */
1.130 + if (n < 0)
1.131 + { ret = 2;
1.132 + goto done;
1.133 + }
1.134 + /* check location A_ptr[1] */
1.135 + if (A_ptr[1] != 1)
1.136 + { ret = 3;
1.137 + goto done;
1.138 + }
1.139 + /* check row patterns */
1.140 + flag = xcalloc(1+n, sizeof(int));
1.141 + for (j = 1; j <= n; j++) flag[j] = 0;
1.142 + for (i = 1; i <= m; i++)
1.143 + { /* check pattern of row i */
1.144 + for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
1.145 + { j = A_ind[ptr];
1.146 + /* check column index */
1.147 + if (!(1 <= j && j <= n))
1.148 + { ret = 4;
1.149 + goto done;
1.150 + }
1.151 + /* check for duplication */
1.152 + if (flag[j])
1.153 + { ret = 5;
1.154 + goto done;
1.155 + }
1.156 + flag[j] = 1;
1.157 + }
1.158 + /* clear flags */
1.159 + for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
1.160 + { j = A_ind[ptr];
1.161 + flag[j] = 0;
1.162 + }
1.163 + }
1.164 + /* the pattern is ok */
1.165 + ret = 0;
1.166 +done: if (flag != NULL) xfree(flag);
1.167 + return ret;
1.168 +}
1.169 +
1.170 +/*----------------------------------------------------------------------
1.171 +-- transpose - transpose sparse matrix.
1.172 +--
1.173 +-- *Synopsis*
1.174 +--
1.175 +-- #include "glpmat.h"
1.176 +-- void transpose(int m, int n, int A_ptr[], int A_ind[],
1.177 +-- double A_val[], int AT_ptr[], int AT_ind[], double AT_val[]);
1.178 +--
1.179 +-- *Description*
1.180 +--
1.181 +-- For a given mxn sparse matrix A the routine transpose builds a nxm
1.182 +-- sparse matrix A' which is a matrix transposed to A.
1.183 +--
1.184 +-- The arrays A_ptr, A_ind, and A_val specify a given mxn matrix A to
1.185 +-- be transposed in storage-by-rows format. The parameter A_val can be
1.186 +-- NULL, in which case numeric values are not copied. The arrays A_ptr,
1.187 +-- A_ind, and A_val are not changed on exit.
1.188 +--
1.189 +-- On entry the arrays AT_ptr, AT_ind, and AT_val must be allocated,
1.190 +-- but their content is ignored. On exit the routine stores a resultant
1.191 +-- nxm matrix A' in these arrays in storage-by-rows format. Note that
1.192 +-- if the parameter A_val is NULL, the array AT_val is not used.
1.193 +--
1.194 +-- The routine transpose has a side effect that elements in rows of the
1.195 +-- resultant matrix A' follow in ascending their column indices. */
1.196 +
1.197 +void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[],
1.198 + int AT_ptr[], int AT_ind[], double AT_val[])
1.199 +{ int i, j, t, beg, end, pos, len;
1.200 + /* determine row lengths of resultant matrix */
1.201 + for (j = 1; j <= n; j++) AT_ptr[j] = 0;
1.202 + for (i = 1; i <= m; i++)
1.203 + { beg = A_ptr[i], end = A_ptr[i+1];
1.204 + for (t = beg; t < end; t++) AT_ptr[A_ind[t]]++;
1.205 + }
1.206 + /* set up row pointers of resultant matrix */
1.207 + pos = 1;
1.208 + for (j = 1; j <= n; j++)
1.209 + len = AT_ptr[j], pos += len, AT_ptr[j] = pos;
1.210 + AT_ptr[n+1] = pos;
1.211 + /* build resultant matrix */
1.212 + for (i = m; i >= 1; i--)
1.213 + { beg = A_ptr[i], end = A_ptr[i+1];
1.214 + for (t = beg; t < end; t++)
1.215 + { pos = --AT_ptr[A_ind[t]];
1.216 + AT_ind[pos] = i;
1.217 + if (A_val != NULL) AT_val[pos] = A_val[t];
1.218 + }
1.219 + }
1.220 + return;
1.221 +}
1.222 +
1.223 +/*----------------------------------------------------------------------
1.224 +-- adat_symbolic - compute S = P*A*D*A'*P' (symbolic phase).
1.225 +--
1.226 +-- *Synopsis*
1.227 +--
1.228 +-- #include "glpmat.h"
1.229 +-- int *adat_symbolic(int m, int n, int P_per[], int A_ptr[],
1.230 +-- int A_ind[], int S_ptr[]);
1.231 +--
1.232 +-- *Description*
1.233 +--
1.234 +-- The routine adat_symbolic implements the symbolic phase to compute
1.235 +-- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
1.236 +-- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
1.237 +-- transposed to A, P' is an inverse of P.
1.238 +--
1.239 +-- The parameter m is the number of rows in A and the order of P.
1.240 +--
1.241 +-- The parameter n is the number of columns in A and the order of D.
1.242 +--
1.243 +-- The array P_per specifies permutation matrix P. It is not changed on
1.244 +-- exit.
1.245 +--
1.246 +-- The arrays A_ptr and A_ind specify the pattern of matrix A. They are
1.247 +-- not changed on exit.
1.248 +--
1.249 +-- On exit the routine stores the pattern of upper triangular part of
1.250 +-- matrix S without diagonal elements in the arrays S_ptr and S_ind in
1.251 +-- storage-by-rows format. The array S_ptr should be allocated on entry,
1.252 +-- however, its content is ignored. The array S_ind is allocated by the
1.253 +-- routine itself which returns a pointer to it.
1.254 +--
1.255 +-- *Returns*
1.256 +--
1.257 +-- The routine returns a pointer to the array S_ind. */
1.258 +
1.259 +int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[],
1.260 + int S_ptr[])
1.261 +{ int i, j, t, ii, jj, tt, k, size, len;
1.262 + int *S_ind, *AT_ptr, *AT_ind, *ind, *map, *temp;
1.263 + /* build the pattern of A', which is a matrix transposed to A, to
1.264 + efficiently access A in column-wise manner */
1.265 + AT_ptr = xcalloc(1+n+1, sizeof(int));
1.266 + AT_ind = xcalloc(A_ptr[m+1], sizeof(int));
1.267 + transpose(m, n, A_ptr, A_ind, NULL, AT_ptr, AT_ind, NULL);
1.268 + /* allocate the array S_ind */
1.269 + size = A_ptr[m+1] - 1;
1.270 + if (size < m) size = m;
1.271 + S_ind = xcalloc(1+size, sizeof(int));
1.272 + /* allocate and initialize working arrays */
1.273 + ind = xcalloc(1+m, sizeof(int));
1.274 + map = xcalloc(1+m, sizeof(int));
1.275 + for (jj = 1; jj <= m; jj++) map[jj] = 0;
1.276 + /* compute pattern of S; note that symbolically S = B*B', where
1.277 + B = P*A, B' is matrix transposed to B */
1.278 + S_ptr[1] = 1;
1.279 + for (ii = 1; ii <= m; ii++)
1.280 + { /* compute pattern of ii-th row of S */
1.281 + len = 0;
1.282 + i = P_per[ii]; /* i-th row of A = ii-th row of B */
1.283 + for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
1.284 + { k = A_ind[t];
1.285 + /* walk through k-th column of A */
1.286 + for (tt = AT_ptr[k]; tt < AT_ptr[k+1]; tt++)
1.287 + { j = AT_ind[tt];
1.288 + jj = P_per[m+j]; /* j-th row of A = jj-th row of B */
1.289 + /* a[i,k] != 0 and a[j,k] != 0 ergo s[ii,jj] != 0 */
1.290 + if (ii < jj && !map[jj]) ind[++len] = jj, map[jj] = 1;
1.291 + }
1.292 + }
1.293 + /* now (ind) is pattern of ii-th row of S */
1.294 + S_ptr[ii+1] = S_ptr[ii] + len;
1.295 + /* at least (S_ptr[ii+1] - 1) locations should be available in
1.296 + the array S_ind */
1.297 + if (S_ptr[ii+1] - 1 > size)
1.298 + { temp = S_ind;
1.299 + size += size;
1.300 + S_ind = xcalloc(1+size, sizeof(int));
1.301 + memcpy(&S_ind[1], &temp[1], (S_ptr[ii] - 1) * sizeof(int));
1.302 + xfree(temp);
1.303 + }
1.304 + xassert(S_ptr[ii+1] - 1 <= size);
1.305 + /* (ii-th row of S) := (ind) */
1.306 + memcpy(&S_ind[S_ptr[ii]], &ind[1], len * sizeof(int));
1.307 + /* clear the row pattern map */
1.308 + for (t = 1; t <= len; t++) map[ind[t]] = 0;
1.309 + }
1.310 + /* free working arrays */
1.311 + xfree(AT_ptr);
1.312 + xfree(AT_ind);
1.313 + xfree(ind);
1.314 + xfree(map);
1.315 + /* reallocate the array S_ind to free unused locations */
1.316 + temp = S_ind;
1.317 + size = S_ptr[m+1] - 1;
1.318 + S_ind = xcalloc(1+size, sizeof(int));
1.319 + memcpy(&S_ind[1], &temp[1], size * sizeof(int));
1.320 + xfree(temp);
1.321 + return S_ind;
1.322 +}
1.323 +
1.324 +/*----------------------------------------------------------------------
1.325 +-- adat_numeric - compute S = P*A*D*A'*P' (numeric phase).
1.326 +--
1.327 +-- *Synopsis*
1.328 +--
1.329 +-- #include "glpmat.h"
1.330 +-- void adat_numeric(int m, int n, int P_per[],
1.331 +-- int A_ptr[], int A_ind[], double A_val[], double D_diag[],
1.332 +-- int S_ptr[], int S_ind[], double S_val[], double S_diag[]);
1.333 +--
1.334 +-- *Description*
1.335 +--
1.336 +-- The routine adat_numeric implements the numeric phase to compute
1.337 +-- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
1.338 +-- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
1.339 +-- transposed to A, P' is an inverse of P.
1.340 +--
1.341 +-- The parameter m is the number of rows in A and the order of P.
1.342 +--
1.343 +-- The parameter n is the number of columns in A and the order of D.
1.344 +--
1.345 +-- The matrix P is specified in the array P_per, which is not changed
1.346 +-- on exit.
1.347 +--
1.348 +-- The matrix A is specified in the arrays A_ptr, A_ind, and A_val in
1.349 +-- storage-by-rows format. These arrays are not changed on exit.
1.350 +--
1.351 +-- Diagonal elements of the matrix D are specified in the array D_diag,
1.352 +-- where D_diag[0] is not used, D_diag[i] = d[i,i] for i = 1, ..., n.
1.353 +-- The array D_diag is not changed on exit.
1.354 +--
1.355 +-- The pattern of the upper triangular part of the matrix S without
1.356 +-- diagonal elements (previously computed by the routine adat_symbolic)
1.357 +-- is specified in the arrays S_ptr and S_ind, which are not changed on
1.358 +-- exit. Numeric values of non-diagonal elements of S are stored in
1.359 +-- corresponding locations of the array S_val, and values of diagonal
1.360 +-- elements of S are stored in locations S_diag[1], ..., S_diag[n]. */
1.361 +
1.362 +void adat_numeric(int m, int n, int P_per[],
1.363 + int A_ptr[], int A_ind[], double A_val[], double D_diag[],
1.364 + int S_ptr[], int S_ind[], double S_val[], double S_diag[])
1.365 +{ int i, j, t, ii, jj, tt, beg, end, beg1, end1, k;
1.366 + double sum, *work;
1.367 + work = xcalloc(1+n, sizeof(double));
1.368 + for (j = 1; j <= n; j++) work[j] = 0.0;
1.369 + /* compute S = B*D*B', where B = P*A, B' is a matrix transposed
1.370 + to B */
1.371 + for (ii = 1; ii <= m; ii++)
1.372 + { i = P_per[ii]; /* i-th row of A = ii-th row of B */
1.373 + /* (work) := (i-th row of A) */
1.374 + beg = A_ptr[i], end = A_ptr[i+1];
1.375 + for (t = beg; t < end; t++)
1.376 + work[A_ind[t]] = A_val[t];
1.377 + /* compute ii-th row of S */
1.378 + beg = S_ptr[ii], end = S_ptr[ii+1];
1.379 + for (t = beg; t < end; t++)
1.380 + { jj = S_ind[t];
1.381 + j = P_per[jj]; /* j-th row of A = jj-th row of B */
1.382 + /* s[ii,jj] := sum a[i,k] * d[k,k] * a[j,k] */
1.383 + sum = 0.0;
1.384 + beg1 = A_ptr[j], end1 = A_ptr[j+1];
1.385 + for (tt = beg1; tt < end1; tt++)
1.386 + { k = A_ind[tt];
1.387 + sum += work[k] * D_diag[k] * A_val[tt];
1.388 + }
1.389 + S_val[t] = sum;
1.390 + }
1.391 + /* s[ii,ii] := sum a[i,k] * d[k,k] * a[i,k] */
1.392 + sum = 0.0;
1.393 + beg = A_ptr[i], end = A_ptr[i+1];
1.394 + for (t = beg; t < end; t++)
1.395 + { k = A_ind[t];
1.396 + sum += A_val[t] * D_diag[k] * A_val[t];
1.397 + work[k] = 0.0;
1.398 + }
1.399 + S_diag[ii] = sum;
1.400 + }
1.401 + xfree(work);
1.402 + return;
1.403 +}
1.404 +
1.405 +/*----------------------------------------------------------------------
1.406 +-- min_degree - minimum degree ordering.
1.407 +--
1.408 +-- *Synopsis*
1.409 +--
1.410 +-- #include "glpmat.h"
1.411 +-- void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]);
1.412 +--
1.413 +-- *Description*
1.414 +--
1.415 +-- The routine min_degree uses the minimum degree ordering algorithm
1.416 +-- to find a permutation matrix P for a given sparse symmetric positive
1.417 +-- matrix A which minimizes the number of non-zeros in upper triangular
1.418 +-- factor U for Cholesky factorization P*A*P' = U'*U.
1.419 +--
1.420 +-- The parameter n is the order of matrices A and P.
1.421 +--
1.422 +-- The pattern of the given matrix A is specified on entry in the arrays
1.423 +-- A_ptr and A_ind in storage-by-rows format. Only the upper triangular
1.424 +-- part without diagonal elements (which all are assumed to be non-zero)
1.425 +-- should be specified as if A were upper triangular. The arrays A_ptr
1.426 +-- and A_ind are not changed on exit.
1.427 +--
1.428 +-- The permutation matrix P is stored by the routine in the array P_per
1.429 +-- on exit.
1.430 +--
1.431 +-- *Algorithm*
1.432 +--
1.433 +-- The routine min_degree is based on some subroutines from the package
1.434 +-- SPARSPAK (see comments in the module glpqmd). */
1.435 +
1.436 +void min_degree(int n, int A_ptr[], int A_ind[], int P_per[])
1.437 +{ int i, j, ne, t, pos, len;
1.438 + int *xadj, *adjncy, *deg, *marker, *rchset, *nbrhd, *qsize,
1.439 + *qlink, nofsub;
1.440 + /* determine number of non-zeros in complete pattern */
1.441 + ne = A_ptr[n+1] - 1;
1.442 + ne += ne;
1.443 + /* allocate working arrays */
1.444 + xadj = xcalloc(1+n+1, sizeof(int));
1.445 + adjncy = xcalloc(1+ne, sizeof(int));
1.446 + deg = xcalloc(1+n, sizeof(int));
1.447 + marker = xcalloc(1+n, sizeof(int));
1.448 + rchset = xcalloc(1+n, sizeof(int));
1.449 + nbrhd = xcalloc(1+n, sizeof(int));
1.450 + qsize = xcalloc(1+n, sizeof(int));
1.451 + qlink = xcalloc(1+n, sizeof(int));
1.452 + /* determine row lengths in complete pattern */
1.453 + for (i = 1; i <= n; i++) xadj[i] = 0;
1.454 + for (i = 1; i <= n; i++)
1.455 + { for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
1.456 + { j = A_ind[t];
1.457 + xassert(i < j && j <= n);
1.458 + xadj[i]++, xadj[j]++;
1.459 + }
1.460 + }
1.461 + /* set up row pointers for complete pattern */
1.462 + pos = 1;
1.463 + for (i = 1; i <= n; i++)
1.464 + len = xadj[i], pos += len, xadj[i] = pos;
1.465 + xadj[n+1] = pos;
1.466 + xassert(pos - 1 == ne);
1.467 + /* construct complete pattern */
1.468 + for (i = 1; i <= n; i++)
1.469 + { for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
1.470 + { j = A_ind[t];
1.471 + adjncy[--xadj[i]] = j, adjncy[--xadj[j]] = i;
1.472 + }
1.473 + }
1.474 + /* call the main minimimum degree ordering routine */
1.475 + genqmd(&n, xadj, adjncy, P_per, P_per + n, deg, marker, rchset,
1.476 + nbrhd, qsize, qlink, &nofsub);
1.477 + /* make sure that permutation matrix P is correct */
1.478 + for (i = 1; i <= n; i++)
1.479 + { j = P_per[i];
1.480 + xassert(1 <= j && j <= n);
1.481 + xassert(P_per[n+j] == i);
1.482 + }
1.483 + /* free working arrays */
1.484 + xfree(xadj);
1.485 + xfree(adjncy);
1.486 + xfree(deg);
1.487 + xfree(marker);
1.488 + xfree(rchset);
1.489 + xfree(nbrhd);
1.490 + xfree(qsize);
1.491 + xfree(qlink);
1.492 + return;
1.493 +}
1.494 +
1.495 +/**********************************************************************/
1.496 +
1.497 +void amd_order1(int n, int A_ptr[], int A_ind[], int P_per[])
1.498 +{ /* approximate minimum degree ordering (AMD) */
1.499 + int k, ret;
1.500 + double Control[AMD_CONTROL], Info[AMD_INFO];
1.501 + /* get the default parameters */
1.502 + amd_defaults(Control);
1.503 +#if 0
1.504 + /* and print them */
1.505 + amd_control(Control);
1.506 +#endif
1.507 + /* make all indices 0-based */
1.508 + for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--;
1.509 + for (k = 1; k <= n+1; k++) A_ptr[k]--;
1.510 + /* call the ordering routine */
1.511 + ret = amd_order(n, &A_ptr[1], &A_ind[1], &P_per[1], Control, Info)
1.512 + ;
1.513 +#if 0
1.514 + amd_info(Info);
1.515 +#endif
1.516 + xassert(ret == AMD_OK || ret == AMD_OK_BUT_JUMBLED);
1.517 + /* retsore 1-based indices */
1.518 + for (k = 1; k <= n+1; k++) A_ptr[k]++;
1.519 + for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++;
1.520 + /* patch up permutation matrix */
1.521 + memset(&P_per[n+1], 0, n * sizeof(int));
1.522 + for (k = 1; k <= n; k++)
1.523 + { P_per[k]++;
1.524 + xassert(1 <= P_per[k] && P_per[k] <= n);
1.525 + xassert(P_per[n+P_per[k]] == 0);
1.526 + P_per[n+P_per[k]] = k;
1.527 + }
1.528 + return;
1.529 +}
1.530 +
1.531 +/**********************************************************************/
1.532 +
1.533 +static void *allocate(size_t n, size_t size)
1.534 +{ void *ptr;
1.535 + ptr = xcalloc(n, size);
1.536 + memset(ptr, 0, n * size);
1.537 + return ptr;
1.538 +}
1.539 +
1.540 +static void release(void *ptr)
1.541 +{ xfree(ptr);
1.542 + return;
1.543 +}
1.544 +
1.545 +void symamd_ord(int n, int A_ptr[], int A_ind[], int P_per[])
1.546 +{ /* approximate minimum degree ordering (SYMAMD) */
1.547 + int k, ok;
1.548 + int stats[COLAMD_STATS];
1.549 + /* make all indices 0-based */
1.550 + for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--;
1.551 + for (k = 1; k <= n+1; k++) A_ptr[k]--;
1.552 + /* call the ordering routine */
1.553 + ok = symamd(n, &A_ind[1], &A_ptr[1], &P_per[1], NULL, stats,
1.554 + allocate, release);
1.555 +#if 0
1.556 + symamd_report(stats);
1.557 +#endif
1.558 + xassert(ok);
1.559 + /* restore 1-based indices */
1.560 + for (k = 1; k <= n+1; k++) A_ptr[k]++;
1.561 + for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++;
1.562 + /* patch up permutation matrix */
1.563 + memset(&P_per[n+1], 0, n * sizeof(int));
1.564 + for (k = 1; k <= n; k++)
1.565 + { P_per[k]++;
1.566 + xassert(1 <= P_per[k] && P_per[k] <= n);
1.567 + xassert(P_per[n+P_per[k]] == 0);
1.568 + P_per[n+P_per[k]] = k;
1.569 + }
1.570 + return;
1.571 +}
1.572 +
1.573 +/*----------------------------------------------------------------------
1.574 +-- chol_symbolic - compute Cholesky factorization (symbolic phase).
1.575 +--
1.576 +-- *Synopsis*
1.577 +--
1.578 +-- #include "glpmat.h"
1.579 +-- int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]);
1.580 +--
1.581 +-- *Description*
1.582 +--
1.583 +-- The routine chol_symbolic implements the symbolic phase of Cholesky
1.584 +-- factorization A = U'*U, where A is a given sparse symmetric positive
1.585 +-- definite matrix, U is a resultant upper triangular factor, U' is a
1.586 +-- matrix transposed to U.
1.587 +--
1.588 +-- The parameter n is the order of matrices A and U.
1.589 +--
1.590 +-- The pattern of the given matrix A is specified on entry in the arrays
1.591 +-- A_ptr and A_ind in storage-by-rows format. Only the upper triangular
1.592 +-- part without diagonal elements (which all are assumed to be non-zero)
1.593 +-- should be specified as if A were upper triangular. The arrays A_ptr
1.594 +-- and A_ind are not changed on exit.
1.595 +--
1.596 +-- The pattern of the matrix U without diagonal elements (which all are
1.597 +-- assumed to be non-zero) is stored on exit from the routine in the
1.598 +-- arrays U_ptr and U_ind in storage-by-rows format. The array U_ptr
1.599 +-- should be allocated on entry, however, its content is ignored. The
1.600 +-- array U_ind is allocated by the routine which returns a pointer to it
1.601 +-- on exit.
1.602 +--
1.603 +-- *Returns*
1.604 +--
1.605 +-- The routine returns a pointer to the array U_ind.
1.606 +--
1.607 +-- *Method*
1.608 +--
1.609 +-- The routine chol_symbolic computes the pattern of the matrix U in a
1.610 +-- row-wise manner. No pivoting is used.
1.611 +--
1.612 +-- It is known that to compute the pattern of row k of the matrix U we
1.613 +-- need to merge the pattern of row k of the matrix A and the patterns
1.614 +-- of each row i of U, where u[i,k] is non-zero (these rows are already
1.615 +-- computed and placed above row k).
1.616 +--
1.617 +-- However, to reduce the number of rows to be merged the routine uses
1.618 +-- an advanced algorithm proposed in:
1.619 +--
1.620 +-- D.J.Rose, R.E.Tarjan, and G.S.Lueker. Algorithmic aspects of vertex
1.621 +-- elimination on graphs. SIAM J. Comput. 5, 1976, 266-83.
1.622 +--
1.623 +-- The authors of the cited paper show that we have the same result if
1.624 +-- we merge row k of the matrix A and such rows of the matrix U (among
1.625 +-- rows 1, ..., k-1) whose leftmost non-diagonal non-zero element is
1.626 +-- placed in k-th column. This feature signficantly reduces the number
1.627 +-- of rows to be merged, especially on the final steps, where rows of
1.628 +-- the matrix U become quite dense.
1.629 +--
1.630 +-- To determine rows, which should be merged on k-th step, for a fixed
1.631 +-- time the routine uses linked lists of row numbers of the matrix U.
1.632 +-- Location head[k] contains the number of a first row, whose leftmost
1.633 +-- non-diagonal non-zero element is placed in column k, and location
1.634 +-- next[i] contains the number of a next row with the same property as
1.635 +-- row i. */
1.636 +
1.637 +int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[])
1.638 +{ int i, j, k, t, len, size, beg, end, min_j, *U_ind, *head, *next,
1.639 + *ind, *map, *temp;
1.640 + /* initially we assume that on computing the pattern of U fill-in
1.641 + will double the number of non-zeros in A */
1.642 + size = A_ptr[n+1] - 1;
1.643 + if (size < n) size = n;
1.644 + size += size;
1.645 + U_ind = xcalloc(1+size, sizeof(int));
1.646 + /* allocate and initialize working arrays */
1.647 + head = xcalloc(1+n, sizeof(int));
1.648 + for (i = 1; i <= n; i++) head[i] = 0;
1.649 + next = xcalloc(1+n, sizeof(int));
1.650 + ind = xcalloc(1+n, sizeof(int));
1.651 + map = xcalloc(1+n, sizeof(int));
1.652 + for (j = 1; j <= n; j++) map[j] = 0;
1.653 + /* compute the pattern of matrix U */
1.654 + U_ptr[1] = 1;
1.655 + for (k = 1; k <= n; k++)
1.656 + { /* compute the pattern of k-th row of U, which is the union of
1.657 + k-th row of A and those rows of U (among 1, ..., k-1) whose
1.658 + leftmost non-diagonal non-zero is placed in k-th column */
1.659 + /* (ind) := (k-th row of A) */
1.660 + len = A_ptr[k+1] - A_ptr[k];
1.661 + memcpy(&ind[1], &A_ind[A_ptr[k]], len * sizeof(int));
1.662 + for (t = 1; t <= len; t++)
1.663 + { j = ind[t];
1.664 + xassert(k < j && j <= n);
1.665 + map[j] = 1;
1.666 + }
1.667 + /* walk through rows of U whose leftmost non-diagonal non-zero
1.668 + is placed in k-th column */
1.669 + for (i = head[k]; i != 0; i = next[i])
1.670 + { /* (ind) := (ind) union (i-th row of U) */
1.671 + beg = U_ptr[i], end = U_ptr[i+1];
1.672 + for (t = beg; t < end; t++)
1.673 + { j = U_ind[t];
1.674 + if (j > k && !map[j]) ind[++len] = j, map[j] = 1;
1.675 + }
1.676 + }
1.677 + /* now (ind) is the pattern of k-th row of U */
1.678 + U_ptr[k+1] = U_ptr[k] + len;
1.679 + /* at least (U_ptr[k+1] - 1) locations should be available in
1.680 + the array U_ind */
1.681 + if (U_ptr[k+1] - 1 > size)
1.682 + { temp = U_ind;
1.683 + size += size;
1.684 + U_ind = xcalloc(1+size, sizeof(int));
1.685 + memcpy(&U_ind[1], &temp[1], (U_ptr[k] - 1) * sizeof(int));
1.686 + xfree(temp);
1.687 + }
1.688 + xassert(U_ptr[k+1] - 1 <= size);
1.689 + /* (k-th row of U) := (ind) */
1.690 + memcpy(&U_ind[U_ptr[k]], &ind[1], len * sizeof(int));
1.691 + /* determine column index of leftmost non-diagonal non-zero in
1.692 + k-th row of U and clear the row pattern map */
1.693 + min_j = n + 1;
1.694 + for (t = 1; t <= len; t++)
1.695 + { j = ind[t], map[j] = 0;
1.696 + if (min_j > j) min_j = j;
1.697 + }
1.698 + /* include k-th row into corresponding linked list */
1.699 + if (min_j <= n) next[k] = head[min_j], head[min_j] = k;
1.700 + }
1.701 + /* free working arrays */
1.702 + xfree(head);
1.703 + xfree(next);
1.704 + xfree(ind);
1.705 + xfree(map);
1.706 + /* reallocate the array U_ind to free unused locations */
1.707 + temp = U_ind;
1.708 + size = U_ptr[n+1] - 1;
1.709 + U_ind = xcalloc(1+size, sizeof(int));
1.710 + memcpy(&U_ind[1], &temp[1], size * sizeof(int));
1.711 + xfree(temp);
1.712 + return U_ind;
1.713 +}
1.714 +
1.715 +/*----------------------------------------------------------------------
1.716 +-- chol_numeric - compute Cholesky factorization (numeric phase).
1.717 +--
1.718 +-- *Synopsis*
1.719 +--
1.720 +-- #include "glpmat.h"
1.721 +-- int chol_numeric(int n,
1.722 +-- int A_ptr[], int A_ind[], double A_val[], double A_diag[],
1.723 +-- int U_ptr[], int U_ind[], double U_val[], double U_diag[]);
1.724 +--
1.725 +-- *Description*
1.726 +--
1.727 +-- The routine chol_symbolic implements the numeric phase of Cholesky
1.728 +-- factorization A = U'*U, where A is a given sparse symmetric positive
1.729 +-- definite matrix, U is a resultant upper triangular factor, U' is a
1.730 +-- matrix transposed to U.
1.731 +--
1.732 +-- The parameter n is the order of matrices A and U.
1.733 +--
1.734 +-- Upper triangular part of the matrix A without diagonal elements is
1.735 +-- specified in the arrays A_ptr, A_ind, and A_val in storage-by-rows
1.736 +-- format. Diagonal elements of A are specified in the array A_diag,
1.737 +-- where A_diag[0] is not used, A_diag[i] = a[i,i] for i = 1, ..., n.
1.738 +-- The arrays A_ptr, A_ind, A_val, and A_diag are not changed on exit.
1.739 +--
1.740 +-- The pattern of the matrix U without diagonal elements (previously
1.741 +-- computed with the routine chol_symbolic) is specified in the arrays
1.742 +-- U_ptr and U_ind, which are not changed on exit. Numeric values of
1.743 +-- non-diagonal elements of U are stored in corresponding locations of
1.744 +-- the array U_val, and values of diagonal elements of U are stored in
1.745 +-- locations U_diag[1], ..., U_diag[n].
1.746 +--
1.747 +-- *Returns*
1.748 +--
1.749 +-- The routine returns the number of non-positive diagonal elements of
1.750 +-- the matrix U which have been replaced by a huge positive number (see
1.751 +-- the method description below). Zero return code means the matrix A
1.752 +-- has been successfully factorized.
1.753 +--
1.754 +-- *Method*
1.755 +--
1.756 +-- The routine chol_numeric computes the matrix U in a row-wise manner
1.757 +-- using standard gaussian elimination technique. No pivoting is used.
1.758 +--
1.759 +-- Initially the routine sets U = A, and before k-th elimination step
1.760 +-- the matrix U is the following:
1.761 +--
1.762 +-- 1 k n
1.763 +-- 1 x x x x x x x x x x
1.764 +-- . x x x x x x x x x
1.765 +-- . . x x x x x x x x
1.766 +-- . . . x x x x x x x
1.767 +-- k . . . . * * * * * *
1.768 +-- . . . . * * * * * *
1.769 +-- . . . . * * * * * *
1.770 +-- . . . . * * * * * *
1.771 +-- . . . . * * * * * *
1.772 +-- n . . . . * * * * * *
1.773 +--
1.774 +-- where 'x' are elements of already computed rows, '*' are elements of
1.775 +-- the active submatrix. (Note that the lower triangular part of the
1.776 +-- active submatrix being symmetric is not stored and diagonal elements
1.777 +-- are stored separately in the array U_diag.)
1.778 +--
1.779 +-- The matrix A is assumed to be positive definite. However, if it is
1.780 +-- close to semi-definite, on some elimination step a pivot u[k,k] may
1.781 +-- happen to be non-positive due to round-off errors. In this case the
1.782 +-- routine uses a technique proposed in:
1.783 +--
1.784 +-- S.J.Wright. The Cholesky factorization in interior-point and barrier
1.785 +-- methods. Preprint MCS-P600-0596, Mathematics and Computer Science
1.786 +-- Division, Argonne National Laboratory, Argonne, Ill., May 1996.
1.787 +--
1.788 +-- The routine just replaces non-positive u[k,k] by a huge positive
1.789 +-- number. This involves non-diagonal elements in k-th row of U to be
1.790 +-- close to zero that, in turn, involves k-th component of a solution
1.791 +-- vector to be close to zero. Note, however, that this technique works
1.792 +-- only if the system A*x = b is consistent. */
1.793 +
1.794 +int chol_numeric(int n,
1.795 + int A_ptr[], int A_ind[], double A_val[], double A_diag[],
1.796 + int U_ptr[], int U_ind[], double U_val[], double U_diag[])
1.797 +{ int i, j, k, t, t1, beg, end, beg1, end1, count = 0;
1.798 + double ukk, uki, *work;
1.799 + work = xcalloc(1+n, sizeof(double));
1.800 + for (j = 1; j <= n; j++) work[j] = 0.0;
1.801 + /* U := (upper triangle of A) */
1.802 + /* note that the upper traingle of A is a subset of U */
1.803 + for (i = 1; i <= n; i++)
1.804 + { beg = A_ptr[i], end = A_ptr[i+1];
1.805 + for (t = beg; t < end; t++)
1.806 + j = A_ind[t], work[j] = A_val[t];
1.807 + beg = U_ptr[i], end = U_ptr[i+1];
1.808 + for (t = beg; t < end; t++)
1.809 + j = U_ind[t], U_val[t] = work[j], work[j] = 0.0;
1.810 + U_diag[i] = A_diag[i];
1.811 + }
1.812 + /* main elimination loop */
1.813 + for (k = 1; k <= n; k++)
1.814 + { /* transform k-th row of U */
1.815 + ukk = U_diag[k];
1.816 + if (ukk > 0.0)
1.817 + U_diag[k] = ukk = sqrt(ukk);
1.818 + else
1.819 + U_diag[k] = ukk = DBL_MAX, count++;
1.820 + /* (work) := (transformed k-th row) */
1.821 + beg = U_ptr[k], end = U_ptr[k+1];
1.822 + for (t = beg; t < end; t++)
1.823 + work[U_ind[t]] = (U_val[t] /= ukk);
1.824 + /* transform other rows of U */
1.825 + for (t = beg; t < end; t++)
1.826 + { i = U_ind[t];
1.827 + xassert(i > k);
1.828 + /* (i-th row) := (i-th row) - u[k,i] * (k-th row) */
1.829 + uki = work[i];
1.830 + beg1 = U_ptr[i], end1 = U_ptr[i+1];
1.831 + for (t1 = beg1; t1 < end1; t1++)
1.832 + U_val[t1] -= uki * work[U_ind[t1]];
1.833 + U_diag[i] -= uki * uki;
1.834 + }
1.835 + /* (work) := 0 */
1.836 + for (t = beg; t < end; t++)
1.837 + work[U_ind[t]] = 0.0;
1.838 + }
1.839 + xfree(work);
1.840 + return count;
1.841 +}
1.842 +
1.843 +/*----------------------------------------------------------------------
1.844 +-- u_solve - solve upper triangular system U*x = b.
1.845 +--
1.846 +-- *Synopsis*
1.847 +--
1.848 +-- #include "glpmat.h"
1.849 +-- void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
1.850 +-- double U_diag[], double x[]);
1.851 +--
1.852 +-- *Description*
1.853 +--
1.854 +-- The routine u_solve solves an linear system U*x = b, where U is an
1.855 +-- upper triangular matrix.
1.856 +--
1.857 +-- The parameter n is the order of matrix U.
1.858 +--
1.859 +-- The matrix U without diagonal elements is specified in the arrays
1.860 +-- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
1.861 +-- of U are specified in the array U_diag, where U_diag[0] is not used,
1.862 +-- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
1.863 +-- changed on exit.
1.864 +--
1.865 +-- The right-hand side vector b is specified on entry in the array x,
1.866 +-- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
1.867 +-- the routine stores computed components of the vector of unknowns x
1.868 +-- in the array x in the same manner. */
1.869 +
1.870 +void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
1.871 + double U_diag[], double x[])
1.872 +{ int i, t, beg, end;
1.873 + double temp;
1.874 + for (i = n; i >= 1; i--)
1.875 + { temp = x[i];
1.876 + beg = U_ptr[i], end = U_ptr[i+1];
1.877 + for (t = beg; t < end; t++)
1.878 + temp -= U_val[t] * x[U_ind[t]];
1.879 + xassert(U_diag[i] != 0.0);
1.880 + x[i] = temp / U_diag[i];
1.881 + }
1.882 + return;
1.883 +}
1.884 +
1.885 +/*----------------------------------------------------------------------
1.886 +-- ut_solve - solve lower triangular system U'*x = b.
1.887 +--
1.888 +-- *Synopsis*
1.889 +--
1.890 +-- #include "glpmat.h"
1.891 +-- void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
1.892 +-- double U_diag[], double x[]);
1.893 +--
1.894 +-- *Description*
1.895 +--
1.896 +-- The routine ut_solve solves an linear system U'*x = b, where U is a
1.897 +-- matrix transposed to an upper triangular matrix.
1.898 +--
1.899 +-- The parameter n is the order of matrix U.
1.900 +--
1.901 +-- The matrix U without diagonal elements is specified in the arrays
1.902 +-- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
1.903 +-- of U are specified in the array U_diag, where U_diag[0] is not used,
1.904 +-- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
1.905 +-- changed on exit.
1.906 +--
1.907 +-- The right-hand side vector b is specified on entry in the array x,
1.908 +-- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
1.909 +-- the routine stores computed components of the vector of unknowns x
1.910 +-- in the array x in the same manner. */
1.911 +
1.912 +void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
1.913 + double U_diag[], double x[])
1.914 +{ int i, t, beg, end;
1.915 + double temp;
1.916 + for (i = 1; i <= n; i++)
1.917 + { xassert(U_diag[i] != 0.0);
1.918 + temp = (x[i] /= U_diag[i]);
1.919 + if (temp == 0.0) continue;
1.920 + beg = U_ptr[i], end = U_ptr[i+1];
1.921 + for (t = beg; t < end; t++)
1.922 + x[U_ind[t]] -= U_val[t] * temp;
1.923 + }
1.924 + return;
1.925 +}
1.926 +
1.927 +/* eof */