1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpspm.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,846 @@
1.4 +/* glpspm.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 "glphbm.h"
1.29 +#include "glprgr.h"
1.30 +#include "glpspm.h"
1.31 +
1.32 +/***********************************************************************
1.33 +* NAME
1.34 +*
1.35 +* spm_create_mat - create general sparse matrix
1.36 +*
1.37 +* SYNOPSIS
1.38 +*
1.39 +* #include "glpspm.h"
1.40 +* SPM *spm_create_mat(int m, int n);
1.41 +*
1.42 +* DESCRIPTION
1.43 +*
1.44 +* The routine spm_create_mat creates a general sparse matrix having
1.45 +* m rows and n columns. Being created the matrix is zero (empty), i.e.
1.46 +* has no elements.
1.47 +*
1.48 +* RETURNS
1.49 +*
1.50 +* The routine returns a pointer to the matrix created. */
1.51 +
1.52 +SPM *spm_create_mat(int m, int n)
1.53 +{ SPM *A;
1.54 + xassert(0 <= m && m < INT_MAX);
1.55 + xassert(0 <= n && n < INT_MAX);
1.56 + A = xmalloc(sizeof(SPM));
1.57 + A->m = m;
1.58 + A->n = n;
1.59 + if (m == 0 || n == 0)
1.60 + { A->pool = NULL;
1.61 + A->row = NULL;
1.62 + A->col = NULL;
1.63 + }
1.64 + else
1.65 + { int i, j;
1.66 + A->pool = dmp_create_pool();
1.67 + A->row = xcalloc(1+m, sizeof(SPME *));
1.68 + for (i = 1; i <= m; i++) A->row[i] = NULL;
1.69 + A->col = xcalloc(1+n, sizeof(SPME *));
1.70 + for (j = 1; j <= n; j++) A->col[j] = NULL;
1.71 + }
1.72 + return A;
1.73 +}
1.74 +
1.75 +/***********************************************************************
1.76 +* NAME
1.77 +*
1.78 +* spm_new_elem - add new element to sparse matrix
1.79 +*
1.80 +* SYNOPSIS
1.81 +*
1.82 +* #include "glpspm.h"
1.83 +* SPME *spm_new_elem(SPM *A, int i, int j, double val);
1.84 +*
1.85 +* DESCRIPTION
1.86 +*
1.87 +* The routine spm_new_elem adds a new element to the specified sparse
1.88 +* matrix. Parameters i, j, and val specify the row number, the column
1.89 +* number, and a numerical value of the element, respectively.
1.90 +*
1.91 +* RETURNS
1.92 +*
1.93 +* The routine returns a pointer to the new element added. */
1.94 +
1.95 +SPME *spm_new_elem(SPM *A, int i, int j, double val)
1.96 +{ SPME *e;
1.97 + xassert(1 <= i && i <= A->m);
1.98 + xassert(1 <= j && j <= A->n);
1.99 + e = dmp_get_atom(A->pool, sizeof(SPME));
1.100 + e->i = i;
1.101 + e->j = j;
1.102 + e->val = val;
1.103 + e->r_prev = NULL;
1.104 + e->r_next = A->row[i];
1.105 + if (e->r_next != NULL) e->r_next->r_prev = e;
1.106 + e->c_prev = NULL;
1.107 + e->c_next = A->col[j];
1.108 + if (e->c_next != NULL) e->c_next->c_prev = e;
1.109 + A->row[i] = A->col[j] = e;
1.110 + return e;
1.111 +}
1.112 +
1.113 +/***********************************************************************
1.114 +* NAME
1.115 +*
1.116 +* spm_delete_mat - delete general sparse matrix
1.117 +*
1.118 +* SYNOPSIS
1.119 +*
1.120 +* #include "glpspm.h"
1.121 +* void spm_delete_mat(SPM *A);
1.122 +*
1.123 +* DESCRIPTION
1.124 +*
1.125 +* The routine deletes the specified general sparse matrix freeing all
1.126 +* the memory allocated to this object. */
1.127 +
1.128 +void spm_delete_mat(SPM *A)
1.129 +{ /* delete sparse matrix */
1.130 + if (A->pool != NULL) dmp_delete_pool(A->pool);
1.131 + if (A->row != NULL) xfree(A->row);
1.132 + if (A->col != NULL) xfree(A->col);
1.133 + xfree(A);
1.134 + return;
1.135 +}
1.136 +
1.137 +/***********************************************************************
1.138 +* NAME
1.139 +*
1.140 +* spm_test_mat_e - create test sparse matrix of E(n,c) class
1.141 +*
1.142 +* SYNOPSIS
1.143 +*
1.144 +* #include "glpspm.h"
1.145 +* SPM *spm_test_mat_e(int n, int c);
1.146 +*
1.147 +* DESCRIPTION
1.148 +*
1.149 +* The routine spm_test_mat_e creates a test sparse matrix of E(n,c)
1.150 +* class as described in the book: Ole 0sterby, Zahari Zlatev. Direct
1.151 +* Methods for Sparse Matrices. Springer-Verlag, 1983.
1.152 +*
1.153 +* Matrix of E(n,c) class is a symmetric positive definite matrix of
1.154 +* the order n. It has the number 4 on its main diagonal and the number
1.155 +* -1 on its four co-diagonals, two of which are neighbour to the main
1.156 +* diagonal and two others are shifted from the main diagonal on the
1.157 +* distance c.
1.158 +*
1.159 +* It is necessary that n >= 3 and 2 <= c <= n-1.
1.160 +*
1.161 +* RETURNS
1.162 +*
1.163 +* The routine returns a pointer to the matrix created. */
1.164 +
1.165 +SPM *spm_test_mat_e(int n, int c)
1.166 +{ SPM *A;
1.167 + int i;
1.168 + xassert(n >= 3 && 2 <= c && c <= n-1);
1.169 + A = spm_create_mat(n, n);
1.170 + for (i = 1; i <= n; i++)
1.171 + spm_new_elem(A, i, i, 4.0);
1.172 + for (i = 1; i <= n-1; i++)
1.173 + { spm_new_elem(A, i, i+1, -1.0);
1.174 + spm_new_elem(A, i+1, i, -1.0);
1.175 + }
1.176 + for (i = 1; i <= n-c; i++)
1.177 + { spm_new_elem(A, i, i+c, -1.0);
1.178 + spm_new_elem(A, i+c, i, -1.0);
1.179 + }
1.180 + return A;
1.181 +}
1.182 +
1.183 +/***********************************************************************
1.184 +* NAME
1.185 +*
1.186 +* spm_test_mat_d - create test sparse matrix of D(n,c) class
1.187 +*
1.188 +* SYNOPSIS
1.189 +*
1.190 +* #include "glpspm.h"
1.191 +* SPM *spm_test_mat_d(int n, int c);
1.192 +*
1.193 +* DESCRIPTION
1.194 +*
1.195 +* The routine spm_test_mat_d creates a test sparse matrix of D(n,c)
1.196 +* class as described in the book: Ole 0sterby, Zahari Zlatev. Direct
1.197 +* Methods for Sparse Matrices. Springer-Verlag, 1983.
1.198 +*
1.199 +* Matrix of D(n,c) class is a non-singular matrix of the order n. It
1.200 +* has unity main diagonal, three co-diagonals above the main diagonal
1.201 +* on the distance c, which are cyclically continued below the main
1.202 +* diagonal, and a triangle block of the size 10x10 in the upper right
1.203 +* corner.
1.204 +*
1.205 +* It is necessary that n >= 14 and 1 <= c <= n-13.
1.206 +*
1.207 +* RETURNS
1.208 +*
1.209 +* The routine returns a pointer to the matrix created. */
1.210 +
1.211 +SPM *spm_test_mat_d(int n, int c)
1.212 +{ SPM *A;
1.213 + int i, j;
1.214 + xassert(n >= 14 && 1 <= c && c <= n-13);
1.215 + A = spm_create_mat(n, n);
1.216 + for (i = 1; i <= n; i++)
1.217 + spm_new_elem(A, i, i, 1.0);
1.218 + for (i = 1; i <= n-c; i++)
1.219 + spm_new_elem(A, i, i+c, (double)(i+1));
1.220 + for (i = n-c+1; i <= n; i++)
1.221 + spm_new_elem(A, i, i-n+c, (double)(i+1));
1.222 + for (i = 1; i <= n-c-1; i++)
1.223 + spm_new_elem(A, i, i+c+1, (double)(-i));
1.224 + for (i = n-c; i <= n; i++)
1.225 + spm_new_elem(A, i, i-n+c+1, (double)(-i));
1.226 + for (i = 1; i <= n-c-2; i++)
1.227 + spm_new_elem(A, i, i+c+2, 16.0);
1.228 + for (i = n-c-1; i <= n; i++)
1.229 + spm_new_elem(A, i, i-n+c+2, 16.0);
1.230 + for (j = 1; j <= 10; j++)
1.231 + for (i = 1; i <= 11-j; i++)
1.232 + spm_new_elem(A, i, n-11+i+j, 100.0 * (double)j);
1.233 + return A;
1.234 +}
1.235 +
1.236 +/***********************************************************************
1.237 +* NAME
1.238 +*
1.239 +* spm_show_mat - write sparse matrix pattern in BMP file format
1.240 +*
1.241 +* SYNOPSIS
1.242 +*
1.243 +* #include "glpspm.h"
1.244 +* int spm_show_mat(const SPM *A, const char *fname);
1.245 +*
1.246 +* DESCRIPTION
1.247 +*
1.248 +* The routine spm_show_mat writes pattern of the specified sparse
1.249 +* matrix in uncompressed BMP file format (Windows bitmap) to a binary
1.250 +* file whose name is specified by the character string fname.
1.251 +*
1.252 +* Each pixel corresponds to one matrix element. The pixel colors have
1.253 +* the following meaning:
1.254 +*
1.255 +* Black structurally zero element
1.256 +* White positive element
1.257 +* Cyan negative element
1.258 +* Green zero element
1.259 +* Red duplicate element
1.260 +*
1.261 +* RETURNS
1.262 +*
1.263 +* If no error occured, the routine returns zero. Otherwise, it prints
1.264 +* an appropriate error message and returns non-zero. */
1.265 +
1.266 +int spm_show_mat(const SPM *A, const char *fname)
1.267 +{ int m = A->m;
1.268 + int n = A->n;
1.269 + int i, j, k, ret;
1.270 + char *map;
1.271 + xprintf("spm_show_mat: writing matrix pattern to `%s'...\n",
1.272 + fname);
1.273 + xassert(1 <= m && m <= 32767);
1.274 + xassert(1 <= n && n <= 32767);
1.275 + map = xmalloc(m * n);
1.276 + memset(map, 0x08, m * n);
1.277 + for (i = 1; i <= m; i++)
1.278 + { SPME *e;
1.279 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.280 + { j = e->j;
1.281 + xassert(1 <= j && j <= n);
1.282 + k = n * (i - 1) + (j - 1);
1.283 + if (map[k] != 0x08)
1.284 + map[k] = 0x0C;
1.285 + else if (e->val > 0.0)
1.286 + map[k] = 0x0F;
1.287 + else if (e->val < 0.0)
1.288 + map[k] = 0x0B;
1.289 + else
1.290 + map[k] = 0x0A;
1.291 + }
1.292 + }
1.293 + ret = rgr_write_bmp16(fname, m, n, map);
1.294 + xfree(map);
1.295 + return ret;
1.296 +}
1.297 +
1.298 +/***********************************************************************
1.299 +* NAME
1.300 +*
1.301 +* spm_read_hbm - read sparse matrix in Harwell-Boeing format
1.302 +*
1.303 +* SYNOPSIS
1.304 +*
1.305 +* #include "glpspm.h"
1.306 +* SPM *spm_read_hbm(const char *fname);
1.307 +*
1.308 +* DESCRIPTION
1.309 +*
1.310 +* The routine spm_read_hbm reads a sparse matrix in the Harwell-Boeing
1.311 +* format from a text file whose name is the character string fname.
1.312 +*
1.313 +* Detailed description of the Harwell-Boeing format recognised by this
1.314 +* routine can be found in the following report:
1.315 +*
1.316 +* I.S.Duff, R.G.Grimes, J.G.Lewis. User's Guide for the Harwell-Boeing
1.317 +* Sparse Matrix Collection (Release I), TR/PA/92/86, October 1992.
1.318 +*
1.319 +* NOTE
1.320 +*
1.321 +* The routine spm_read_hbm reads the matrix "as is", due to which zero
1.322 +* and/or duplicate elements can appear in the matrix.
1.323 +*
1.324 +* RETURNS
1.325 +*
1.326 +* If no error occured, the routine returns a pointer to the matrix
1.327 +* created. Otherwise, the routine prints an appropriate error message
1.328 +* and returns NULL. */
1.329 +
1.330 +SPM *spm_read_hbm(const char *fname)
1.331 +{ SPM *A = NULL;
1.332 + HBM *hbm;
1.333 + int nrow, ncol, nnzero, i, j, beg, end, ptr, *colptr, *rowind;
1.334 + double val, *values;
1.335 + char *mxtype;
1.336 + hbm = hbm_read_mat(fname);
1.337 + if (hbm == NULL)
1.338 + { xprintf("spm_read_hbm: unable to read matrix\n");
1.339 + goto fini;
1.340 + }
1.341 + mxtype = hbm->mxtype;
1.342 + nrow = hbm->nrow;
1.343 + ncol = hbm->ncol;
1.344 + nnzero = hbm->nnzero;
1.345 + colptr = hbm->colptr;
1.346 + rowind = hbm->rowind;
1.347 + values = hbm->values;
1.348 + if (!(strcmp(mxtype, "RSA") == 0 || strcmp(mxtype, "PSA") == 0 ||
1.349 + strcmp(mxtype, "RUA") == 0 || strcmp(mxtype, "PUA") == 0 ||
1.350 + strcmp(mxtype, "RRA") == 0 || strcmp(mxtype, "PRA") == 0))
1.351 + { xprintf("spm_read_hbm: matrix type `%s' not supported\n",
1.352 + mxtype);
1.353 + goto fini;
1.354 + }
1.355 + A = spm_create_mat(nrow, ncol);
1.356 + if (mxtype[1] == 'S' || mxtype[1] == 'U')
1.357 + xassert(nrow == ncol);
1.358 + for (j = 1; j <= ncol; j++)
1.359 + { beg = colptr[j];
1.360 + end = colptr[j+1];
1.361 + xassert(1 <= beg && beg <= end && end <= nnzero + 1);
1.362 + for (ptr = beg; ptr < end; ptr++)
1.363 + { i = rowind[ptr];
1.364 + xassert(1 <= i && i <= nrow);
1.365 + if (mxtype[0] == 'R')
1.366 + val = values[ptr];
1.367 + else
1.368 + val = 1.0;
1.369 + spm_new_elem(A, i, j, val);
1.370 + if (mxtype[1] == 'S' && i != j)
1.371 + spm_new_elem(A, j, i, val);
1.372 + }
1.373 + }
1.374 +fini: if (hbm != NULL) hbm_free_mat(hbm);
1.375 + return A;
1.376 +}
1.377 +
1.378 +/***********************************************************************
1.379 +* NAME
1.380 +*
1.381 +* spm_count_nnz - determine number of non-zeros in sparse matrix
1.382 +*
1.383 +* SYNOPSIS
1.384 +*
1.385 +* #include "glpspm.h"
1.386 +* int spm_count_nnz(const SPM *A);
1.387 +*
1.388 +* RETURNS
1.389 +*
1.390 +* The routine spm_count_nnz returns the number of structural non-zero
1.391 +* elements in the specified sparse matrix. */
1.392 +
1.393 +int spm_count_nnz(const SPM *A)
1.394 +{ SPME *e;
1.395 + int i, nnz = 0;
1.396 + for (i = 1; i <= A->m; i++)
1.397 + for (e = A->row[i]; e != NULL; e = e->r_next) nnz++;
1.398 + return nnz;
1.399 +}
1.400 +
1.401 +/***********************************************************************
1.402 +* NAME
1.403 +*
1.404 +* spm_drop_zeros - remove zero elements from sparse matrix
1.405 +*
1.406 +* SYNOPSIS
1.407 +*
1.408 +* #include "glpspm.h"
1.409 +* int spm_drop_zeros(SPM *A, double eps);
1.410 +*
1.411 +* DESCRIPTION
1.412 +*
1.413 +* The routine spm_drop_zeros removes all elements from the specified
1.414 +* sparse matrix, whose absolute value is less than eps.
1.415 +*
1.416 +* If the parameter eps is 0, only zero elements are removed from the
1.417 +* matrix.
1.418 +*
1.419 +* RETURNS
1.420 +*
1.421 +* The routine returns the number of elements removed. */
1.422 +
1.423 +int spm_drop_zeros(SPM *A, double eps)
1.424 +{ SPME *e, *next;
1.425 + int i, count = 0;
1.426 + for (i = 1; i <= A->m; i++)
1.427 + { for (e = A->row[i]; e != NULL; e = next)
1.428 + { next = e->r_next;
1.429 + if (e->val == 0.0 || fabs(e->val) < eps)
1.430 + { /* remove element from the row list */
1.431 + if (e->r_prev == NULL)
1.432 + A->row[e->i] = e->r_next;
1.433 + else
1.434 + e->r_prev->r_next = e->r_next;
1.435 + if (e->r_next == NULL)
1.436 + ;
1.437 + else
1.438 + e->r_next->r_prev = e->r_prev;
1.439 + /* remove element from the column list */
1.440 + if (e->c_prev == NULL)
1.441 + A->col[e->j] = e->c_next;
1.442 + else
1.443 + e->c_prev->c_next = e->c_next;
1.444 + if (e->c_next == NULL)
1.445 + ;
1.446 + else
1.447 + e->c_next->c_prev = e->c_prev;
1.448 + /* return element to the memory pool */
1.449 + dmp_free_atom(A->pool, e, sizeof(SPME));
1.450 + count++;
1.451 + }
1.452 + }
1.453 + }
1.454 + return count;
1.455 +}
1.456 +
1.457 +/***********************************************************************
1.458 +* NAME
1.459 +*
1.460 +* spm_read_mat - read sparse matrix from text file
1.461 +*
1.462 +* SYNOPSIS
1.463 +*
1.464 +* #include "glpspm.h"
1.465 +* SPM *spm_read_mat(const char *fname);
1.466 +*
1.467 +* DESCRIPTION
1.468 +*
1.469 +* The routine reads a sparse matrix from a text file whose name is
1.470 +* specified by the parameter fname.
1.471 +*
1.472 +* For the file format see description of the routine spm_write_mat.
1.473 +*
1.474 +* RETURNS
1.475 +*
1.476 +* On success the routine returns a pointer to the matrix created,
1.477 +* otherwise NULL. */
1.478 +
1.479 +#if 1
1.480 +SPM *spm_read_mat(const char *fname)
1.481 +{ xassert(fname != fname);
1.482 + return NULL;
1.483 +}
1.484 +#else
1.485 +SPM *spm_read_mat(const char *fname)
1.486 +{ SPM *A = NULL;
1.487 + PDS *pds;
1.488 + jmp_buf jump;
1.489 + int i, j, k, m, n, nnz, fail = 0;
1.490 + double val;
1.491 + xprintf("spm_read_mat: reading matrix from `%s'...\n", fname);
1.492 + pds = pds_open_file(fname);
1.493 + if (pds == NULL)
1.494 + { xprintf("spm_read_mat: unable to open `%s' - %s\n", fname,
1.495 + strerror(errno));
1.496 + fail = 1;
1.497 + goto done;
1.498 + }
1.499 + if (setjmp(jump))
1.500 + { fail = 1;
1.501 + goto done;
1.502 + }
1.503 + pds_set_jump(pds, jump);
1.504 + /* number of rows, number of columns, number of non-zeros */
1.505 + m = pds_scan_int(pds);
1.506 + if (m < 0)
1.507 + pds_error(pds, "invalid number of rows\n");
1.508 + n = pds_scan_int(pds);
1.509 + if (n < 0)
1.510 + pds_error(pds, "invalid number of columns\n");
1.511 + nnz = pds_scan_int(pds);
1.512 + if (nnz < 0)
1.513 + pds_error(pds, "invalid number of non-zeros\n");
1.514 + /* create matrix */
1.515 + xprintf("spm_read_mat: %d rows, %d columns, %d non-zeros\n",
1.516 + m, n, nnz);
1.517 + A = spm_create_mat(m, n);
1.518 + /* read matrix elements */
1.519 + for (k = 1; k <= nnz; k++)
1.520 + { /* row index, column index, element value */
1.521 + i = pds_scan_int(pds);
1.522 + if (!(1 <= i && i <= m))
1.523 + pds_error(pds, "row index out of range\n");
1.524 + j = pds_scan_int(pds);
1.525 + if (!(1 <= j && j <= n))
1.526 + pds_error(pds, "column index out of range\n");
1.527 + val = pds_scan_num(pds);
1.528 + /* add new element to the matrix */
1.529 + spm_new_elem(A, i, j, val);
1.530 + }
1.531 + xprintf("spm_read_mat: %d lines were read\n", pds->count);
1.532 +done: if (pds != NULL) pds_close_file(pds);
1.533 + if (fail && A != NULL) spm_delete_mat(A), A = NULL;
1.534 + return A;
1.535 +}
1.536 +#endif
1.537 +
1.538 +/***********************************************************************
1.539 +* NAME
1.540 +*
1.541 +* spm_write_mat - write sparse matrix to text file
1.542 +*
1.543 +* SYNOPSIS
1.544 +*
1.545 +* #include "glpspm.h"
1.546 +* int spm_write_mat(const SPM *A, const char *fname);
1.547 +*
1.548 +* DESCRIPTION
1.549 +*
1.550 +* The routine spm_write_mat writes the specified sparse matrix to a
1.551 +* text file whose name is specified by the parameter fname. This file
1.552 +* can be read back with the routine spm_read_mat.
1.553 +*
1.554 +* RETURNS
1.555 +*
1.556 +* On success the routine returns zero, otherwise non-zero.
1.557 +*
1.558 +* FILE FORMAT
1.559 +*
1.560 +* The file created by the routine spm_write_mat is a plain text file,
1.561 +* which contains the following information:
1.562 +*
1.563 +* m n nnz
1.564 +* row[1] col[1] val[1]
1.565 +* row[2] col[2] val[2]
1.566 +* . . .
1.567 +* row[nnz] col[nnz] val[nnz]
1.568 +*
1.569 +* where:
1.570 +* m is the number of rows;
1.571 +* n is the number of columns;
1.572 +* nnz is the number of non-zeros;
1.573 +* row[k], k = 1,...,nnz, are row indices;
1.574 +* col[k], k = 1,...,nnz, are column indices;
1.575 +* val[k], k = 1,...,nnz, are element values. */
1.576 +
1.577 +#if 1
1.578 +int spm_write_mat(const SPM *A, const char *fname)
1.579 +{ xassert(A != A);
1.580 + xassert(fname != fname);
1.581 + return 0;
1.582 +}
1.583 +#else
1.584 +int spm_write_mat(const SPM *A, const char *fname)
1.585 +{ FILE *fp;
1.586 + int i, nnz, ret = 0;
1.587 + xprintf("spm_write_mat: writing matrix to `%s'...\n", fname);
1.588 + fp = fopen(fname, "w");
1.589 + if (fp == NULL)
1.590 + { xprintf("spm_write_mat: unable to create `%s' - %s\n", fname,
1.591 + strerror(errno));
1.592 + ret = 1;
1.593 + goto done;
1.594 + }
1.595 + /* number of rows, number of columns, number of non-zeros */
1.596 + nnz = spm_count_nnz(A);
1.597 + fprintf(fp, "%d %d %d\n", A->m, A->n, nnz);
1.598 + /* walk through rows of the matrix */
1.599 + for (i = 1; i <= A->m; i++)
1.600 + { SPME *e;
1.601 + /* walk through elements of i-th row */
1.602 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.603 + { /* row index, column index, element value */
1.604 + fprintf(fp, "%d %d %.*g\n", e->i, e->j, DBL_DIG, e->val);
1.605 + }
1.606 + }
1.607 + fflush(fp);
1.608 + if (ferror(fp))
1.609 + { xprintf("spm_write_mat: writing error on `%s' - %s\n", fname,
1.610 + strerror(errno));
1.611 + ret = 1;
1.612 + goto done;
1.613 + }
1.614 + xprintf("spm_write_mat: %d lines were written\n", 1 + nnz);
1.615 +done: if (fp != NULL) fclose(fp);
1.616 + return ret;
1.617 +}
1.618 +#endif
1.619 +
1.620 +/***********************************************************************
1.621 +* NAME
1.622 +*
1.623 +* spm_transpose - transpose sparse matrix
1.624 +*
1.625 +* SYNOPSIS
1.626 +*
1.627 +* #include "glpspm.h"
1.628 +* SPM *spm_transpose(const SPM *A);
1.629 +*
1.630 +* RETURNS
1.631 +*
1.632 +* The routine computes and returns sparse matrix B, which is a matrix
1.633 +* transposed to sparse matrix A. */
1.634 +
1.635 +SPM *spm_transpose(const SPM *A)
1.636 +{ SPM *B;
1.637 + int i;
1.638 + B = spm_create_mat(A->n, A->m);
1.639 + for (i = 1; i <= A->m; i++)
1.640 + { SPME *e;
1.641 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.642 + spm_new_elem(B, e->j, i, e->val);
1.643 + }
1.644 + return B;
1.645 +}
1.646 +
1.647 +SPM *spm_add_sym(const SPM *A, const SPM *B)
1.648 +{ /* add two sparse matrices (symbolic phase) */
1.649 + SPM *C;
1.650 + int i, j, *flag;
1.651 + xassert(A->m == B->m);
1.652 + xassert(A->n == B->n);
1.653 + /* create resultant matrix */
1.654 + C = spm_create_mat(A->m, A->n);
1.655 + /* allocate and clear the flag array */
1.656 + flag = xcalloc(1+C->n, sizeof(int));
1.657 + for (j = 1; j <= C->n; j++)
1.658 + flag[j] = 0;
1.659 + /* compute pattern of C = A + B */
1.660 + for (i = 1; i <= C->m; i++)
1.661 + { SPME *e;
1.662 + /* at the beginning i-th row of C is empty */
1.663 + /* (i-th row of C) := (i-th row of C) union (i-th row of A) */
1.664 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.665 + { /* (note that i-th row of A may have duplicate elements) */
1.666 + j = e->j;
1.667 + if (!flag[j])
1.668 + { spm_new_elem(C, i, j, 0.0);
1.669 + flag[j] = 1;
1.670 + }
1.671 + }
1.672 + /* (i-th row of C) := (i-th row of C) union (i-th row of B) */
1.673 + for (e = B->row[i]; e != NULL; e = e->r_next)
1.674 + { /* (note that i-th row of B may have duplicate elements) */
1.675 + j = e->j;
1.676 + if (!flag[j])
1.677 + { spm_new_elem(C, i, j, 0.0);
1.678 + flag[j] = 1;
1.679 + }
1.680 + }
1.681 + /* reset the flag array */
1.682 + for (e = C->row[i]; e != NULL; e = e->r_next)
1.683 + flag[e->j] = 0;
1.684 + }
1.685 + /* check and deallocate the flag array */
1.686 + for (j = 1; j <= C->n; j++)
1.687 + xassert(!flag[j]);
1.688 + xfree(flag);
1.689 + return C;
1.690 +}
1.691 +
1.692 +void spm_add_num(SPM *C, double alfa, const SPM *A, double beta,
1.693 + const SPM *B)
1.694 +{ /* add two sparse matrices (numeric phase) */
1.695 + int i, j;
1.696 + double *work;
1.697 + /* allocate and clear the working array */
1.698 + work = xcalloc(1+C->n, sizeof(double));
1.699 + for (j = 1; j <= C->n; j++)
1.700 + work[j] = 0.0;
1.701 + /* compute matrix C = alfa * A + beta * B */
1.702 + for (i = 1; i <= C->n; i++)
1.703 + { SPME *e;
1.704 + /* work := alfa * (i-th row of A) + beta * (i-th row of B) */
1.705 + /* (note that A and/or B may have duplicate elements) */
1.706 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.707 + work[e->j] += alfa * e->val;
1.708 + for (e = B->row[i]; e != NULL; e = e->r_next)
1.709 + work[e->j] += beta * e->val;
1.710 + /* (i-th row of C) := work, work := 0 */
1.711 + for (e = C->row[i]; e != NULL; e = e->r_next)
1.712 + { j = e->j;
1.713 + e->val = work[j];
1.714 + work[j] = 0.0;
1.715 + }
1.716 + }
1.717 + /* check and deallocate the working array */
1.718 + for (j = 1; j <= C->n; j++)
1.719 + xassert(work[j] == 0.0);
1.720 + xfree(work);
1.721 + return;
1.722 +}
1.723 +
1.724 +SPM *spm_add_mat(double alfa, const SPM *A, double beta, const SPM *B)
1.725 +{ /* add two sparse matrices (driver routine) */
1.726 + SPM *C;
1.727 + C = spm_add_sym(A, B);
1.728 + spm_add_num(C, alfa, A, beta, B);
1.729 + return C;
1.730 +}
1.731 +
1.732 +SPM *spm_mul_sym(const SPM *A, const SPM *B)
1.733 +{ /* multiply two sparse matrices (symbolic phase) */
1.734 + int i, j, k, *flag;
1.735 + SPM *C;
1.736 + xassert(A->n == B->m);
1.737 + /* create resultant matrix */
1.738 + C = spm_create_mat(A->m, B->n);
1.739 + /* allocate and clear the flag array */
1.740 + flag = xcalloc(1+C->n, sizeof(int));
1.741 + for (j = 1; j <= C->n; j++)
1.742 + flag[j] = 0;
1.743 + /* compute pattern of C = A * B */
1.744 + for (i = 1; i <= C->m; i++)
1.745 + { SPME *e, *ee;
1.746 + /* compute pattern of i-th row of C */
1.747 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.748 + { k = e->j;
1.749 + for (ee = B->row[k]; ee != NULL; ee = ee->r_next)
1.750 + { j = ee->j;
1.751 + /* if a[i,k] != 0 and b[k,j] != 0 then c[i,j] != 0 */
1.752 + if (!flag[j])
1.753 + { /* c[i,j] does not exist, so create it */
1.754 + spm_new_elem(C, i, j, 0.0);
1.755 + flag[j] = 1;
1.756 + }
1.757 + }
1.758 + }
1.759 + /* reset the flag array */
1.760 + for (e = C->row[i]; e != NULL; e = e->r_next)
1.761 + flag[e->j] = 0;
1.762 + }
1.763 + /* check and deallocate the flag array */
1.764 + for (j = 1; j <= C->n; j++)
1.765 + xassert(!flag[j]);
1.766 + xfree(flag);
1.767 + return C;
1.768 +}
1.769 +
1.770 +void spm_mul_num(SPM *C, const SPM *A, const SPM *B)
1.771 +{ /* multiply two sparse matrices (numeric phase) */
1.772 + int i, j;
1.773 + double *work;
1.774 + /* allocate and clear the working array */
1.775 + work = xcalloc(1+A->n, sizeof(double));
1.776 + for (j = 1; j <= A->n; j++)
1.777 + work[j] = 0.0;
1.778 + /* compute matrix C = A * B */
1.779 + for (i = 1; i <= C->m; i++)
1.780 + { SPME *e, *ee;
1.781 + double temp;
1.782 + /* work := (i-th row of A) */
1.783 + /* (note that A may have duplicate elements) */
1.784 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.785 + work[e->j] += e->val;
1.786 + /* compute i-th row of C */
1.787 + for (e = C->row[i]; e != NULL; e = e->r_next)
1.788 + { j = e->j;
1.789 + /* c[i,j] := work * (j-th column of B) */
1.790 + temp = 0.0;
1.791 + for (ee = B->col[j]; ee != NULL; ee = ee->c_next)
1.792 + temp += work[ee->i] * ee->val;
1.793 + e->val = temp;
1.794 + }
1.795 + /* reset the working array */
1.796 + for (e = A->row[i]; e != NULL; e = e->r_next)
1.797 + work[e->j] = 0.0;
1.798 + }
1.799 + /* check and deallocate the working array */
1.800 + for (j = 1; j <= A->n; j++)
1.801 + xassert(work[j] == 0.0);
1.802 + xfree(work);
1.803 + return;
1.804 +}
1.805 +
1.806 +SPM *spm_mul_mat(const SPM *A, const SPM *B)
1.807 +{ /* multiply two sparse matrices (driver routine) */
1.808 + SPM *C;
1.809 + C = spm_mul_sym(A, B);
1.810 + spm_mul_num(C, A, B);
1.811 + return C;
1.812 +}
1.813 +
1.814 +PER *spm_create_per(int n)
1.815 +{ /* create permutation matrix */
1.816 + PER *P;
1.817 + int k;
1.818 + xassert(n >= 0);
1.819 + P = xmalloc(sizeof(PER));
1.820 + P->n = n;
1.821 + P->row = xcalloc(1+n, sizeof(int));
1.822 + P->col = xcalloc(1+n, sizeof(int));
1.823 + /* initially it is identity matrix */
1.824 + for (k = 1; k <= n; k++)
1.825 + P->row[k] = P->col[k] = k;
1.826 + return P;
1.827 +}
1.828 +
1.829 +void spm_check_per(PER *P)
1.830 +{ /* check permutation matrix for correctness */
1.831 + int i, j;
1.832 + xassert(P->n >= 0);
1.833 + for (i = 1; i <= P->n; i++)
1.834 + { j = P->row[i];
1.835 + xassert(1 <= j && j <= P->n);
1.836 + xassert(P->col[j] == i);
1.837 + }
1.838 + return;
1.839 +}
1.840 +
1.841 +void spm_delete_per(PER *P)
1.842 +{ /* delete permutation matrix */
1.843 + xfree(P->row);
1.844 + xfree(P->col);
1.845 + xfree(P);
1.846 + return;
1.847 +}
1.848 +
1.849 +/* eof */