lemon-project-template-glpk

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