COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpspm.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 24.3 KB
Line 
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 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
49SPM *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
92SPME *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
125void 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
162SPM *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
208SPM *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
263int 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
327SPM *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      }
371fini: 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
390int 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
420int 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
477SPM *spm_read_mat(const char *fname)
478{     xassert(fname != fname);
479      return NULL;
480}
481#else
482SPM *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);
529done: 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
575int spm_write_mat(const SPM *A, const char *fname)
576{     xassert(A != A);
577      xassert(fname != fname);
578      return 0;
579}
580#else
581int 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);
612done: 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
632SPM *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
644SPM *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
689void 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
721SPM *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
729SPM *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
767void 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
803SPM *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
811PER *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
826void 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
838void 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 */
Note: See TracBrowser for help on using the repository browser.