/* MOD2SPARSE.C - Procedures for handling sparse mod2 matrices. *//* Copyright (c) 2000, 2001 by Radford M. Neal  * * Permission is granted for anyone to copy, use, or modify this program  * for purposes of research or education, provided this copyright notice  * is retained, and note is made of any changes that have been made.  * * This program is distributed without any warranty, express or implied. * As this program was written for research purposes only, it has not been * tested to the degree that would be advisable in any important application. * All use of this program is entirely at the user's own risk. *//* NOTE:  See mod2sparse.html for documentation on these procedures. */#include <stdlib.h>#include <stdio.h>#include <math.h>#include "alloc.h"#include "intio.h"#include "mod2sparse.h"/* ALLOCATE AN ENTRY WITHIN A MATRIX.  This local procedure is used to   allocate a new entry, representing a non-zero element, within a given   matrix.  Entries in this matrix that were previously allocated and   then freed are re-used.  If there are no such entries, a new block   of entries is allocated. */static mod2entry *alloc_entry( mod2sparse *m){   mod2block *b;  mod2entry *e;  int k;  if (m->next_free==0)  {     b = chk_alloc (1, sizeof *b);    b->next = m->blocks;    m->blocks = b;    for (k = 0; k<Mod2sparse_block; k++)    { b->entry[k].left = m->next_free;      m->next_free = &b->entry[k];    }  }  e = m->next_free;  m->next_free = e->left;  e->pr = 0;  e->lr = 0;  return e;}/* ALLOCATE SPACE FOR A SPARSE MOD2 MATRIX.  */mod2sparse *mod2sparse_allocate( int n_rows, 		/* Number of rows in matrix */  int n_cols		/* Number of columns in matrix */){  mod2sparse *m;  mod2entry *e;  int i, j;  if (n_rows<=0 || n_cols<=0)  { fprintf(stderr,"mod2sparse_allocate: Invalid number of rows or columns\n");    exit(1);  }  m = chk_alloc (1, sizeof *m);  m->n_rows = n_rows;  m->n_cols = n_cols;  m->rows = chk_alloc (n_rows, sizeof *m->rows);  m->cols = chk_alloc (n_cols, sizeof *m->cols);  m->blocks = 0;  m->next_free = 0;  for (i = 0; i<n_rows; i++)  { e = &m->rows[i];    e->left = e->right = e->up = e->down = e;    e->row = e->col = -1;  }  for (j = 0; j<n_cols; j++)  { e = &m->cols[j];    e->left = e->right = e->up = e->down = e;    e->row = e->col = -1;  }  return m;}/* FREE SPACE OCCUPIED BY A SPARSE MOD2 MATRIX. */void mod2sparse_free( mod2sparse *m		/* Matrix to free */){   mod2block *b;  free(m->rows);  free(m->cols);  while (m->blocks!=0)  { b = m->blocks;    m->blocks = b->next;    free(b);  }}/* CLEAR A SPARSE MATRIX TO ALL ZEROS. */void mod2sparse_clear( mod2sparse *r){  mod2block *b;  mod2entry *e;  int i, j;  for (i = 0; i<mod2sparse_rows(r); i++)  { e = &r->rows[i];    e->left = e->right = e->up = e->down = e;  }  for (j = 0; j<mod2sparse_cols(r); j++)  { e = &r->cols[j];    e->left = e->right = e->up = e->down = e;  }  while (r->blocks!=0)  { b = r->blocks;    r->blocks = b->next;    free(b);  }}/* COPY A SPARSE MATRIX. */void mod2sparse_copy( mod2sparse *m,	/* Matrix to copy */  mod2sparse *r		/* Place to store copy of matrix */){  mod2entry *e, *f;  int i;  if (mod2sparse_rows(m)>mod2sparse_rows(r)    || mod2sparse_cols(m)>mod2sparse_cols(r))  { fprintf(stderr,"mod2sparse_copy: Destination matrix is too small\n");    exit(1);  }  mod2sparse_clear(r);  for (i = 0; i<mod2sparse_rows(m); i++)  {    e = mod2sparse_first_in_row(m,i);     while (!mod2sparse_at_end(e))    { f = mod2sparse_insert(r,e->row,e->col);      f->lr = e->lr;      f->pr = e->pr;      e = mod2sparse_next_in_row(e);    }  }}/* COPY ROWS OF A SPARSE MOD2 MATRIX. */void mod2sparse_copyrows( mod2sparse *m,	/* Matrix to copy */  mod2sparse *r,	/* Place to store copy of matrix */  int *rows		/* Indexes of rows to copy, from 0 */){   mod2entry *e;  int i;  if (mod2sparse_cols(m)>mod2sparse_cols(r))  { fprintf(stderr,     "mod2sparse_copyrows: Destination matrix has fewer columns than source\n");    exit(1);  }  mod2sparse_clear(r);  for (i = 0; i<mod2sparse_rows(r); i++)  { if (rows[i]<0 || rows[i]>=mod2sparse_rows(m))    { fprintf(stderr,"mod2sparse_copyrows: Row index out of range\n");      exit(1);    }    e = mod2sparse_first_in_row(m,rows[i]);    while (!mod2sparse_at_end(e))    { mod2sparse_insert(r,i,e->col);      e = mod2sparse_next_in_row(e);    }  }}/* COPY COLUMNS OF A SPARSE MOD2 MATRIX. */void mod2sparse_copycols( mod2sparse *m,	/* Matrix to copy */  mod2sparse *r,	/* Place to store copy of matrix */  int *cols		/* Indexes of columns to copy, from 0 */){   mod2entry *e;  int j;  if (mod2sparse_rows(m)>mod2sparse_rows(r))  { fprintf(stderr,      "mod2sparse_copycols: Destination matrix has fewer rows than source\n");    exit(1);  }  mod2sparse_clear(r);  for (j = 0; j<mod2sparse_cols(r); j++)  { if (cols[j]<0 || cols[j]>=mod2sparse_cols(m))    { fprintf(stderr,"mod2sparse_copycols: Column index out of range\n");      exit(1);    }    e = mod2sparse_first_in_col(m,cols[j]);    while (!mod2sparse_at_end(e))    { mod2sparse_insert(r,e->row,j);      e = mod2sparse_next_in_col(e);    }  }}/* PRINT A SPARSE MOD2 MATRIX IN HUMAN-READABLE FORM. */void mod2sparse_print( FILE *f,  mod2sparse *m){   int rdigits, cdigits;  mod2entry *e;  int i;  rdigits = mod2sparse_rows(m)<=10 ? 1           : mod2sparse_rows(m)<=100 ? 2          : mod2sparse_rows(m)<=1000 ? 3          : mod2sparse_rows(m)<=10000 ? 4          : mod2sparse_rows(m)<=100000 ? 5          : 6;  cdigits = mod2sparse_cols(m)<=10 ? 1           : mod2sparse_cols(m)<=100 ? 2          : mod2sparse_cols(m)<=1000 ? 3          : mod2sparse_cols(m)<=10000 ? 4          : mod2sparse_cols(m)<=100000 ? 5          : 6;  for (i = 0; i<mod2sparse_rows(m); i++)  {     fprintf(f,"%*d:",rdigits,i);    e = mod2sparse_first_in_row(m,i);    while (!mod2sparse_at_end(e))    { fprintf(f," %*d",cdigits,mod2sparse_col(e));      e = mod2sparse_next_in_row(e);    }    fprintf(f,"\n");  }}/* WRITE A SPARSE MOD2 MATRIX TO A FILE IN MACHINE-READABLE FORM. */int mod2sparse_write( FILE *f,  mod2sparse *m){  mod2entry *e;  int i;  intio_write(f,m->n_rows);  if (ferror(f)) return 0;  intio_write(f,m->n_cols);  if (ferror(f)) return 0;  for (i = 0; i<mod2sparse_rows(m); i++)  {     e = mod2sparse_first_in_row(m,i);    if (!mod2sparse_at_end(e))    {      intio_write (f, -(i+1));      if (ferror(f)) return 0;      while (!mod2sparse_at_end(e))      {         intio_write (f, mod2sparse_col(e)+1);        if (ferror(f)) return 0;        e = mod2sparse_next_in_row(e);      }    }  }  intio_write(f,0);  if (ferror(f)) return 0;  return 1;}/* READ A SPARSE MOD2 MATRIX STORED IN MACHINE-READABLE FORM FROM A FILE. */mod2sparse *mod2sparse_read( FILE *f){  int n_rows, n_cols;  mod2sparse *m;  int v, row, col;  n_rows = intio_read(f);  if (feof(f) || ferror(f) || n_rows<=0) return 0;  n_cols = intio_read(f);  if (feof(f) || ferror(f) || n_cols<=0) return 0;  m = mod2sparse_allocate(n_rows,n_cols);  row = -1;  for (;;)  {    v = intio_read(f);    if (feof(f) || ferror(f)) break;    if (v==0)    { return m;    }    else if (v<0)     { row = -v-1;      if (row>=n_rows) break;    }    else     { col = v-1;      if (col>=n_cols) break;      if (row==-1) break;      mod2sparse_insert(m,row,col);    }  }  /* Error if we get here. */  mod2sparse_free(m);  return 0;   }/* LOOK FOR AN ENTRY WITH GIVEN ROW AND COLUMN. */mod2entry *mod2sparse_find( mod2sparse *m,  int row,  int col){   mod2entry *re, *ce;  if (row<0 || row>=mod2sparse_rows(m) || col<0 || col>=mod2sparse_cols(m))  { fprintf(stderr,"mod2sparse_find: row or column index out of bounds\n");    exit(1);  }  /* Check last entries in row and column. */  re = mod2sparse_last_in_row(m,row);  if (mod2sparse_at_end(re) || mod2sparse_col(re)<col)   { return 0;  }  if (mod2sparse_col(re)==col)   { return re;  }  ce = mod2sparse_last_in_col(m,col);  if (mod2sparse_at_end(ce) || mod2sparse_row(ce)<row)  { return 0;  }  if (mod2sparse_row(ce)==row)  { return ce;  }  /* Search row and column in parallel, from the front. */  re = mod2sparse_first_in_row(m,row);  ce = mod2sparse_first_in_col(m,col);  for (;;)  {     if (mod2sparse_at_end(re) || mod2sparse_col(re)>col)    { return 0;    }     if (mod2sparse_col(re)==col)     { return re;    }    if (mod2sparse_at_end(ce) || mod2sparse_row(ce)>row)    { return 0;    }     if (mod2sparse_row(ce)==row)    { return ce;    }    re = mod2sparse_next_in_row(re);    ce = mod2sparse_next_in_col(ce);  }}/* INSERT AN ENTRY WITH GIVEN ROW AND COLUMN. */mod2entry *mod2sparse_insert( mod2sparse *m,  int row,  int col){  mod2entry *re, *ce, *ne;  if (row<0 || row>=mod2sparse_rows(m) || col<0 || col>=mod2sparse_cols(m))  { fprintf(stderr,"mod2sparse_insert: row or column index out of bounds\n");    exit(1);  }  /* Find old entry and return it, or allocate new entry and insert into row. */  re = mod2sparse_last_in_row(m,row);  if (!mod2sparse_at_end(re) && mod2sparse_col(re)==col)   { return re;  }  if (mod2sparse_at_end(re) || mod2sparse_col(re)<col)   { re = re->right;  }  else  {    re = mod2sparse_first_in_row(m,row);    for (;;)    {       if (!mod2sparse_at_end(re) && mod2sparse_col(re)==col)       { return re;      }      if (mod2sparse_at_end(re) || mod2sparse_col(re)>col)      { break;      }       re = mod2sparse_next_in_row(re);    }  }  ne = alloc_entry(m);  ne->row = row;  ne->col = col;  ne->left = re->left;  ne->right = re;  ne->left->right = ne;  ne->right->left = ne;  /* Insert new entry into column.  If we find an existing entry here,     the matrix must be garbled, since we didn't find it in the row. */  ce = mod2sparse_last_in_col(m,col);  if (!mod2sparse_at_end(ce) && mod2sparse_row(ce)==row)   { fprintf(stderr,"mod2sparse_insert: Garbled matrix\n");    exit(1);  }  if (mod2sparse_at_end(ce) || mod2sparse_row(ce)<row)   { ce = ce->down;  }  else  {    ce = mod2sparse_first_in_col(m,col);    for (;;)    {       if (!mod2sparse_at_end(ce) && mod2sparse_row(ce)==row)       { fprintf(stderr,"mod2sparse_insert: Garbled matrix\n");        exit(1);      }      if (mod2sparse_at_end(ce) || mod2sparse_row(ce)>row)      { break;      }       ce = mod2sparse_next_in_col(ce);    }  }      ne->up = ce->up;  ne->down = ce;  ne->up->down = ne;  ne->down->up = ne;  /* Return the new entry. */  return ne;}/* DELETE AN ENTRY FROM A SPARSE MATRIX. */void mod2sparse_delete( mod2sparse *m,  mod2entry *e){   if (e==0)  { fprintf(stderr,"mod2sparse_delete: Trying to delete a null entry\n");    exit(1);  }  if (e->row<0 || e->col<0)  { fprintf(stderr,"mod2sparse_delete: Trying to delete a header entry\n");    exit(1);  }  e->left->right = e->right;  e->right->left = e->left;  e->up->down = e->down;  e->down->up = e->up;   e->left = m->next_free;  m->next_free = e;}/* TEST WHETHER TWO SPARSE MATRICES ARE EQUAL. */int mod2sparse_equal( mod2sparse *m1,  mod2sparse *m2){  mod2entry *e1, *e2;  int i;  if (mod2sparse_rows(m1)!=mod2sparse_rows(m2)    || mod2sparse_cols(m1)!=mod2sparse_cols(m2))  { fprintf(stderr,"mod2sparse_equal: Matrices have different dimensions\n");    exit(1);  }    for (i = 0; i<mod2sparse_rows(m1); i++)  {     e1 = mod2sparse_first_in_row(m1,i);    e2 = mod2sparse_first_in_row(m2,i);    while (!mod2sparse_at_end(e1) && !mod2sparse_at_end(e2))    {        if (mod2sparse_col(e1)!=mod2sparse_col(e2))      { return 0;      }      e1 = mod2sparse_next_in_row(e1);      e2 = mod2sparse_next_in_row(e2);    }    if (!mod2sparse_at_end(e1) || !mod2sparse_at_end(e2))     { return 0;    }  }  return 1;}/* COMPUTE THE TRANSPOSE OF A SPARSE MOD2 MATRIX. */void mod2sparse_transpose( mod2sparse *m,	/* Matrix to compute transpose of (left unchanged) */  mod2sparse *r		/* Result of transpose operation */){  mod2entry *e;  int i;  if (mod2sparse_rows(m)!=mod2sparse_cols(r)    || mod2sparse_cols(m)!=mod2sparse_rows(r))  { fprintf(stderr,     "mod2sparse_transpose: Matrices have incompatible dimensions\n");    exit(1);  }  if (r==m)  { fprintf(stderr,      "mod2sparse_transpose: Result matrix is the same as the operand\n");    exit(1);  }  mod2sparse_clear(r);  for (i = 0; i<mod2sparse_rows(m); i++)  {    e = mod2sparse_first_in_row(m,i);    while (!mod2sparse_at_end(e))    { mod2sparse_insert(r,mod2sparse_col(e),i);      e = mod2sparse_next_in_row(e);    }  }}/* ADD TWO SPARSE MOD2 MATRICES. */void mod2sparse_add( mod2sparse *m1,	/* Left operand of add */  mod2sparse *m2,	/* Right operand of add */  mod2sparse *r		/* Place to store result of add */){  mod2entry *e1, *e2;  int i;  if (mod2sparse_rows(m1)!=mod2sparse_rows(r)    || mod2sparse_cols(m1)!=mod2sparse_cols(r)    || mod2sparse_rows(m2)!=mod2sparse_rows(r)   || mod2sparse_cols(m2)!=mod2sparse_cols(r))   { fprintf(stderr,"mod2sparse_add: Matrices have different dimensions\n");    exit(1);  }  if (r==m1 || r==m2)  { fprintf(stderr,     "mod2sparse_add: Result matrix is the same as one of the operands\n");    exit(1);  }  mod2sparse_clear(r);  for (i = 0; i<mod2sparse_rows(r); i++)  {     e1 = mod2sparse_first_in_row(m1,i);    e2 = mod2sparse_first_in_row(m2,i);    while (!mod2sparse_at_end(e1) && !mod2sparse_at_end(e2))    {       if (mod2sparse_col(e1)==mod2sparse_col(e2))      { e1 = mod2sparse_next_in_row(e1);        e2 = mod2sparse_next_in_row(e2);       }