lemon-project-template-glpk

comparison deps/glpk/src/glpspm.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:4c32059463f8
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 */