[9] | 1 | /* glpspm.c */ |
---|
| 2 | |
---|
| 3 | /*********************************************************************** |
---|
| 4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
| 5 | * |
---|
| 6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
| 7 | * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, |
---|
| 8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
| 9 | * E-mail: <mao@gnu.org>. |
---|
| 10 | * |
---|
| 11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
| 12 | * under the terms of the GNU General Public License as published by |
---|
| 13 | * the Free Software Foundation, either version 3 of the License, or |
---|
| 14 | * (at your option) any later version. |
---|
| 15 | * |
---|
| 16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
| 17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
| 18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
| 19 | * License for more details. |
---|
| 20 | * |
---|
| 21 | * You should have received a copy of the GNU General Public License |
---|
| 22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
| 23 | ***********************************************************************/ |
---|
| 24 | |
---|
| 25 | #include "glphbm.h" |
---|
| 26 | #include "glprgr.h" |
---|
| 27 | #include "glpspm.h" |
---|
| 28 | |
---|
| 29 | /*********************************************************************** |
---|
| 30 | * NAME |
---|
| 31 | * |
---|
| 32 | * spm_create_mat - create general sparse matrix |
---|
| 33 | * |
---|
| 34 | * SYNOPSIS |
---|
| 35 | * |
---|
| 36 | * #include "glpspm.h" |
---|
| 37 | * SPM *spm_create_mat(int m, int n); |
---|
| 38 | * |
---|
| 39 | * DESCRIPTION |
---|
| 40 | * |
---|
| 41 | * The routine spm_create_mat creates a general sparse matrix having |
---|
| 42 | * m rows and n columns. Being created the matrix is zero (empty), i.e. |
---|
| 43 | * has no elements. |
---|
| 44 | * |
---|
| 45 | * RETURNS |
---|
| 46 | * |
---|
| 47 | * The routine returns a pointer to the matrix created. */ |
---|
| 48 | |
---|
| 49 | SPM *spm_create_mat(int m, int n) |
---|
| 50 | { SPM *A; |
---|
| 51 | xassert(0 <= m && m < INT_MAX); |
---|
| 52 | xassert(0 <= n && n < INT_MAX); |
---|
| 53 | A = xmalloc(sizeof(SPM)); |
---|
| 54 | A->m = m; |
---|
| 55 | A->n = n; |
---|
| 56 | if (m == 0 || n == 0) |
---|
| 57 | { A->pool = NULL; |
---|
| 58 | A->row = NULL; |
---|
| 59 | A->col = NULL; |
---|
| 60 | } |
---|
| 61 | else |
---|
| 62 | { int i, j; |
---|
| 63 | A->pool = dmp_create_pool(); |
---|
| 64 | A->row = xcalloc(1+m, sizeof(SPME *)); |
---|
| 65 | for (i = 1; i <= m; i++) A->row[i] = NULL; |
---|
| 66 | A->col = xcalloc(1+n, sizeof(SPME *)); |
---|
| 67 | for (j = 1; j <= n; j++) A->col[j] = NULL; |
---|
| 68 | } |
---|
| 69 | return A; |
---|
| 70 | } |
---|
| 71 | |
---|
| 72 | /*********************************************************************** |
---|
| 73 | * NAME |
---|
| 74 | * |
---|
| 75 | * spm_new_elem - add new element to sparse matrix |
---|
| 76 | * |
---|
| 77 | * SYNOPSIS |
---|
| 78 | * |
---|
| 79 | * #include "glpspm.h" |
---|
| 80 | * SPME *spm_new_elem(SPM *A, int i, int j, double val); |
---|
| 81 | * |
---|
| 82 | * DESCRIPTION |
---|
| 83 | * |
---|
| 84 | * The routine spm_new_elem adds a new element to the specified sparse |
---|
| 85 | * matrix. Parameters i, j, and val specify the row number, the column |
---|
| 86 | * number, and a numerical value of the element, respectively. |
---|
| 87 | * |
---|
| 88 | * RETURNS |
---|
| 89 | * |
---|
| 90 | * The routine returns a pointer to the new element added. */ |
---|
| 91 | |
---|
| 92 | SPME *spm_new_elem(SPM *A, int i, int j, double val) |
---|
| 93 | { SPME *e; |
---|
| 94 | xassert(1 <= i && i <= A->m); |
---|
| 95 | xassert(1 <= j && j <= A->n); |
---|
| 96 | e = dmp_get_atom(A->pool, sizeof(SPME)); |
---|
| 97 | e->i = i; |
---|
| 98 | e->j = j; |
---|
| 99 | e->val = val; |
---|
| 100 | e->r_prev = NULL; |
---|
| 101 | e->r_next = A->row[i]; |
---|
| 102 | if (e->r_next != NULL) e->r_next->r_prev = e; |
---|
| 103 | e->c_prev = NULL; |
---|
| 104 | e->c_next = A->col[j]; |
---|
| 105 | if (e->c_next != NULL) e->c_next->c_prev = e; |
---|
| 106 | A->row[i] = A->col[j] = e; |
---|
| 107 | return e; |
---|
| 108 | } |
---|
| 109 | |
---|
| 110 | /*********************************************************************** |
---|
| 111 | * NAME |
---|
| 112 | * |
---|
| 113 | * spm_delete_mat - delete general sparse matrix |
---|
| 114 | * |
---|
| 115 | * SYNOPSIS |
---|
| 116 | * |
---|
| 117 | * #include "glpspm.h" |
---|
| 118 | * void spm_delete_mat(SPM *A); |
---|
| 119 | * |
---|
| 120 | * DESCRIPTION |
---|
| 121 | * |
---|
| 122 | * The routine deletes the specified general sparse matrix freeing all |
---|
| 123 | * the memory allocated to this object. */ |
---|
| 124 | |
---|
| 125 | void spm_delete_mat(SPM *A) |
---|
| 126 | { /* delete sparse matrix */ |
---|
| 127 | if (A->pool != NULL) dmp_delete_pool(A->pool); |
---|
| 128 | if (A->row != NULL) xfree(A->row); |
---|
| 129 | if (A->col != NULL) xfree(A->col); |
---|
| 130 | xfree(A); |
---|
| 131 | return; |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | /*********************************************************************** |
---|
| 135 | * NAME |
---|
| 136 | * |
---|
| 137 | * spm_test_mat_e - create test sparse matrix of E(n,c) class |
---|
| 138 | * |
---|
| 139 | * SYNOPSIS |
---|
| 140 | * |
---|
| 141 | * #include "glpspm.h" |
---|
| 142 | * SPM *spm_test_mat_e(int n, int c); |
---|
| 143 | * |
---|
| 144 | * DESCRIPTION |
---|
| 145 | * |
---|
| 146 | * The routine spm_test_mat_e creates a test sparse matrix of E(n,c) |
---|
| 147 | * class as described in the book: Ole 0sterby, Zahari Zlatev. Direct |
---|
| 148 | * Methods for Sparse Matrices. Springer-Verlag, 1983. |
---|
| 149 | * |
---|
| 150 | * Matrix of E(n,c) class is a symmetric positive definite matrix of |
---|
| 151 | * the order n. It has the number 4 on its main diagonal and the number |
---|
| 152 | * -1 on its four co-diagonals, two of which are neighbour to the main |
---|
| 153 | * diagonal and two others are shifted from the main diagonal on the |
---|
| 154 | * distance c. |
---|
| 155 | * |
---|
| 156 | * It is necessary that n >= 3 and 2 <= c <= n-1. |
---|
| 157 | * |
---|
| 158 | * RETURNS |
---|
| 159 | * |
---|
| 160 | * The routine returns a pointer to the matrix created. */ |
---|
| 161 | |
---|
| 162 | SPM *spm_test_mat_e(int n, int c) |
---|
| 163 | { SPM *A; |
---|
| 164 | int i; |
---|
| 165 | xassert(n >= 3 && 2 <= c && c <= n-1); |
---|
| 166 | A = spm_create_mat(n, n); |
---|
| 167 | for (i = 1; i <= n; i++) |
---|
| 168 | spm_new_elem(A, i, i, 4.0); |
---|
| 169 | for (i = 1; i <= n-1; i++) |
---|
| 170 | { spm_new_elem(A, i, i+1, -1.0); |
---|
| 171 | spm_new_elem(A, i+1, i, -1.0); |
---|
| 172 | } |
---|
| 173 | for (i = 1; i <= n-c; i++) |
---|
| 174 | { spm_new_elem(A, i, i+c, -1.0); |
---|
| 175 | spm_new_elem(A, i+c, i, -1.0); |
---|
| 176 | } |
---|
| 177 | return A; |
---|
| 178 | } |
---|
| 179 | |
---|
| 180 | /*********************************************************************** |
---|
| 181 | * NAME |
---|
| 182 | * |
---|
| 183 | * spm_test_mat_d - create test sparse matrix of D(n,c) class |
---|
| 184 | * |
---|
| 185 | * SYNOPSIS |
---|
| 186 | * |
---|
| 187 | * #include "glpspm.h" |
---|
| 188 | * SPM *spm_test_mat_d(int n, int c); |
---|
| 189 | * |
---|
| 190 | * DESCRIPTION |
---|
| 191 | * |
---|
| 192 | * The routine spm_test_mat_d creates a test sparse matrix of D(n,c) |
---|
| 193 | * class as described in the book: Ole 0sterby, Zahari Zlatev. Direct |
---|
| 194 | * Methods for Sparse Matrices. Springer-Verlag, 1983. |
---|
| 195 | * |
---|
| 196 | * Matrix of D(n,c) class is a non-singular matrix of the order n. It |
---|
| 197 | * has unity main diagonal, three co-diagonals above the main diagonal |
---|
| 198 | * on the distance c, which are cyclically continued below the main |
---|
| 199 | * diagonal, and a triangle block of the size 10x10 in the upper right |
---|
| 200 | * corner. |
---|
| 201 | * |
---|
| 202 | * It is necessary that n >= 14 and 1 <= c <= n-13. |
---|
| 203 | * |
---|
| 204 | * RETURNS |
---|
| 205 | * |
---|
| 206 | * The routine returns a pointer to the matrix created. */ |
---|
| 207 | |
---|
| 208 | SPM *spm_test_mat_d(int n, int c) |
---|
| 209 | { SPM *A; |
---|
| 210 | int i, j; |
---|
| 211 | xassert(n >= 14 && 1 <= c && c <= n-13); |
---|
| 212 | A = spm_create_mat(n, n); |
---|
| 213 | for (i = 1; i <= n; i++) |
---|
| 214 | spm_new_elem(A, i, i, 1.0); |
---|
| 215 | for (i = 1; i <= n-c; i++) |
---|
| 216 | spm_new_elem(A, i, i+c, (double)(i+1)); |
---|
| 217 | for (i = n-c+1; i <= n; i++) |
---|
| 218 | spm_new_elem(A, i, i-n+c, (double)(i+1)); |
---|
| 219 | for (i = 1; i <= n-c-1; i++) |
---|
| 220 | spm_new_elem(A, i, i+c+1, (double)(-i)); |
---|
| 221 | for (i = n-c; i <= n; i++) |
---|
| 222 | spm_new_elem(A, i, i-n+c+1, (double)(-i)); |
---|
| 223 | for (i = 1; i <= n-c-2; i++) |
---|
| 224 | spm_new_elem(A, i, i+c+2, 16.0); |
---|
| 225 | for (i = n-c-1; i <= n; i++) |
---|
| 226 | spm_new_elem(A, i, i-n+c+2, 16.0); |
---|
| 227 | for (j = 1; j <= 10; j++) |
---|
| 228 | for (i = 1; i <= 11-j; i++) |
---|
| 229 | spm_new_elem(A, i, n-11+i+j, 100.0 * (double)j); |
---|
| 230 | return A; |
---|
| 231 | } |
---|
| 232 | |
---|
| 233 | /*********************************************************************** |
---|
| 234 | * NAME |
---|
| 235 | * |
---|
| 236 | * spm_show_mat - write sparse matrix pattern in BMP file format |
---|
| 237 | * |
---|
| 238 | * SYNOPSIS |
---|
| 239 | * |
---|
| 240 | * #include "glpspm.h" |
---|
| 241 | * int spm_show_mat(const SPM *A, const char *fname); |
---|
| 242 | * |
---|
| 243 | * DESCRIPTION |
---|
| 244 | * |
---|
| 245 | * The routine spm_show_mat writes pattern of the specified sparse |
---|
| 246 | * matrix in uncompressed BMP file format (Windows bitmap) to a binary |
---|
| 247 | * file whose name is specified by the character string fname. |
---|
| 248 | * |
---|
| 249 | * Each pixel corresponds to one matrix element. The pixel colors have |
---|
| 250 | * the following meaning: |
---|
| 251 | * |
---|
| 252 | * Black structurally zero element |
---|
| 253 | * White positive element |
---|
| 254 | * Cyan negative element |
---|
| 255 | * Green zero element |
---|
| 256 | * Red duplicate element |
---|
| 257 | * |
---|
| 258 | * RETURNS |
---|
| 259 | * |
---|
| 260 | * If no error occured, the routine returns zero. Otherwise, it prints |
---|
| 261 | * an appropriate error message and returns non-zero. */ |
---|
| 262 | |
---|
| 263 | int spm_show_mat(const SPM *A, const char *fname) |
---|
| 264 | { int m = A->m; |
---|
| 265 | int n = A->n; |
---|
| 266 | int i, j, k, ret; |
---|
| 267 | char *map; |
---|
| 268 | xprintf("spm_show_mat: writing matrix pattern to `%s'...\n", |
---|
| 269 | fname); |
---|
| 270 | xassert(1 <= m && m <= 32767); |
---|
| 271 | xassert(1 <= n && n <= 32767); |
---|
| 272 | map = xmalloc(m * n); |
---|
| 273 | memset(map, 0x08, m * n); |
---|
| 274 | for (i = 1; i <= m; i++) |
---|
| 275 | { SPME *e; |
---|
| 276 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 277 | { j = e->j; |
---|
| 278 | xassert(1 <= j && j <= n); |
---|
| 279 | k = n * (i - 1) + (j - 1); |
---|
| 280 | if (map[k] != 0x08) |
---|
| 281 | map[k] = 0x0C; |
---|
| 282 | else if (e->val > 0.0) |
---|
| 283 | map[k] = 0x0F; |
---|
| 284 | else if (e->val < 0.0) |
---|
| 285 | map[k] = 0x0B; |
---|
| 286 | else |
---|
| 287 | map[k] = 0x0A; |
---|
| 288 | } |
---|
| 289 | } |
---|
| 290 | ret = rgr_write_bmp16(fname, m, n, map); |
---|
| 291 | xfree(map); |
---|
| 292 | return ret; |
---|
| 293 | } |
---|
| 294 | |
---|
| 295 | /*********************************************************************** |
---|
| 296 | * NAME |
---|
| 297 | * |
---|
| 298 | * spm_read_hbm - read sparse matrix in Harwell-Boeing format |
---|
| 299 | * |
---|
| 300 | * SYNOPSIS |
---|
| 301 | * |
---|
| 302 | * #include "glpspm.h" |
---|
| 303 | * SPM *spm_read_hbm(const char *fname); |
---|
| 304 | * |
---|
| 305 | * DESCRIPTION |
---|
| 306 | * |
---|
| 307 | * The routine spm_read_hbm reads a sparse matrix in the Harwell-Boeing |
---|
| 308 | * format from a text file whose name is the character string fname. |
---|
| 309 | * |
---|
| 310 | * Detailed description of the Harwell-Boeing format recognised by this |
---|
| 311 | * routine can be found in the following report: |
---|
| 312 | * |
---|
| 313 | * I.S.Duff, R.G.Grimes, J.G.Lewis. User's Guide for the Harwell-Boeing |
---|
| 314 | * Sparse Matrix Collection (Release I), TR/PA/92/86, October 1992. |
---|
| 315 | * |
---|
| 316 | * NOTE |
---|
| 317 | * |
---|
| 318 | * The routine spm_read_hbm reads the matrix "as is", due to which zero |
---|
| 319 | * and/or duplicate elements can appear in the matrix. |
---|
| 320 | * |
---|
| 321 | * RETURNS |
---|
| 322 | * |
---|
| 323 | * If no error occured, the routine returns a pointer to the matrix |
---|
| 324 | * created. Otherwise, the routine prints an appropriate error message |
---|
| 325 | * and returns NULL. */ |
---|
| 326 | |
---|
| 327 | SPM *spm_read_hbm(const char *fname) |
---|
| 328 | { SPM *A = NULL; |
---|
| 329 | HBM *hbm; |
---|
| 330 | int nrow, ncol, nnzero, i, j, beg, end, ptr, *colptr, *rowind; |
---|
| 331 | double val, *values; |
---|
| 332 | char *mxtype; |
---|
| 333 | hbm = hbm_read_mat(fname); |
---|
| 334 | if (hbm == NULL) |
---|
| 335 | { xprintf("spm_read_hbm: unable to read matrix\n"); |
---|
| 336 | goto fini; |
---|
| 337 | } |
---|
| 338 | mxtype = hbm->mxtype; |
---|
| 339 | nrow = hbm->nrow; |
---|
| 340 | ncol = hbm->ncol; |
---|
| 341 | nnzero = hbm->nnzero; |
---|
| 342 | colptr = hbm->colptr; |
---|
| 343 | rowind = hbm->rowind; |
---|
| 344 | values = hbm->values; |
---|
| 345 | if (!(strcmp(mxtype, "RSA") == 0 || strcmp(mxtype, "PSA") == 0 || |
---|
| 346 | strcmp(mxtype, "RUA") == 0 || strcmp(mxtype, "PUA") == 0 || |
---|
| 347 | strcmp(mxtype, "RRA") == 0 || strcmp(mxtype, "PRA") == 0)) |
---|
| 348 | { xprintf("spm_read_hbm: matrix type `%s' not supported\n", |
---|
| 349 | mxtype); |
---|
| 350 | goto fini; |
---|
| 351 | } |
---|
| 352 | A = spm_create_mat(nrow, ncol); |
---|
| 353 | if (mxtype[1] == 'S' || mxtype[1] == 'U') |
---|
| 354 | xassert(nrow == ncol); |
---|
| 355 | for (j = 1; j <= ncol; j++) |
---|
| 356 | { beg = colptr[j]; |
---|
| 357 | end = colptr[j+1]; |
---|
| 358 | xassert(1 <= beg && beg <= end && end <= nnzero + 1); |
---|
| 359 | for (ptr = beg; ptr < end; ptr++) |
---|
| 360 | { i = rowind[ptr]; |
---|
| 361 | xassert(1 <= i && i <= nrow); |
---|
| 362 | if (mxtype[0] == 'R') |
---|
| 363 | val = values[ptr]; |
---|
| 364 | else |
---|
| 365 | val = 1.0; |
---|
| 366 | spm_new_elem(A, i, j, val); |
---|
| 367 | if (mxtype[1] == 'S' && i != j) |
---|
| 368 | spm_new_elem(A, j, i, val); |
---|
| 369 | } |
---|
| 370 | } |
---|
| 371 | fini: if (hbm != NULL) hbm_free_mat(hbm); |
---|
| 372 | return A; |
---|
| 373 | } |
---|
| 374 | |
---|
| 375 | /*********************************************************************** |
---|
| 376 | * NAME |
---|
| 377 | * |
---|
| 378 | * spm_count_nnz - determine number of non-zeros in sparse matrix |
---|
| 379 | * |
---|
| 380 | * SYNOPSIS |
---|
| 381 | * |
---|
| 382 | * #include "glpspm.h" |
---|
| 383 | * int spm_count_nnz(const SPM *A); |
---|
| 384 | * |
---|
| 385 | * RETURNS |
---|
| 386 | * |
---|
| 387 | * The routine spm_count_nnz returns the number of structural non-zero |
---|
| 388 | * elements in the specified sparse matrix. */ |
---|
| 389 | |
---|
| 390 | int spm_count_nnz(const SPM *A) |
---|
| 391 | { SPME *e; |
---|
| 392 | int i, nnz = 0; |
---|
| 393 | for (i = 1; i <= A->m; i++) |
---|
| 394 | for (e = A->row[i]; e != NULL; e = e->r_next) nnz++; |
---|
| 395 | return nnz; |
---|
| 396 | } |
---|
| 397 | |
---|
| 398 | /*********************************************************************** |
---|
| 399 | * NAME |
---|
| 400 | * |
---|
| 401 | * spm_drop_zeros - remove zero elements from sparse matrix |
---|
| 402 | * |
---|
| 403 | * SYNOPSIS |
---|
| 404 | * |
---|
| 405 | * #include "glpspm.h" |
---|
| 406 | * int spm_drop_zeros(SPM *A, double eps); |
---|
| 407 | * |
---|
| 408 | * DESCRIPTION |
---|
| 409 | * |
---|
| 410 | * The routine spm_drop_zeros removes all elements from the specified |
---|
| 411 | * sparse matrix, whose absolute value is less than eps. |
---|
| 412 | * |
---|
| 413 | * If the parameter eps is 0, only zero elements are removed from the |
---|
| 414 | * matrix. |
---|
| 415 | * |
---|
| 416 | * RETURNS |
---|
| 417 | * |
---|
| 418 | * The routine returns the number of elements removed. */ |
---|
| 419 | |
---|
| 420 | int spm_drop_zeros(SPM *A, double eps) |
---|
| 421 | { SPME *e, *next; |
---|
| 422 | int i, count = 0; |
---|
| 423 | for (i = 1; i <= A->m; i++) |
---|
| 424 | { for (e = A->row[i]; e != NULL; e = next) |
---|
| 425 | { next = e->r_next; |
---|
| 426 | if (e->val == 0.0 || fabs(e->val) < eps) |
---|
| 427 | { /* remove element from the row list */ |
---|
| 428 | if (e->r_prev == NULL) |
---|
| 429 | A->row[e->i] = e->r_next; |
---|
| 430 | else |
---|
| 431 | e->r_prev->r_next = e->r_next; |
---|
| 432 | if (e->r_next == NULL) |
---|
| 433 | ; |
---|
| 434 | else |
---|
| 435 | e->r_next->r_prev = e->r_prev; |
---|
| 436 | /* remove element from the column list */ |
---|
| 437 | if (e->c_prev == NULL) |
---|
| 438 | A->col[e->j] = e->c_next; |
---|
| 439 | else |
---|
| 440 | e->c_prev->c_next = e->c_next; |
---|
| 441 | if (e->c_next == NULL) |
---|
| 442 | ; |
---|
| 443 | else |
---|
| 444 | e->c_next->c_prev = e->c_prev; |
---|
| 445 | /* return element to the memory pool */ |
---|
| 446 | dmp_free_atom(A->pool, e, sizeof(SPME)); |
---|
| 447 | count++; |
---|
| 448 | } |
---|
| 449 | } |
---|
| 450 | } |
---|
| 451 | return count; |
---|
| 452 | } |
---|
| 453 | |
---|
| 454 | /*********************************************************************** |
---|
| 455 | * NAME |
---|
| 456 | * |
---|
| 457 | * spm_read_mat - read sparse matrix from text file |
---|
| 458 | * |
---|
| 459 | * SYNOPSIS |
---|
| 460 | * |
---|
| 461 | * #include "glpspm.h" |
---|
| 462 | * SPM *spm_read_mat(const char *fname); |
---|
| 463 | * |
---|
| 464 | * DESCRIPTION |
---|
| 465 | * |
---|
| 466 | * The routine reads a sparse matrix from a text file whose name is |
---|
| 467 | * specified by the parameter fname. |
---|
| 468 | * |
---|
| 469 | * For the file format see description of the routine spm_write_mat. |
---|
| 470 | * |
---|
| 471 | * RETURNS |
---|
| 472 | * |
---|
| 473 | * On success the routine returns a pointer to the matrix created, |
---|
| 474 | * otherwise NULL. */ |
---|
| 475 | |
---|
| 476 | #if 1 |
---|
| 477 | SPM *spm_read_mat(const char *fname) |
---|
| 478 | { xassert(fname != fname); |
---|
| 479 | return NULL; |
---|
| 480 | } |
---|
| 481 | #else |
---|
| 482 | SPM *spm_read_mat(const char *fname) |
---|
| 483 | { SPM *A = NULL; |
---|
| 484 | PDS *pds; |
---|
| 485 | jmp_buf jump; |
---|
| 486 | int i, j, k, m, n, nnz, fail = 0; |
---|
| 487 | double val; |
---|
| 488 | xprintf("spm_read_mat: reading matrix from `%s'...\n", fname); |
---|
| 489 | pds = pds_open_file(fname); |
---|
| 490 | if (pds == NULL) |
---|
| 491 | { xprintf("spm_read_mat: unable to open `%s' - %s\n", fname, |
---|
| 492 | strerror(errno)); |
---|
| 493 | fail = 1; |
---|
| 494 | goto done; |
---|
| 495 | } |
---|
| 496 | if (setjmp(jump)) |
---|
| 497 | { fail = 1; |
---|
| 498 | goto done; |
---|
| 499 | } |
---|
| 500 | pds_set_jump(pds, jump); |
---|
| 501 | /* number of rows, number of columns, number of non-zeros */ |
---|
| 502 | m = pds_scan_int(pds); |
---|
| 503 | if (m < 0) |
---|
| 504 | pds_error(pds, "invalid number of rows\n"); |
---|
| 505 | n = pds_scan_int(pds); |
---|
| 506 | if (n < 0) |
---|
| 507 | pds_error(pds, "invalid number of columns\n"); |
---|
| 508 | nnz = pds_scan_int(pds); |
---|
| 509 | if (nnz < 0) |
---|
| 510 | pds_error(pds, "invalid number of non-zeros\n"); |
---|
| 511 | /* create matrix */ |
---|
| 512 | xprintf("spm_read_mat: %d rows, %d columns, %d non-zeros\n", |
---|
| 513 | m, n, nnz); |
---|
| 514 | A = spm_create_mat(m, n); |
---|
| 515 | /* read matrix elements */ |
---|
| 516 | for (k = 1; k <= nnz; k++) |
---|
| 517 | { /* row index, column index, element value */ |
---|
| 518 | i = pds_scan_int(pds); |
---|
| 519 | if (!(1 <= i && i <= m)) |
---|
| 520 | pds_error(pds, "row index out of range\n"); |
---|
| 521 | j = pds_scan_int(pds); |
---|
| 522 | if (!(1 <= j && j <= n)) |
---|
| 523 | pds_error(pds, "column index out of range\n"); |
---|
| 524 | val = pds_scan_num(pds); |
---|
| 525 | /* add new element to the matrix */ |
---|
| 526 | spm_new_elem(A, i, j, val); |
---|
| 527 | } |
---|
| 528 | xprintf("spm_read_mat: %d lines were read\n", pds->count); |
---|
| 529 | done: if (pds != NULL) pds_close_file(pds); |
---|
| 530 | if (fail && A != NULL) spm_delete_mat(A), A = NULL; |
---|
| 531 | return A; |
---|
| 532 | } |
---|
| 533 | #endif |
---|
| 534 | |
---|
| 535 | /*********************************************************************** |
---|
| 536 | * NAME |
---|
| 537 | * |
---|
| 538 | * spm_write_mat - write sparse matrix to text file |
---|
| 539 | * |
---|
| 540 | * SYNOPSIS |
---|
| 541 | * |
---|
| 542 | * #include "glpspm.h" |
---|
| 543 | * int spm_write_mat(const SPM *A, const char *fname); |
---|
| 544 | * |
---|
| 545 | * DESCRIPTION |
---|
| 546 | * |
---|
| 547 | * The routine spm_write_mat writes the specified sparse matrix to a |
---|
| 548 | * text file whose name is specified by the parameter fname. This file |
---|
| 549 | * can be read back with the routine spm_read_mat. |
---|
| 550 | * |
---|
| 551 | * RETURNS |
---|
| 552 | * |
---|
| 553 | * On success the routine returns zero, otherwise non-zero. |
---|
| 554 | * |
---|
| 555 | * FILE FORMAT |
---|
| 556 | * |
---|
| 557 | * The file created by the routine spm_write_mat is a plain text file, |
---|
| 558 | * which contains the following information: |
---|
| 559 | * |
---|
| 560 | * m n nnz |
---|
| 561 | * row[1] col[1] val[1] |
---|
| 562 | * row[2] col[2] val[2] |
---|
| 563 | * . . . |
---|
| 564 | * row[nnz] col[nnz] val[nnz] |
---|
| 565 | * |
---|
| 566 | * where: |
---|
| 567 | * m is the number of rows; |
---|
| 568 | * n is the number of columns; |
---|
| 569 | * nnz is the number of non-zeros; |
---|
| 570 | * row[k], k = 1,...,nnz, are row indices; |
---|
| 571 | * col[k], k = 1,...,nnz, are column indices; |
---|
| 572 | * val[k], k = 1,...,nnz, are element values. */ |
---|
| 573 | |
---|
| 574 | #if 1 |
---|
| 575 | int spm_write_mat(const SPM *A, const char *fname) |
---|
| 576 | { xassert(A != A); |
---|
| 577 | xassert(fname != fname); |
---|
| 578 | return 0; |
---|
| 579 | } |
---|
| 580 | #else |
---|
| 581 | int spm_write_mat(const SPM *A, const char *fname) |
---|
| 582 | { FILE *fp; |
---|
| 583 | int i, nnz, ret = 0; |
---|
| 584 | xprintf("spm_write_mat: writing matrix to `%s'...\n", fname); |
---|
| 585 | fp = fopen(fname, "w"); |
---|
| 586 | if (fp == NULL) |
---|
| 587 | { xprintf("spm_write_mat: unable to create `%s' - %s\n", fname, |
---|
| 588 | strerror(errno)); |
---|
| 589 | ret = 1; |
---|
| 590 | goto done; |
---|
| 591 | } |
---|
| 592 | /* number of rows, number of columns, number of non-zeros */ |
---|
| 593 | nnz = spm_count_nnz(A); |
---|
| 594 | fprintf(fp, "%d %d %d\n", A->m, A->n, nnz); |
---|
| 595 | /* walk through rows of the matrix */ |
---|
| 596 | for (i = 1; i <= A->m; i++) |
---|
| 597 | { SPME *e; |
---|
| 598 | /* walk through elements of i-th row */ |
---|
| 599 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 600 | { /* row index, column index, element value */ |
---|
| 601 | fprintf(fp, "%d %d %.*g\n", e->i, e->j, DBL_DIG, e->val); |
---|
| 602 | } |
---|
| 603 | } |
---|
| 604 | fflush(fp); |
---|
| 605 | if (ferror(fp)) |
---|
| 606 | { xprintf("spm_write_mat: writing error on `%s' - %s\n", fname, |
---|
| 607 | strerror(errno)); |
---|
| 608 | ret = 1; |
---|
| 609 | goto done; |
---|
| 610 | } |
---|
| 611 | xprintf("spm_write_mat: %d lines were written\n", 1 + nnz); |
---|
| 612 | done: if (fp != NULL) fclose(fp); |
---|
| 613 | return ret; |
---|
| 614 | } |
---|
| 615 | #endif |
---|
| 616 | |
---|
| 617 | /*********************************************************************** |
---|
| 618 | * NAME |
---|
| 619 | * |
---|
| 620 | * spm_transpose - transpose sparse matrix |
---|
| 621 | * |
---|
| 622 | * SYNOPSIS |
---|
| 623 | * |
---|
| 624 | * #include "glpspm.h" |
---|
| 625 | * SPM *spm_transpose(const SPM *A); |
---|
| 626 | * |
---|
| 627 | * RETURNS |
---|
| 628 | * |
---|
| 629 | * The routine computes and returns sparse matrix B, which is a matrix |
---|
| 630 | * transposed to sparse matrix A. */ |
---|
| 631 | |
---|
| 632 | SPM *spm_transpose(const SPM *A) |
---|
| 633 | { SPM *B; |
---|
| 634 | int i; |
---|
| 635 | B = spm_create_mat(A->n, A->m); |
---|
| 636 | for (i = 1; i <= A->m; i++) |
---|
| 637 | { SPME *e; |
---|
| 638 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 639 | spm_new_elem(B, e->j, i, e->val); |
---|
| 640 | } |
---|
| 641 | return B; |
---|
| 642 | } |
---|
| 643 | |
---|
| 644 | SPM *spm_add_sym(const SPM *A, const SPM *B) |
---|
| 645 | { /* add two sparse matrices (symbolic phase) */ |
---|
| 646 | SPM *C; |
---|
| 647 | int i, j, *flag; |
---|
| 648 | xassert(A->m == B->m); |
---|
| 649 | xassert(A->n == B->n); |
---|
| 650 | /* create resultant matrix */ |
---|
| 651 | C = spm_create_mat(A->m, A->n); |
---|
| 652 | /* allocate and clear the flag array */ |
---|
| 653 | flag = xcalloc(1+C->n, sizeof(int)); |
---|
| 654 | for (j = 1; j <= C->n; j++) |
---|
| 655 | flag[j] = 0; |
---|
| 656 | /* compute pattern of C = A + B */ |
---|
| 657 | for (i = 1; i <= C->m; i++) |
---|
| 658 | { SPME *e; |
---|
| 659 | /* at the beginning i-th row of C is empty */ |
---|
| 660 | /* (i-th row of C) := (i-th row of C) union (i-th row of A) */ |
---|
| 661 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 662 | { /* (note that i-th row of A may have duplicate elements) */ |
---|
| 663 | j = e->j; |
---|
| 664 | if (!flag[j]) |
---|
| 665 | { spm_new_elem(C, i, j, 0.0); |
---|
| 666 | flag[j] = 1; |
---|
| 667 | } |
---|
| 668 | } |
---|
| 669 | /* (i-th row of C) := (i-th row of C) union (i-th row of B) */ |
---|
| 670 | for (e = B->row[i]; e != NULL; e = e->r_next) |
---|
| 671 | { /* (note that i-th row of B may have duplicate elements) */ |
---|
| 672 | j = e->j; |
---|
| 673 | if (!flag[j]) |
---|
| 674 | { spm_new_elem(C, i, j, 0.0); |
---|
| 675 | flag[j] = 1; |
---|
| 676 | } |
---|
| 677 | } |
---|
| 678 | /* reset the flag array */ |
---|
| 679 | for (e = C->row[i]; e != NULL; e = e->r_next) |
---|
| 680 | flag[e->j] = 0; |
---|
| 681 | } |
---|
| 682 | /* check and deallocate the flag array */ |
---|
| 683 | for (j = 1; j <= C->n; j++) |
---|
| 684 | xassert(!flag[j]); |
---|
| 685 | xfree(flag); |
---|
| 686 | return C; |
---|
| 687 | } |
---|
| 688 | |
---|
| 689 | void spm_add_num(SPM *C, double alfa, const SPM *A, double beta, |
---|
| 690 | const SPM *B) |
---|
| 691 | { /* add two sparse matrices (numeric phase) */ |
---|
| 692 | int i, j; |
---|
| 693 | double *work; |
---|
| 694 | /* allocate and clear the working array */ |
---|
| 695 | work = xcalloc(1+C->n, sizeof(double)); |
---|
| 696 | for (j = 1; j <= C->n; j++) |
---|
| 697 | work[j] = 0.0; |
---|
| 698 | /* compute matrix C = alfa * A + beta * B */ |
---|
| 699 | for (i = 1; i <= C->n; i++) |
---|
| 700 | { SPME *e; |
---|
| 701 | /* work := alfa * (i-th row of A) + beta * (i-th row of B) */ |
---|
| 702 | /* (note that A and/or B may have duplicate elements) */ |
---|
| 703 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 704 | work[e->j] += alfa * e->val; |
---|
| 705 | for (e = B->row[i]; e != NULL; e = e->r_next) |
---|
| 706 | work[e->j] += beta * e->val; |
---|
| 707 | /* (i-th row of C) := work, work := 0 */ |
---|
| 708 | for (e = C->row[i]; e != NULL; e = e->r_next) |
---|
| 709 | { j = e->j; |
---|
| 710 | e->val = work[j]; |
---|
| 711 | work[j] = 0.0; |
---|
| 712 | } |
---|
| 713 | } |
---|
| 714 | /* check and deallocate the working array */ |
---|
| 715 | for (j = 1; j <= C->n; j++) |
---|
| 716 | xassert(work[j] == 0.0); |
---|
| 717 | xfree(work); |
---|
| 718 | return; |
---|
| 719 | } |
---|
| 720 | |
---|
| 721 | SPM *spm_add_mat(double alfa, const SPM *A, double beta, const SPM *B) |
---|
| 722 | { /* add two sparse matrices (driver routine) */ |
---|
| 723 | SPM *C; |
---|
| 724 | C = spm_add_sym(A, B); |
---|
| 725 | spm_add_num(C, alfa, A, beta, B); |
---|
| 726 | return C; |
---|
| 727 | } |
---|
| 728 | |
---|
| 729 | SPM *spm_mul_sym(const SPM *A, const SPM *B) |
---|
| 730 | { /* multiply two sparse matrices (symbolic phase) */ |
---|
| 731 | int i, j, k, *flag; |
---|
| 732 | SPM *C; |
---|
| 733 | xassert(A->n == B->m); |
---|
| 734 | /* create resultant matrix */ |
---|
| 735 | C = spm_create_mat(A->m, B->n); |
---|
| 736 | /* allocate and clear the flag array */ |
---|
| 737 | flag = xcalloc(1+C->n, sizeof(int)); |
---|
| 738 | for (j = 1; j <= C->n; j++) |
---|
| 739 | flag[j] = 0; |
---|
| 740 | /* compute pattern of C = A * B */ |
---|
| 741 | for (i = 1; i <= C->m; i++) |
---|
| 742 | { SPME *e, *ee; |
---|
| 743 | /* compute pattern of i-th row of C */ |
---|
| 744 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 745 | { k = e->j; |
---|
| 746 | for (ee = B->row[k]; ee != NULL; ee = ee->r_next) |
---|
| 747 | { j = ee->j; |
---|
| 748 | /* if a[i,k] != 0 and b[k,j] != 0 then c[i,j] != 0 */ |
---|
| 749 | if (!flag[j]) |
---|
| 750 | { /* c[i,j] does not exist, so create it */ |
---|
| 751 | spm_new_elem(C, i, j, 0.0); |
---|
| 752 | flag[j] = 1; |
---|
| 753 | } |
---|
| 754 | } |
---|
| 755 | } |
---|
| 756 | /* reset the flag array */ |
---|
| 757 | for (e = C->row[i]; e != NULL; e = e->r_next) |
---|
| 758 | flag[e->j] = 0; |
---|
| 759 | } |
---|
| 760 | /* check and deallocate the flag array */ |
---|
| 761 | for (j = 1; j <= C->n; j++) |
---|
| 762 | xassert(!flag[j]); |
---|
| 763 | xfree(flag); |
---|
| 764 | return C; |
---|
| 765 | } |
---|
| 766 | |
---|
| 767 | void spm_mul_num(SPM *C, const SPM *A, const SPM *B) |
---|
| 768 | { /* multiply two sparse matrices (numeric phase) */ |
---|
| 769 | int i, j; |
---|
| 770 | double *work; |
---|
| 771 | /* allocate and clear the working array */ |
---|
| 772 | work = xcalloc(1+A->n, sizeof(double)); |
---|
| 773 | for (j = 1; j <= A->n; j++) |
---|
| 774 | work[j] = 0.0; |
---|
| 775 | /* compute matrix C = A * B */ |
---|
| 776 | for (i = 1; i <= C->m; i++) |
---|
| 777 | { SPME *e, *ee; |
---|
| 778 | double temp; |
---|
| 779 | /* work := (i-th row of A) */ |
---|
| 780 | /* (note that A may have duplicate elements) */ |
---|
| 781 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 782 | work[e->j] += e->val; |
---|
| 783 | /* compute i-th row of C */ |
---|
| 784 | for (e = C->row[i]; e != NULL; e = e->r_next) |
---|
| 785 | { j = e->j; |
---|
| 786 | /* c[i,j] := work * (j-th column of B) */ |
---|
| 787 | temp = 0.0; |
---|
| 788 | for (ee = B->col[j]; ee != NULL; ee = ee->c_next) |
---|
| 789 | temp += work[ee->i] * ee->val; |
---|
| 790 | e->val = temp; |
---|
| 791 | } |
---|
| 792 | /* reset the working array */ |
---|
| 793 | for (e = A->row[i]; e != NULL; e = e->r_next) |
---|
| 794 | work[e->j] = 0.0; |
---|
| 795 | } |
---|
| 796 | /* check and deallocate the working array */ |
---|
| 797 | for (j = 1; j <= A->n; j++) |
---|
| 798 | xassert(work[j] == 0.0); |
---|
| 799 | xfree(work); |
---|
| 800 | return; |
---|
| 801 | } |
---|
| 802 | |
---|
| 803 | SPM *spm_mul_mat(const SPM *A, const SPM *B) |
---|
| 804 | { /* multiply two sparse matrices (driver routine) */ |
---|
| 805 | SPM *C; |
---|
| 806 | C = spm_mul_sym(A, B); |
---|
| 807 | spm_mul_num(C, A, B); |
---|
| 808 | return C; |
---|
| 809 | } |
---|
| 810 | |
---|
| 811 | PER *spm_create_per(int n) |
---|
| 812 | { /* create permutation matrix */ |
---|
| 813 | PER *P; |
---|
| 814 | int k; |
---|
| 815 | xassert(n >= 0); |
---|
| 816 | P = xmalloc(sizeof(PER)); |
---|
| 817 | P->n = n; |
---|
| 818 | P->row = xcalloc(1+n, sizeof(int)); |
---|
| 819 | P->col = xcalloc(1+n, sizeof(int)); |
---|
| 820 | /* initially it is identity matrix */ |
---|
| 821 | for (k = 1; k <= n; k++) |
---|
| 822 | P->row[k] = P->col[k] = k; |
---|
| 823 | return P; |
---|
| 824 | } |
---|
| 825 | |
---|
| 826 | void spm_check_per(PER *P) |
---|
| 827 | { /* check permutation matrix for correctness */ |
---|
| 828 | int i, j; |
---|
| 829 | xassert(P->n >= 0); |
---|
| 830 | for (i = 1; i <= P->n; i++) |
---|
| 831 | { j = P->row[i]; |
---|
| 832 | xassert(1 <= j && j <= P->n); |
---|
| 833 | xassert(P->col[j] == i); |
---|
| 834 | } |
---|
| 835 | return; |
---|
| 836 | } |
---|
| 837 | |
---|
| 838 | void spm_delete_per(PER *P) |
---|
| 839 | { /* delete permutation matrix */ |
---|
| 840 | xfree(P->row); |
---|
| 841 | xfree(P->col); |
---|
| 842 | xfree(P); |
---|
| 843 | return; |
---|
| 844 | } |
---|
| 845 | |
---|
| 846 | /* eof */ |
---|