lemon-project-template-glpk

annotate deps/glpk/src/glpmat.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
rev   line source
alpar@9 1 /* glpmat.c */
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 #include "glpenv.h"
alpar@9 26 #include "glpmat.h"
alpar@9 27 #include "glpqmd.h"
alpar@9 28 #include "amd/amd.h"
alpar@9 29 #include "colamd/colamd.h"
alpar@9 30
alpar@9 31 /*----------------------------------------------------------------------
alpar@9 32 -- check_fvs - check sparse vector in full-vector storage format.
alpar@9 33 --
alpar@9 34 -- SYNOPSIS
alpar@9 35 --
alpar@9 36 -- #include "glpmat.h"
alpar@9 37 -- int check_fvs(int n, int nnz, int ind[], double vec[]);
alpar@9 38 --
alpar@9 39 -- DESCRIPTION
alpar@9 40 --
alpar@9 41 -- The routine check_fvs checks if a given vector of dimension n in
alpar@9 42 -- full-vector storage format has correct representation.
alpar@9 43 --
alpar@9 44 -- RETURNS
alpar@9 45 --
alpar@9 46 -- The routine returns one of the following codes:
alpar@9 47 --
alpar@9 48 -- 0 - the vector is correct;
alpar@9 49 -- 1 - the number of elements (n) is negative;
alpar@9 50 -- 2 - the number of non-zero elements (nnz) is negative;
alpar@9 51 -- 3 - some element index is out of range;
alpar@9 52 -- 4 - some element index is duplicate;
alpar@9 53 -- 5 - some non-zero element is out of pattern. */
alpar@9 54
alpar@9 55 int check_fvs(int n, int nnz, int ind[], double vec[])
alpar@9 56 { int i, t, ret, *flag = NULL;
alpar@9 57 /* check the number of elements */
alpar@9 58 if (n < 0)
alpar@9 59 { ret = 1;
alpar@9 60 goto done;
alpar@9 61 }
alpar@9 62 /* check the number of non-zero elements */
alpar@9 63 if (nnz < 0)
alpar@9 64 { ret = 2;
alpar@9 65 goto done;
alpar@9 66 }
alpar@9 67 /* check vector indices */
alpar@9 68 flag = xcalloc(1+n, sizeof(int));
alpar@9 69 for (i = 1; i <= n; i++) flag[i] = 0;
alpar@9 70 for (t = 1; t <= nnz; t++)
alpar@9 71 { i = ind[t];
alpar@9 72 if (!(1 <= i && i <= n))
alpar@9 73 { ret = 3;
alpar@9 74 goto done;
alpar@9 75 }
alpar@9 76 if (flag[i])
alpar@9 77 { ret = 4;
alpar@9 78 goto done;
alpar@9 79 }
alpar@9 80 flag[i] = 1;
alpar@9 81 }
alpar@9 82 /* check vector elements */
alpar@9 83 for (i = 1; i <= n; i++)
alpar@9 84 { if (!flag[i] && vec[i] != 0.0)
alpar@9 85 { ret = 5;
alpar@9 86 goto done;
alpar@9 87 }
alpar@9 88 }
alpar@9 89 /* the vector is ok */
alpar@9 90 ret = 0;
alpar@9 91 done: if (flag != NULL) xfree(flag);
alpar@9 92 return ret;
alpar@9 93 }
alpar@9 94
alpar@9 95 /*----------------------------------------------------------------------
alpar@9 96 -- check_pattern - check pattern of sparse matrix.
alpar@9 97 --
alpar@9 98 -- SYNOPSIS
alpar@9 99 --
alpar@9 100 -- #include "glpmat.h"
alpar@9 101 -- int check_pattern(int m, int n, int A_ptr[], int A_ind[]);
alpar@9 102 --
alpar@9 103 -- DESCRIPTION
alpar@9 104 --
alpar@9 105 -- The routine check_pattern checks the pattern of a given mxn matrix
alpar@9 106 -- in storage-by-rows format.
alpar@9 107 --
alpar@9 108 -- RETURNS
alpar@9 109 --
alpar@9 110 -- The routine returns one of the following codes:
alpar@9 111 --
alpar@9 112 -- 0 - the pattern is correct;
alpar@9 113 -- 1 - the number of rows (m) is negative;
alpar@9 114 -- 2 - the number of columns (n) is negative;
alpar@9 115 -- 3 - A_ptr[1] is not 1;
alpar@9 116 -- 4 - some column index is out of range;
alpar@9 117 -- 5 - some column indices are duplicate. */
alpar@9 118
alpar@9 119 int check_pattern(int m, int n, int A_ptr[], int A_ind[])
alpar@9 120 { int i, j, ptr, ret, *flag = NULL;
alpar@9 121 /* check the number of rows */
alpar@9 122 if (m < 0)
alpar@9 123 { ret = 1;
alpar@9 124 goto done;
alpar@9 125 }
alpar@9 126 /* check the number of columns */
alpar@9 127 if (n < 0)
alpar@9 128 { ret = 2;
alpar@9 129 goto done;
alpar@9 130 }
alpar@9 131 /* check location A_ptr[1] */
alpar@9 132 if (A_ptr[1] != 1)
alpar@9 133 { ret = 3;
alpar@9 134 goto done;
alpar@9 135 }
alpar@9 136 /* check row patterns */
alpar@9 137 flag = xcalloc(1+n, sizeof(int));
alpar@9 138 for (j = 1; j <= n; j++) flag[j] = 0;
alpar@9 139 for (i = 1; i <= m; i++)
alpar@9 140 { /* check pattern of row i */
alpar@9 141 for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
alpar@9 142 { j = A_ind[ptr];
alpar@9 143 /* check column index */
alpar@9 144 if (!(1 <= j && j <= n))
alpar@9 145 { ret = 4;
alpar@9 146 goto done;
alpar@9 147 }
alpar@9 148 /* check for duplication */
alpar@9 149 if (flag[j])
alpar@9 150 { ret = 5;
alpar@9 151 goto done;
alpar@9 152 }
alpar@9 153 flag[j] = 1;
alpar@9 154 }
alpar@9 155 /* clear flags */
alpar@9 156 for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
alpar@9 157 { j = A_ind[ptr];
alpar@9 158 flag[j] = 0;
alpar@9 159 }
alpar@9 160 }
alpar@9 161 /* the pattern is ok */
alpar@9 162 ret = 0;
alpar@9 163 done: if (flag != NULL) xfree(flag);
alpar@9 164 return ret;
alpar@9 165 }
alpar@9 166
alpar@9 167 /*----------------------------------------------------------------------
alpar@9 168 -- transpose - transpose sparse matrix.
alpar@9 169 --
alpar@9 170 -- *Synopsis*
alpar@9 171 --
alpar@9 172 -- #include "glpmat.h"
alpar@9 173 -- void transpose(int m, int n, int A_ptr[], int A_ind[],
alpar@9 174 -- double A_val[], int AT_ptr[], int AT_ind[], double AT_val[]);
alpar@9 175 --
alpar@9 176 -- *Description*
alpar@9 177 --
alpar@9 178 -- For a given mxn sparse matrix A the routine transpose builds a nxm
alpar@9 179 -- sparse matrix A' which is a matrix transposed to A.
alpar@9 180 --
alpar@9 181 -- The arrays A_ptr, A_ind, and A_val specify a given mxn matrix A to
alpar@9 182 -- be transposed in storage-by-rows format. The parameter A_val can be
alpar@9 183 -- NULL, in which case numeric values are not copied. The arrays A_ptr,
alpar@9 184 -- A_ind, and A_val are not changed on exit.
alpar@9 185 --
alpar@9 186 -- On entry the arrays AT_ptr, AT_ind, and AT_val must be allocated,
alpar@9 187 -- but their content is ignored. On exit the routine stores a resultant
alpar@9 188 -- nxm matrix A' in these arrays in storage-by-rows format. Note that
alpar@9 189 -- if the parameter A_val is NULL, the array AT_val is not used.
alpar@9 190 --
alpar@9 191 -- The routine transpose has a side effect that elements in rows of the
alpar@9 192 -- resultant matrix A' follow in ascending their column indices. */
alpar@9 193
alpar@9 194 void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[],
alpar@9 195 int AT_ptr[], int AT_ind[], double AT_val[])
alpar@9 196 { int i, j, t, beg, end, pos, len;
alpar@9 197 /* determine row lengths of resultant matrix */
alpar@9 198 for (j = 1; j <= n; j++) AT_ptr[j] = 0;
alpar@9 199 for (i = 1; i <= m; i++)
alpar@9 200 { beg = A_ptr[i], end = A_ptr[i+1];
alpar@9 201 for (t = beg; t < end; t++) AT_ptr[A_ind[t]]++;
alpar@9 202 }
alpar@9 203 /* set up row pointers of resultant matrix */
alpar@9 204 pos = 1;
alpar@9 205 for (j = 1; j <= n; j++)
alpar@9 206 len = AT_ptr[j], pos += len, AT_ptr[j] = pos;
alpar@9 207 AT_ptr[n+1] = pos;
alpar@9 208 /* build resultant matrix */
alpar@9 209 for (i = m; i >= 1; i--)
alpar@9 210 { beg = A_ptr[i], end = A_ptr[i+1];
alpar@9 211 for (t = beg; t < end; t++)
alpar@9 212 { pos = --AT_ptr[A_ind[t]];
alpar@9 213 AT_ind[pos] = i;
alpar@9 214 if (A_val != NULL) AT_val[pos] = A_val[t];
alpar@9 215 }
alpar@9 216 }
alpar@9 217 return;
alpar@9 218 }
alpar@9 219
alpar@9 220 /*----------------------------------------------------------------------
alpar@9 221 -- adat_symbolic - compute S = P*A*D*A'*P' (symbolic phase).
alpar@9 222 --
alpar@9 223 -- *Synopsis*
alpar@9 224 --
alpar@9 225 -- #include "glpmat.h"
alpar@9 226 -- int *adat_symbolic(int m, int n, int P_per[], int A_ptr[],
alpar@9 227 -- int A_ind[], int S_ptr[]);
alpar@9 228 --
alpar@9 229 -- *Description*
alpar@9 230 --
alpar@9 231 -- The routine adat_symbolic implements the symbolic phase to compute
alpar@9 232 -- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
alpar@9 233 -- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
alpar@9 234 -- transposed to A, P' is an inverse of P.
alpar@9 235 --
alpar@9 236 -- The parameter m is the number of rows in A and the order of P.
alpar@9 237 --
alpar@9 238 -- The parameter n is the number of columns in A and the order of D.
alpar@9 239 --
alpar@9 240 -- The array P_per specifies permutation matrix P. It is not changed on
alpar@9 241 -- exit.
alpar@9 242 --
alpar@9 243 -- The arrays A_ptr and A_ind specify the pattern of matrix A. They are
alpar@9 244 -- not changed on exit.
alpar@9 245 --
alpar@9 246 -- On exit the routine stores the pattern of upper triangular part of
alpar@9 247 -- matrix S without diagonal elements in the arrays S_ptr and S_ind in
alpar@9 248 -- storage-by-rows format. The array S_ptr should be allocated on entry,
alpar@9 249 -- however, its content is ignored. The array S_ind is allocated by the
alpar@9 250 -- routine itself which returns a pointer to it.
alpar@9 251 --
alpar@9 252 -- *Returns*
alpar@9 253 --
alpar@9 254 -- The routine returns a pointer to the array S_ind. */
alpar@9 255
alpar@9 256 int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[],
alpar@9 257 int S_ptr[])
alpar@9 258 { int i, j, t, ii, jj, tt, k, size, len;
alpar@9 259 int *S_ind, *AT_ptr, *AT_ind, *ind, *map, *temp;
alpar@9 260 /* build the pattern of A', which is a matrix transposed to A, to
alpar@9 261 efficiently access A in column-wise manner */
alpar@9 262 AT_ptr = xcalloc(1+n+1, sizeof(int));
alpar@9 263 AT_ind = xcalloc(A_ptr[m+1], sizeof(int));
alpar@9 264 transpose(m, n, A_ptr, A_ind, NULL, AT_ptr, AT_ind, NULL);
alpar@9 265 /* allocate the array S_ind */
alpar@9 266 size = A_ptr[m+1] - 1;
alpar@9 267 if (size < m) size = m;
alpar@9 268 S_ind = xcalloc(1+size, sizeof(int));
alpar@9 269 /* allocate and initialize working arrays */
alpar@9 270 ind = xcalloc(1+m, sizeof(int));
alpar@9 271 map = xcalloc(1+m, sizeof(int));
alpar@9 272 for (jj = 1; jj <= m; jj++) map[jj] = 0;
alpar@9 273 /* compute pattern of S; note that symbolically S = B*B', where
alpar@9 274 B = P*A, B' is matrix transposed to B */
alpar@9 275 S_ptr[1] = 1;
alpar@9 276 for (ii = 1; ii <= m; ii++)
alpar@9 277 { /* compute pattern of ii-th row of S */
alpar@9 278 len = 0;
alpar@9 279 i = P_per[ii]; /* i-th row of A = ii-th row of B */
alpar@9 280 for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
alpar@9 281 { k = A_ind[t];
alpar@9 282 /* walk through k-th column of A */
alpar@9 283 for (tt = AT_ptr[k]; tt < AT_ptr[k+1]; tt++)
alpar@9 284 { j = AT_ind[tt];
alpar@9 285 jj = P_per[m+j]; /* j-th row of A = jj-th row of B */
alpar@9 286 /* a[i,k] != 0 and a[j,k] != 0 ergo s[ii,jj] != 0 */
alpar@9 287 if (ii < jj && !map[jj]) ind[++len] = jj, map[jj] = 1;
alpar@9 288 }
alpar@9 289 }
alpar@9 290 /* now (ind) is pattern of ii-th row of S */
alpar@9 291 S_ptr[ii+1] = S_ptr[ii] + len;
alpar@9 292 /* at least (S_ptr[ii+1] - 1) locations should be available in
alpar@9 293 the array S_ind */
alpar@9 294 if (S_ptr[ii+1] - 1 > size)
alpar@9 295 { temp = S_ind;
alpar@9 296 size += size;
alpar@9 297 S_ind = xcalloc(1+size, sizeof(int));
alpar@9 298 memcpy(&S_ind[1], &temp[1], (S_ptr[ii] - 1) * sizeof(int));
alpar@9 299 xfree(temp);
alpar@9 300 }
alpar@9 301 xassert(S_ptr[ii+1] - 1 <= size);
alpar@9 302 /* (ii-th row of S) := (ind) */
alpar@9 303 memcpy(&S_ind[S_ptr[ii]], &ind[1], len * sizeof(int));
alpar@9 304 /* clear the row pattern map */
alpar@9 305 for (t = 1; t <= len; t++) map[ind[t]] = 0;
alpar@9 306 }
alpar@9 307 /* free working arrays */
alpar@9 308 xfree(AT_ptr);
alpar@9 309 xfree(AT_ind);
alpar@9 310 xfree(ind);
alpar@9 311 xfree(map);
alpar@9 312 /* reallocate the array S_ind to free unused locations */
alpar@9 313 temp = S_ind;
alpar@9 314 size = S_ptr[m+1] - 1;
alpar@9 315 S_ind = xcalloc(1+size, sizeof(int));
alpar@9 316 memcpy(&S_ind[1], &temp[1], size * sizeof(int));
alpar@9 317 xfree(temp);
alpar@9 318 return S_ind;
alpar@9 319 }
alpar@9 320
alpar@9 321 /*----------------------------------------------------------------------
alpar@9 322 -- adat_numeric - compute S = P*A*D*A'*P' (numeric phase).
alpar@9 323 --
alpar@9 324 -- *Synopsis*
alpar@9 325 --
alpar@9 326 -- #include "glpmat.h"
alpar@9 327 -- void adat_numeric(int m, int n, int P_per[],
alpar@9 328 -- int A_ptr[], int A_ind[], double A_val[], double D_diag[],
alpar@9 329 -- int S_ptr[], int S_ind[], double S_val[], double S_diag[]);
alpar@9 330 --
alpar@9 331 -- *Description*
alpar@9 332 --
alpar@9 333 -- The routine adat_numeric implements the numeric phase to compute
alpar@9 334 -- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
alpar@9 335 -- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
alpar@9 336 -- transposed to A, P' is an inverse of P.
alpar@9 337 --
alpar@9 338 -- The parameter m is the number of rows in A and the order of P.
alpar@9 339 --
alpar@9 340 -- The parameter n is the number of columns in A and the order of D.
alpar@9 341 --
alpar@9 342 -- The matrix P is specified in the array P_per, which is not changed
alpar@9 343 -- on exit.
alpar@9 344 --
alpar@9 345 -- The matrix A is specified in the arrays A_ptr, A_ind, and A_val in
alpar@9 346 -- storage-by-rows format. These arrays are not changed on exit.
alpar@9 347 --
alpar@9 348 -- Diagonal elements of the matrix D are specified in the array D_diag,
alpar@9 349 -- where D_diag[0] is not used, D_diag[i] = d[i,i] for i = 1, ..., n.
alpar@9 350 -- The array D_diag is not changed on exit.
alpar@9 351 --
alpar@9 352 -- The pattern of the upper triangular part of the matrix S without
alpar@9 353 -- diagonal elements (previously computed by the routine adat_symbolic)
alpar@9 354 -- is specified in the arrays S_ptr and S_ind, which are not changed on
alpar@9 355 -- exit. Numeric values of non-diagonal elements of S are stored in
alpar@9 356 -- corresponding locations of the array S_val, and values of diagonal
alpar@9 357 -- elements of S are stored in locations S_diag[1], ..., S_diag[n]. */
alpar@9 358
alpar@9 359 void adat_numeric(int m, int n, int P_per[],
alpar@9 360 int A_ptr[], int A_ind[], double A_val[], double D_diag[],
alpar@9 361 int S_ptr[], int S_ind[], double S_val[], double S_diag[])
alpar@9 362 { int i, j, t, ii, jj, tt, beg, end, beg1, end1, k;
alpar@9 363 double sum, *work;
alpar@9 364 work = xcalloc(1+n, sizeof(double));
alpar@9 365 for (j = 1; j <= n; j++) work[j] = 0.0;
alpar@9 366 /* compute S = B*D*B', where B = P*A, B' is a matrix transposed
alpar@9 367 to B */
alpar@9 368 for (ii = 1; ii <= m; ii++)
alpar@9 369 { i = P_per[ii]; /* i-th row of A = ii-th row of B */
alpar@9 370 /* (work) := (i-th row of A) */
alpar@9 371 beg = A_ptr[i], end = A_ptr[i+1];
alpar@9 372 for (t = beg; t < end; t++)
alpar@9 373 work[A_ind[t]] = A_val[t];
alpar@9 374 /* compute ii-th row of S */
alpar@9 375 beg = S_ptr[ii], end = S_ptr[ii+1];
alpar@9 376 for (t = beg; t < end; t++)
alpar@9 377 { jj = S_ind[t];
alpar@9 378 j = P_per[jj]; /* j-th row of A = jj-th row of B */
alpar@9 379 /* s[ii,jj] := sum a[i,k] * d[k,k] * a[j,k] */
alpar@9 380 sum = 0.0;
alpar@9 381 beg1 = A_ptr[j], end1 = A_ptr[j+1];
alpar@9 382 for (tt = beg1; tt < end1; tt++)
alpar@9 383 { k = A_ind[tt];
alpar@9 384 sum += work[k] * D_diag[k] * A_val[tt];
alpar@9 385 }
alpar@9 386 S_val[t] = sum;
alpar@9 387 }
alpar@9 388 /* s[ii,ii] := sum a[i,k] * d[k,k] * a[i,k] */
alpar@9 389 sum = 0.0;
alpar@9 390 beg = A_ptr[i], end = A_ptr[i+1];
alpar@9 391 for (t = beg; t < end; t++)
alpar@9 392 { k = A_ind[t];
alpar@9 393 sum += A_val[t] * D_diag[k] * A_val[t];
alpar@9 394 work[k] = 0.0;
alpar@9 395 }
alpar@9 396 S_diag[ii] = sum;
alpar@9 397 }
alpar@9 398 xfree(work);
alpar@9 399 return;
alpar@9 400 }
alpar@9 401
alpar@9 402 /*----------------------------------------------------------------------
alpar@9 403 -- min_degree - minimum degree ordering.
alpar@9 404 --
alpar@9 405 -- *Synopsis*
alpar@9 406 --
alpar@9 407 -- #include "glpmat.h"
alpar@9 408 -- void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]);
alpar@9 409 --
alpar@9 410 -- *Description*
alpar@9 411 --
alpar@9 412 -- The routine min_degree uses the minimum degree ordering algorithm
alpar@9 413 -- to find a permutation matrix P for a given sparse symmetric positive
alpar@9 414 -- matrix A which minimizes the number of non-zeros in upper triangular
alpar@9 415 -- factor U for Cholesky factorization P*A*P' = U'*U.
alpar@9 416 --
alpar@9 417 -- The parameter n is the order of matrices A and P.
alpar@9 418 --
alpar@9 419 -- The pattern of the given matrix A is specified on entry in the arrays
alpar@9 420 -- A_ptr and A_ind in storage-by-rows format. Only the upper triangular
alpar@9 421 -- part without diagonal elements (which all are assumed to be non-zero)
alpar@9 422 -- should be specified as if A were upper triangular. The arrays A_ptr
alpar@9 423 -- and A_ind are not changed on exit.
alpar@9 424 --
alpar@9 425 -- The permutation matrix P is stored by the routine in the array P_per
alpar@9 426 -- on exit.
alpar@9 427 --
alpar@9 428 -- *Algorithm*
alpar@9 429 --
alpar@9 430 -- The routine min_degree is based on some subroutines from the package
alpar@9 431 -- SPARSPAK (see comments in the module glpqmd). */
alpar@9 432
alpar@9 433 void min_degree(int n, int A_ptr[], int A_ind[], int P_per[])
alpar@9 434 { int i, j, ne, t, pos, len;
alpar@9 435 int *xadj, *adjncy, *deg, *marker, *rchset, *nbrhd, *qsize,
alpar@9 436 *qlink, nofsub;
alpar@9 437 /* determine number of non-zeros in complete pattern */
alpar@9 438 ne = A_ptr[n+1] - 1;
alpar@9 439 ne += ne;
alpar@9 440 /* allocate working arrays */
alpar@9 441 xadj = xcalloc(1+n+1, sizeof(int));
alpar@9 442 adjncy = xcalloc(1+ne, sizeof(int));
alpar@9 443 deg = xcalloc(1+n, sizeof(int));
alpar@9 444 marker = xcalloc(1+n, sizeof(int));
alpar@9 445 rchset = xcalloc(1+n, sizeof(int));
alpar@9 446 nbrhd = xcalloc(1+n, sizeof(int));
alpar@9 447 qsize = xcalloc(1+n, sizeof(int));
alpar@9 448 qlink = xcalloc(1+n, sizeof(int));
alpar@9 449 /* determine row lengths in complete pattern */
alpar@9 450 for (i = 1; i <= n; i++) xadj[i] = 0;
alpar@9 451 for (i = 1; i <= n; i++)
alpar@9 452 { for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
alpar@9 453 { j = A_ind[t];
alpar@9 454 xassert(i < j && j <= n);
alpar@9 455 xadj[i]++, xadj[j]++;
alpar@9 456 }
alpar@9 457 }
alpar@9 458 /* set up row pointers for complete pattern */
alpar@9 459 pos = 1;
alpar@9 460 for (i = 1; i <= n; i++)
alpar@9 461 len = xadj[i], pos += len, xadj[i] = pos;
alpar@9 462 xadj[n+1] = pos;
alpar@9 463 xassert(pos - 1 == ne);
alpar@9 464 /* construct complete pattern */
alpar@9 465 for (i = 1; i <= n; i++)
alpar@9 466 { for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
alpar@9 467 { j = A_ind[t];
alpar@9 468 adjncy[--xadj[i]] = j, adjncy[--xadj[j]] = i;
alpar@9 469 }
alpar@9 470 }
alpar@9 471 /* call the main minimimum degree ordering routine */
alpar@9 472 genqmd(&n, xadj, adjncy, P_per, P_per + n, deg, marker, rchset,
alpar@9 473 nbrhd, qsize, qlink, &nofsub);
alpar@9 474 /* make sure that permutation matrix P is correct */
alpar@9 475 for (i = 1; i <= n; i++)
alpar@9 476 { j = P_per[i];
alpar@9 477 xassert(1 <= j && j <= n);
alpar@9 478 xassert(P_per[n+j] == i);
alpar@9 479 }
alpar@9 480 /* free working arrays */
alpar@9 481 xfree(xadj);
alpar@9 482 xfree(adjncy);
alpar@9 483 xfree(deg);
alpar@9 484 xfree(marker);
alpar@9 485 xfree(rchset);
alpar@9 486 xfree(nbrhd);
alpar@9 487 xfree(qsize);
alpar@9 488 xfree(qlink);
alpar@9 489 return;
alpar@9 490 }
alpar@9 491
alpar@9 492 /**********************************************************************/
alpar@9 493
alpar@9 494 void amd_order1(int n, int A_ptr[], int A_ind[], int P_per[])
alpar@9 495 { /* approximate minimum degree ordering (AMD) */
alpar@9 496 int k, ret;
alpar@9 497 double Control[AMD_CONTROL], Info[AMD_INFO];
alpar@9 498 /* get the default parameters */
alpar@9 499 amd_defaults(Control);
alpar@9 500 #if 0
alpar@9 501 /* and print them */
alpar@9 502 amd_control(Control);
alpar@9 503 #endif
alpar@9 504 /* make all indices 0-based */
alpar@9 505 for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--;
alpar@9 506 for (k = 1; k <= n+1; k++) A_ptr[k]--;
alpar@9 507 /* call the ordering routine */
alpar@9 508 ret = amd_order(n, &A_ptr[1], &A_ind[1], &P_per[1], Control, Info)
alpar@9 509 ;
alpar@9 510 #if 0
alpar@9 511 amd_info(Info);
alpar@9 512 #endif
alpar@9 513 xassert(ret == AMD_OK || ret == AMD_OK_BUT_JUMBLED);
alpar@9 514 /* retsore 1-based indices */
alpar@9 515 for (k = 1; k <= n+1; k++) A_ptr[k]++;
alpar@9 516 for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++;
alpar@9 517 /* patch up permutation matrix */
alpar@9 518 memset(&P_per[n+1], 0, n * sizeof(int));
alpar@9 519 for (k = 1; k <= n; k++)
alpar@9 520 { P_per[k]++;
alpar@9 521 xassert(1 <= P_per[k] && P_per[k] <= n);
alpar@9 522 xassert(P_per[n+P_per[k]] == 0);
alpar@9 523 P_per[n+P_per[k]] = k;
alpar@9 524 }
alpar@9 525 return;
alpar@9 526 }
alpar@9 527
alpar@9 528 /**********************************************************************/
alpar@9 529
alpar@9 530 static void *allocate(size_t n, size_t size)
alpar@9 531 { void *ptr;
alpar@9 532 ptr = xcalloc(n, size);
alpar@9 533 memset(ptr, 0, n * size);
alpar@9 534 return ptr;
alpar@9 535 }
alpar@9 536
alpar@9 537 static void release(void *ptr)
alpar@9 538 { xfree(ptr);
alpar@9 539 return;
alpar@9 540 }
alpar@9 541
alpar@9 542 void symamd_ord(int n, int A_ptr[], int A_ind[], int P_per[])
alpar@9 543 { /* approximate minimum degree ordering (SYMAMD) */
alpar@9 544 int k, ok;
alpar@9 545 int stats[COLAMD_STATS];
alpar@9 546 /* make all indices 0-based */
alpar@9 547 for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--;
alpar@9 548 for (k = 1; k <= n+1; k++) A_ptr[k]--;
alpar@9 549 /* call the ordering routine */
alpar@9 550 ok = symamd(n, &A_ind[1], &A_ptr[1], &P_per[1], NULL, stats,
alpar@9 551 allocate, release);
alpar@9 552 #if 0
alpar@9 553 symamd_report(stats);
alpar@9 554 #endif
alpar@9 555 xassert(ok);
alpar@9 556 /* restore 1-based indices */
alpar@9 557 for (k = 1; k <= n+1; k++) A_ptr[k]++;
alpar@9 558 for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++;
alpar@9 559 /* patch up permutation matrix */
alpar@9 560 memset(&P_per[n+1], 0, n * sizeof(int));
alpar@9 561 for (k = 1; k <= n; k++)
alpar@9 562 { P_per[k]++;
alpar@9 563 xassert(1 <= P_per[k] && P_per[k] <= n);
alpar@9 564 xassert(P_per[n+P_per[k]] == 0);
alpar@9 565 P_per[n+P_per[k]] = k;
alpar@9 566 }
alpar@9 567 return;
alpar@9 568 }
alpar@9 569
alpar@9 570 /*----------------------------------------------------------------------
alpar@9 571 -- chol_symbolic - compute Cholesky factorization (symbolic phase).
alpar@9 572 --
alpar@9 573 -- *Synopsis*
alpar@9 574 --
alpar@9 575 -- #include "glpmat.h"
alpar@9 576 -- int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]);
alpar@9 577 --
alpar@9 578 -- *Description*
alpar@9 579 --
alpar@9 580 -- The routine chol_symbolic implements the symbolic phase of Cholesky
alpar@9 581 -- factorization A = U'*U, where A is a given sparse symmetric positive
alpar@9 582 -- definite matrix, U is a resultant upper triangular factor, U' is a
alpar@9 583 -- matrix transposed to U.
alpar@9 584 --
alpar@9 585 -- The parameter n is the order of matrices A and U.
alpar@9 586 --
alpar@9 587 -- The pattern of the given matrix A is specified on entry in the arrays
alpar@9 588 -- A_ptr and A_ind in storage-by-rows format. Only the upper triangular
alpar@9 589 -- part without diagonal elements (which all are assumed to be non-zero)
alpar@9 590 -- should be specified as if A were upper triangular. The arrays A_ptr
alpar@9 591 -- and A_ind are not changed on exit.
alpar@9 592 --
alpar@9 593 -- The pattern of the matrix U without diagonal elements (which all are
alpar@9 594 -- assumed to be non-zero) is stored on exit from the routine in the
alpar@9 595 -- arrays U_ptr and U_ind in storage-by-rows format. The array U_ptr
alpar@9 596 -- should be allocated on entry, however, its content is ignored. The
alpar@9 597 -- array U_ind is allocated by the routine which returns a pointer to it
alpar@9 598 -- on exit.
alpar@9 599 --
alpar@9 600 -- *Returns*
alpar@9 601 --
alpar@9 602 -- The routine returns a pointer to the array U_ind.
alpar@9 603 --
alpar@9 604 -- *Method*
alpar@9 605 --
alpar@9 606 -- The routine chol_symbolic computes the pattern of the matrix U in a
alpar@9 607 -- row-wise manner. No pivoting is used.
alpar@9 608 --
alpar@9 609 -- It is known that to compute the pattern of row k of the matrix U we
alpar@9 610 -- need to merge the pattern of row k of the matrix A and the patterns
alpar@9 611 -- of each row i of U, where u[i,k] is non-zero (these rows are already
alpar@9 612 -- computed and placed above row k).
alpar@9 613 --
alpar@9 614 -- However, to reduce the number of rows to be merged the routine uses
alpar@9 615 -- an advanced algorithm proposed in:
alpar@9 616 --
alpar@9 617 -- D.J.Rose, R.E.Tarjan, and G.S.Lueker. Algorithmic aspects of vertex
alpar@9 618 -- elimination on graphs. SIAM J. Comput. 5, 1976, 266-83.
alpar@9 619 --
alpar@9 620 -- The authors of the cited paper show that we have the same result if
alpar@9 621 -- we merge row k of the matrix A and such rows of the matrix U (among
alpar@9 622 -- rows 1, ..., k-1) whose leftmost non-diagonal non-zero element is
alpar@9 623 -- placed in k-th column. This feature signficantly reduces the number
alpar@9 624 -- of rows to be merged, especially on the final steps, where rows of
alpar@9 625 -- the matrix U become quite dense.
alpar@9 626 --
alpar@9 627 -- To determine rows, which should be merged on k-th step, for a fixed
alpar@9 628 -- time the routine uses linked lists of row numbers of the matrix U.
alpar@9 629 -- Location head[k] contains the number of a first row, whose leftmost
alpar@9 630 -- non-diagonal non-zero element is placed in column k, and location
alpar@9 631 -- next[i] contains the number of a next row with the same property as
alpar@9 632 -- row i. */
alpar@9 633
alpar@9 634 int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[])
alpar@9 635 { int i, j, k, t, len, size, beg, end, min_j, *U_ind, *head, *next,
alpar@9 636 *ind, *map, *temp;
alpar@9 637 /* initially we assume that on computing the pattern of U fill-in
alpar@9 638 will double the number of non-zeros in A */
alpar@9 639 size = A_ptr[n+1] - 1;
alpar@9 640 if (size < n) size = n;
alpar@9 641 size += size;
alpar@9 642 U_ind = xcalloc(1+size, sizeof(int));
alpar@9 643 /* allocate and initialize working arrays */
alpar@9 644 head = xcalloc(1+n, sizeof(int));
alpar@9 645 for (i = 1; i <= n; i++) head[i] = 0;
alpar@9 646 next = xcalloc(1+n, sizeof(int));
alpar@9 647 ind = xcalloc(1+n, sizeof(int));
alpar@9 648 map = xcalloc(1+n, sizeof(int));
alpar@9 649 for (j = 1; j <= n; j++) map[j] = 0;
alpar@9 650 /* compute the pattern of matrix U */
alpar@9 651 U_ptr[1] = 1;
alpar@9 652 for (k = 1; k <= n; k++)
alpar@9 653 { /* compute the pattern of k-th row of U, which is the union of
alpar@9 654 k-th row of A and those rows of U (among 1, ..., k-1) whose
alpar@9 655 leftmost non-diagonal non-zero is placed in k-th column */
alpar@9 656 /* (ind) := (k-th row of A) */
alpar@9 657 len = A_ptr[k+1] - A_ptr[k];
alpar@9 658 memcpy(&ind[1], &A_ind[A_ptr[k]], len * sizeof(int));
alpar@9 659 for (t = 1; t <= len; t++)
alpar@9 660 { j = ind[t];
alpar@9 661 xassert(k < j && j <= n);
alpar@9 662 map[j] = 1;
alpar@9 663 }
alpar@9 664 /* walk through rows of U whose leftmost non-diagonal non-zero
alpar@9 665 is placed in k-th column */
alpar@9 666 for (i = head[k]; i != 0; i = next[i])
alpar@9 667 { /* (ind) := (ind) union (i-th row of U) */
alpar@9 668 beg = U_ptr[i], end = U_ptr[i+1];
alpar@9 669 for (t = beg; t < end; t++)
alpar@9 670 { j = U_ind[t];
alpar@9 671 if (j > k && !map[j]) ind[++len] = j, map[j] = 1;
alpar@9 672 }
alpar@9 673 }
alpar@9 674 /* now (ind) is the pattern of k-th row of U */
alpar@9 675 U_ptr[k+1] = U_ptr[k] + len;
alpar@9 676 /* at least (U_ptr[k+1] - 1) locations should be available in
alpar@9 677 the array U_ind */
alpar@9 678 if (U_ptr[k+1] - 1 > size)
alpar@9 679 { temp = U_ind;
alpar@9 680 size += size;
alpar@9 681 U_ind = xcalloc(1+size, sizeof(int));
alpar@9 682 memcpy(&U_ind[1], &temp[1], (U_ptr[k] - 1) * sizeof(int));
alpar@9 683 xfree(temp);
alpar@9 684 }
alpar@9 685 xassert(U_ptr[k+1] - 1 <= size);
alpar@9 686 /* (k-th row of U) := (ind) */
alpar@9 687 memcpy(&U_ind[U_ptr[k]], &ind[1], len * sizeof(int));
alpar@9 688 /* determine column index of leftmost non-diagonal non-zero in
alpar@9 689 k-th row of U and clear the row pattern map */
alpar@9 690 min_j = n + 1;
alpar@9 691 for (t = 1; t <= len; t++)
alpar@9 692 { j = ind[t], map[j] = 0;
alpar@9 693 if (min_j > j) min_j = j;
alpar@9 694 }
alpar@9 695 /* include k-th row into corresponding linked list */
alpar@9 696 if (min_j <= n) next[k] = head[min_j], head[min_j] = k;
alpar@9 697 }
alpar@9 698 /* free working arrays */
alpar@9 699 xfree(head);
alpar@9 700 xfree(next);
alpar@9 701 xfree(ind);
alpar@9 702 xfree(map);
alpar@9 703 /* reallocate the array U_ind to free unused locations */
alpar@9 704 temp = U_ind;
alpar@9 705 size = U_ptr[n+1] - 1;
alpar@9 706 U_ind = xcalloc(1+size, sizeof(int));
alpar@9 707 memcpy(&U_ind[1], &temp[1], size * sizeof(int));
alpar@9 708 xfree(temp);
alpar@9 709 return U_ind;
alpar@9 710 }
alpar@9 711
alpar@9 712 /*----------------------------------------------------------------------
alpar@9 713 -- chol_numeric - compute Cholesky factorization (numeric phase).
alpar@9 714 --
alpar@9 715 -- *Synopsis*
alpar@9 716 --
alpar@9 717 -- #include "glpmat.h"
alpar@9 718 -- int chol_numeric(int n,
alpar@9 719 -- int A_ptr[], int A_ind[], double A_val[], double A_diag[],
alpar@9 720 -- int U_ptr[], int U_ind[], double U_val[], double U_diag[]);
alpar@9 721 --
alpar@9 722 -- *Description*
alpar@9 723 --
alpar@9 724 -- The routine chol_symbolic implements the numeric phase of Cholesky
alpar@9 725 -- factorization A = U'*U, where A is a given sparse symmetric positive
alpar@9 726 -- definite matrix, U is a resultant upper triangular factor, U' is a
alpar@9 727 -- matrix transposed to U.
alpar@9 728 --
alpar@9 729 -- The parameter n is the order of matrices A and U.
alpar@9 730 --
alpar@9 731 -- Upper triangular part of the matrix A without diagonal elements is
alpar@9 732 -- specified in the arrays A_ptr, A_ind, and A_val in storage-by-rows
alpar@9 733 -- format. Diagonal elements of A are specified in the array A_diag,
alpar@9 734 -- where A_diag[0] is not used, A_diag[i] = a[i,i] for i = 1, ..., n.
alpar@9 735 -- The arrays A_ptr, A_ind, A_val, and A_diag are not changed on exit.
alpar@9 736 --
alpar@9 737 -- The pattern of the matrix U without diagonal elements (previously
alpar@9 738 -- computed with the routine chol_symbolic) is specified in the arrays
alpar@9 739 -- U_ptr and U_ind, which are not changed on exit. Numeric values of
alpar@9 740 -- non-diagonal elements of U are stored in corresponding locations of
alpar@9 741 -- the array U_val, and values of diagonal elements of U are stored in
alpar@9 742 -- locations U_diag[1], ..., U_diag[n].
alpar@9 743 --
alpar@9 744 -- *Returns*
alpar@9 745 --
alpar@9 746 -- The routine returns the number of non-positive diagonal elements of
alpar@9 747 -- the matrix U which have been replaced by a huge positive number (see
alpar@9 748 -- the method description below). Zero return code means the matrix A
alpar@9 749 -- has been successfully factorized.
alpar@9 750 --
alpar@9 751 -- *Method*
alpar@9 752 --
alpar@9 753 -- The routine chol_numeric computes the matrix U in a row-wise manner
alpar@9 754 -- using standard gaussian elimination technique. No pivoting is used.
alpar@9 755 --
alpar@9 756 -- Initially the routine sets U = A, and before k-th elimination step
alpar@9 757 -- the matrix U is the following:
alpar@9 758 --
alpar@9 759 -- 1 k n
alpar@9 760 -- 1 x x x x x x x x x x
alpar@9 761 -- . x x x x x x x x x
alpar@9 762 -- . . x x x x x x x x
alpar@9 763 -- . . . x x x x x x x
alpar@9 764 -- k . . . . * * * * * *
alpar@9 765 -- . . . . * * * * * *
alpar@9 766 -- . . . . * * * * * *
alpar@9 767 -- . . . . * * * * * *
alpar@9 768 -- . . . . * * * * * *
alpar@9 769 -- n . . . . * * * * * *
alpar@9 770 --
alpar@9 771 -- where 'x' are elements of already computed rows, '*' are elements of
alpar@9 772 -- the active submatrix. (Note that the lower triangular part of the
alpar@9 773 -- active submatrix being symmetric is not stored and diagonal elements
alpar@9 774 -- are stored separately in the array U_diag.)
alpar@9 775 --
alpar@9 776 -- The matrix A is assumed to be positive definite. However, if it is
alpar@9 777 -- close to semi-definite, on some elimination step a pivot u[k,k] may
alpar@9 778 -- happen to be non-positive due to round-off errors. In this case the
alpar@9 779 -- routine uses a technique proposed in:
alpar@9 780 --
alpar@9 781 -- S.J.Wright. The Cholesky factorization in interior-point and barrier
alpar@9 782 -- methods. Preprint MCS-P600-0596, Mathematics and Computer Science
alpar@9 783 -- Division, Argonne National Laboratory, Argonne, Ill., May 1996.
alpar@9 784 --
alpar@9 785 -- The routine just replaces non-positive u[k,k] by a huge positive
alpar@9 786 -- number. This involves non-diagonal elements in k-th row of U to be
alpar@9 787 -- close to zero that, in turn, involves k-th component of a solution
alpar@9 788 -- vector to be close to zero. Note, however, that this technique works
alpar@9 789 -- only if the system A*x = b is consistent. */
alpar@9 790
alpar@9 791 int chol_numeric(int n,
alpar@9 792 int A_ptr[], int A_ind[], double A_val[], double A_diag[],
alpar@9 793 int U_ptr[], int U_ind[], double U_val[], double U_diag[])
alpar@9 794 { int i, j, k, t, t1, beg, end, beg1, end1, count = 0;
alpar@9 795 double ukk, uki, *work;
alpar@9 796 work = xcalloc(1+n, sizeof(double));
alpar@9 797 for (j = 1; j <= n; j++) work[j] = 0.0;
alpar@9 798 /* U := (upper triangle of A) */
alpar@9 799 /* note that the upper traingle of A is a subset of U */
alpar@9 800 for (i = 1; i <= n; i++)
alpar@9 801 { beg = A_ptr[i], end = A_ptr[i+1];
alpar@9 802 for (t = beg; t < end; t++)
alpar@9 803 j = A_ind[t], work[j] = A_val[t];
alpar@9 804 beg = U_ptr[i], end = U_ptr[i+1];
alpar@9 805 for (t = beg; t < end; t++)
alpar@9 806 j = U_ind[t], U_val[t] = work[j], work[j] = 0.0;
alpar@9 807 U_diag[i] = A_diag[i];
alpar@9 808 }
alpar@9 809 /* main elimination loop */
alpar@9 810 for (k = 1; k <= n; k++)
alpar@9 811 { /* transform k-th row of U */
alpar@9 812 ukk = U_diag[k];
alpar@9 813 if (ukk > 0.0)
alpar@9 814 U_diag[k] = ukk = sqrt(ukk);
alpar@9 815 else
alpar@9 816 U_diag[k] = ukk = DBL_MAX, count++;
alpar@9 817 /* (work) := (transformed k-th row) */
alpar@9 818 beg = U_ptr[k], end = U_ptr[k+1];
alpar@9 819 for (t = beg; t < end; t++)
alpar@9 820 work[U_ind[t]] = (U_val[t] /= ukk);
alpar@9 821 /* transform other rows of U */
alpar@9 822 for (t = beg; t < end; t++)
alpar@9 823 { i = U_ind[t];
alpar@9 824 xassert(i > k);
alpar@9 825 /* (i-th row) := (i-th row) - u[k,i] * (k-th row) */
alpar@9 826 uki = work[i];
alpar@9 827 beg1 = U_ptr[i], end1 = U_ptr[i+1];
alpar@9 828 for (t1 = beg1; t1 < end1; t1++)
alpar@9 829 U_val[t1] -= uki * work[U_ind[t1]];
alpar@9 830 U_diag[i] -= uki * uki;
alpar@9 831 }
alpar@9 832 /* (work) := 0 */
alpar@9 833 for (t = beg; t < end; t++)
alpar@9 834 work[U_ind[t]] = 0.0;
alpar@9 835 }
alpar@9 836 xfree(work);
alpar@9 837 return count;
alpar@9 838 }
alpar@9 839
alpar@9 840 /*----------------------------------------------------------------------
alpar@9 841 -- u_solve - solve upper triangular system U*x = b.
alpar@9 842 --
alpar@9 843 -- *Synopsis*
alpar@9 844 --
alpar@9 845 -- #include "glpmat.h"
alpar@9 846 -- void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
alpar@9 847 -- double U_diag[], double x[]);
alpar@9 848 --
alpar@9 849 -- *Description*
alpar@9 850 --
alpar@9 851 -- The routine u_solve solves an linear system U*x = b, where U is an
alpar@9 852 -- upper triangular matrix.
alpar@9 853 --
alpar@9 854 -- The parameter n is the order of matrix U.
alpar@9 855 --
alpar@9 856 -- The matrix U without diagonal elements is specified in the arrays
alpar@9 857 -- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
alpar@9 858 -- of U are specified in the array U_diag, where U_diag[0] is not used,
alpar@9 859 -- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
alpar@9 860 -- changed on exit.
alpar@9 861 --
alpar@9 862 -- The right-hand side vector b is specified on entry in the array x,
alpar@9 863 -- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
alpar@9 864 -- the routine stores computed components of the vector of unknowns x
alpar@9 865 -- in the array x in the same manner. */
alpar@9 866
alpar@9 867 void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
alpar@9 868 double U_diag[], double x[])
alpar@9 869 { int i, t, beg, end;
alpar@9 870 double temp;
alpar@9 871 for (i = n; i >= 1; i--)
alpar@9 872 { temp = x[i];
alpar@9 873 beg = U_ptr[i], end = U_ptr[i+1];
alpar@9 874 for (t = beg; t < end; t++)
alpar@9 875 temp -= U_val[t] * x[U_ind[t]];
alpar@9 876 xassert(U_diag[i] != 0.0);
alpar@9 877 x[i] = temp / U_diag[i];
alpar@9 878 }
alpar@9 879 return;
alpar@9 880 }
alpar@9 881
alpar@9 882 /*----------------------------------------------------------------------
alpar@9 883 -- ut_solve - solve lower triangular system U'*x = b.
alpar@9 884 --
alpar@9 885 -- *Synopsis*
alpar@9 886 --
alpar@9 887 -- #include "glpmat.h"
alpar@9 888 -- void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
alpar@9 889 -- double U_diag[], double x[]);
alpar@9 890 --
alpar@9 891 -- *Description*
alpar@9 892 --
alpar@9 893 -- The routine ut_solve solves an linear system U'*x = b, where U is a
alpar@9 894 -- matrix transposed to an upper triangular matrix.
alpar@9 895 --
alpar@9 896 -- The parameter n is the order of matrix U.
alpar@9 897 --
alpar@9 898 -- The matrix U without diagonal elements is specified in the arrays
alpar@9 899 -- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
alpar@9 900 -- of U are specified in the array U_diag, where U_diag[0] is not used,
alpar@9 901 -- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
alpar@9 902 -- changed on exit.
alpar@9 903 --
alpar@9 904 -- The right-hand side vector b is specified on entry in the array x,
alpar@9 905 -- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
alpar@9 906 -- the routine stores computed components of the vector of unknowns x
alpar@9 907 -- in the array x in the same manner. */
alpar@9 908
alpar@9 909 void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
alpar@9 910 double U_diag[], double x[])
alpar@9 911 { int i, t, beg, end;
alpar@9 912 double temp;
alpar@9 913 for (i = 1; i <= n; i++)
alpar@9 914 { xassert(U_diag[i] != 0.0);
alpar@9 915 temp = (x[i] /= U_diag[i]);
alpar@9 916 if (temp == 0.0) continue;
alpar@9 917 beg = U_ptr[i], end = U_ptr[i+1];
alpar@9 918 for (t = beg; t < end; t++)
alpar@9 919 x[U_ind[t]] -= U_val[t] * temp;
alpar@9 920 }
alpar@9 921 return;
alpar@9 922 }
alpar@9 923
alpar@9 924 /* eof */