1366 lines
28 KiB
C
Executable File
1366 lines
28 KiB
C
Executable File
/* MOD2SPARSE.C - Procedures for handling sparse mod2 matrices. */
|
|
|
|
/* Copyright (c) 1995-2012 by Radford M. Neal.
|
|
*
|
|
* Permission is granted for anyone to copy, use, modify, and distribute
|
|
* these programs and accompanying documents for any purpose, provided
|
|
* this copyright notice is retained and prominently displayed, and note
|
|
* is made of any changes made to these programs. These programs and
|
|
* documents are distributed without any warranty, express or implied.
|
|
* As the programs were written for research purposes only, they have not
|
|
* been tested to the degree that would be advisable in any important
|
|
* application. All use of these programs 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;
|
|
|
|
if (m)
|
|
{ free(m->rows);
|
|
free(m->cols);
|
|
|
|
while (m->blocks!=0)
|
|
{ b = m->blocks;
|
|
m->blocks = b->next;
|
|
free(b);
|
|
}
|
|
|
|
free (m);
|
|
}
|
|
}
|
|
|
|
|
|
/* 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);
|
|
}
|
|
|
|
else if (mod2sparse_col(e1)<mod2sparse_col(e2))
|
|
{ mod2sparse_insert(r,i,mod2sparse_col(e1));
|
|
e1 = mod2sparse_next_in_row(e1);
|
|
}
|
|
|
|
else
|
|
{ mod2sparse_insert(r,i,mod2sparse_col(e2));
|
|
e2 = mod2sparse_next_in_row(e2);
|
|
}
|
|
}
|
|
|
|
while (!mod2sparse_at_end(e1))
|
|
{ mod2sparse_insert(r,i,mod2sparse_col(e1));
|
|
e1 = mod2sparse_next_in_row(e1);
|
|
}
|
|
|
|
while (!mod2sparse_at_end(e2))
|
|
{ mod2sparse_insert(r,i,mod2sparse_col(e2));
|
|
e2 = mod2sparse_next_in_row(e2);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/* MULTIPLY TWO SPARSE MOD2 MATRICES. */
|
|
|
|
void mod2sparse_multiply
|
|
( mod2sparse *m1, /* Left operand of multiply */
|
|
mod2sparse *m2, /* Right operand of multiply */
|
|
mod2sparse *r /* Place to store result of multiply */
|
|
)
|
|
{
|
|
mod2entry *e1, *e2;
|
|
int i, j, b;
|
|
|
|
if (mod2sparse_cols(m1)!=mod2sparse_rows(m2)
|
|
|| mod2sparse_rows(m1)!=mod2sparse_rows(r)
|
|
|| mod2sparse_cols(m2)!=mod2sparse_cols(r))
|
|
{ fprintf (stderr,
|
|
"mod2sparse_multiply: Matrices have incompatible dimensions\n");
|
|
exit(1);
|
|
}
|
|
|
|
if (r==m1 || r==m2)
|
|
{ fprintf(stderr,
|
|
"mod2sparse_multiply: Result matrix is the same as one of the operands\n");
|
|
exit(1);
|
|
}
|
|
|
|
mod2sparse_clear(r);
|
|
|
|
for (i = 0; i<mod2sparse_rows(m1); i++)
|
|
{
|
|
if (mod2sparse_at_end(mod2sparse_first_in_row(m1,i)))
|
|
{ continue;
|
|
}
|
|
|
|
for (j = 0; j<mod2sparse_cols(m2); j++)
|
|
{
|
|
b = 0;
|
|
|
|
e1 = mod2sparse_first_in_row(m1,i);
|
|
e2 = mod2sparse_first_in_col(m2,j);
|
|
|
|
while (!mod2sparse_at_end(e1) && !mod2sparse_at_end(e2))
|
|
{
|
|
if (mod2sparse_col(e1)==mod2sparse_row(e2))
|
|
{ b ^= 1;
|
|
e1 = mod2sparse_next_in_row(e1);
|
|
e2 = mod2sparse_next_in_col(e2);
|
|
}
|
|
|
|
else if (mod2sparse_col(e1)<mod2sparse_row(e2))
|
|
{ e1 = mod2sparse_next_in_row(e1);
|
|
}
|
|
|
|
else
|
|
{ e2 = mod2sparse_next_in_col(e2);
|
|
}
|
|
}
|
|
|
|
if (b)
|
|
{ mod2sparse_insert(r,i,j);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/* MULTIPLY VECTOR BY SPARSE MATRIX. */
|
|
|
|
void mod2sparse_mulvec
|
|
( mod2sparse *m, /* The sparse matrix, with M rows and N columns */
|
|
char *u, /* The input vector, N long */
|
|
char *v /* Place to store the result, M long */
|
|
)
|
|
{
|
|
mod2entry *e;
|
|
int M, N;
|
|
int i, j;
|
|
|
|
M = mod2sparse_rows(m);
|
|
N = mod2sparse_cols(m);
|
|
|
|
for (i = 0; i<M; i++) v[i] = 0;
|
|
|
|
for (j = 0; j<N; j++)
|
|
{ if (u[j])
|
|
{ for (e = mod2sparse_first_in_col(m,j);
|
|
!mod2sparse_at_end(e);
|
|
e = mod2sparse_next_in_col(e))
|
|
{ v[mod2sparse_row(e)] ^= 1;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
/* COUNT ENTRIES IN A ROW. */
|
|
|
|
int mod2sparse_count_row
|
|
( mod2sparse *m,
|
|
int row
|
|
)
|
|
{
|
|
mod2entry *e;
|
|
int count;
|
|
|
|
if (row<0 || row>=mod2sparse_rows(m))
|
|
{ fprintf(stderr,"mod2sparse_count_row: row index out of bounds\n");
|
|
exit(1);
|
|
}
|
|
|
|
count = 0;
|
|
|
|
for (e = mod2sparse_first_in_row(m,row);
|
|
!mod2sparse_at_end(e);
|
|
e = mod2sparse_next_in_row(e))
|
|
{ count += 1;
|
|
}
|
|
|
|
return count;
|
|
}
|
|
|
|
|
|
/* COUNT ENTRIES IN A COLUMN. */
|
|
|
|
int mod2sparse_count_col
|
|
( mod2sparse *m,
|
|
int col
|
|
)
|
|
{
|
|
mod2entry *e;
|
|
int count;
|
|
|
|
if (col<0 || col>=mod2sparse_cols(m))
|
|
{ fprintf(stderr,"mod2sparse_count_col: column index out of bounds\n");
|
|
exit(1);
|
|
}
|
|
|
|
count = 0;
|
|
|
|
for (e = mod2sparse_first_in_col(m,col);
|
|
!mod2sparse_at_end(e);
|
|
e = mod2sparse_next_in_col(e))
|
|
{ count += 1;
|
|
}
|
|
|
|
return count;
|
|
}
|
|
|
|
|
|
/* ADD TO A ROW. */
|
|
|
|
void mod2sparse_add_row
|
|
( mod2sparse *m1, /* Matrix containing row to add to */
|
|
int row1, /* Index in this matrix of row to add to */
|
|
mod2sparse *m2, /* Matrix containing row to add from */
|
|
int row2 /* Index in this matrix of row to add from */
|
|
)
|
|
{
|
|
mod2entry *f1, *f2, *ft;
|
|
|
|
if (mod2sparse_cols(m1)<mod2sparse_cols(m2))
|
|
{ fprintf (stderr,
|
|
"mod2sparse_add_row: row added to is shorter than row added from\n");
|
|
exit(1);
|
|
}
|
|
|
|
if (row1<0 || row1>=mod2sparse_rows(m1)
|
|
|| row2<0 || row2>=mod2sparse_rows(m2))
|
|
{ fprintf (stderr,"mod2sparse_add_row: row index out of range\n");
|
|
exit(1);
|
|
}
|
|
|
|
f1 = mod2sparse_first_in_row(m1,row1);
|
|
f2 = mod2sparse_first_in_row(m2,row2);
|
|
|
|
while (!mod2sparse_at_end(f1) && !mod2sparse_at_end(f2))
|
|
{ if (mod2sparse_col(f1)>mod2sparse_col(f2))
|
|
{ mod2sparse_insert(m1,row1,mod2sparse_col(f2));
|
|
f2 = mod2sparse_next_in_row(f2);
|
|
}
|
|
else
|
|
{ ft = mod2sparse_next_in_row(f1);
|
|
if (mod2sparse_col(f1)==mod2sparse_col(f2))
|
|
{ mod2sparse_delete(m1,f1);
|
|
f2 = mod2sparse_next_in_row(f2);
|
|
}
|
|
f1 = ft;
|
|
}
|
|
}
|
|
|
|
while (!mod2sparse_at_end(f2))
|
|
{ mod2sparse_insert(m1,row1,mod2sparse_col(f2));
|
|
f2 = mod2sparse_next_in_row(f2);
|
|
}
|
|
}
|
|
|
|
|
|
/* ADD TO A COLUMN. */
|
|
|
|
void mod2sparse_add_col
|
|
( mod2sparse *m1, /* Matrix containing column to add to */
|
|
int col1, /* Index in this matrix of column to add to */
|
|
mod2sparse *m2, /* Matrix containing column to add from */
|
|
int col2 /* Index in this matrix of column to add from */
|
|
)
|
|
{
|
|
mod2entry *f1, *f2, *ft;
|
|
|
|
if (mod2sparse_rows(m1)<mod2sparse_rows(m2))
|
|
{ fprintf (stderr,
|
|
"mod2sparse_add_col: Column added to is shorter than column added from\n");
|
|
exit(1);
|
|
}
|
|
|
|
if (col1<0 || col1>=mod2sparse_cols(m1)
|
|
|| col2<0 || col2>=mod2sparse_cols(m2))
|
|
{ fprintf (stderr,"mod2sparse_add_col: Column index out of range\n");
|
|
exit(1);
|
|
}
|
|
|
|
f1 = mod2sparse_first_in_col(m1,col1);
|
|
f2 = mod2sparse_first_in_col(m2,col2);
|
|
|
|
while (!mod2sparse_at_end(f1) && !mod2sparse_at_end(f2))
|
|
{ if (mod2sparse_row(f1)>mod2sparse_row(f2))
|
|
{ mod2sparse_insert(m1,mod2sparse_row(f2),col1);
|
|
f2 = mod2sparse_next_in_col(f2);
|
|
}
|
|
else
|
|
{ ft = mod2sparse_next_in_col(f1);
|
|
if (mod2sparse_row(f1)==mod2sparse_row(f2))
|
|
{ mod2sparse_delete(m1,f1);
|
|
f2 = mod2sparse_next_in_col(f2);
|
|
}
|
|
f1 = ft;
|
|
}
|
|
}
|
|
|
|
while (!mod2sparse_at_end(f2))
|
|
{ mod2sparse_insert(m1,mod2sparse_row(f2),col1);
|
|
f2 = mod2sparse_next_in_col(f2);
|
|
}
|
|
}
|
|
|
|
|
|
/* FIND AN LU DECOMPOSITION OF A SPARSE MATRIX. */
|
|
|
|
int mod2sparse_decomp
|
|
( mod2sparse *A, /* Input matrix, M by N */
|
|
int K, /* Size of sub-matrix to find LU decomposition of */
|
|
mod2sparse *L, /* Matrix in which L is stored, M by K */
|
|
mod2sparse *U, /* Matrix in which U is stored, K by N */
|
|
int *rows, /* Array where row indexes are stored, M long */
|
|
int *cols, /* Array where column indexes are stored, N long */
|
|
mod2sparse_strategy strategy, /* Strategy to follow in picking rows/columns */
|
|
int abandon_number, /* Number of columns to abandon at some point */
|
|
int abandon_when /* When to abandon these columns */
|
|
)
|
|
{
|
|
int *rinv, *cinv, *acnt, *rcnt;
|
|
mod2sparse *B;
|
|
int M, N;
|
|
|
|
mod2entry *e, *f, *fn, *e2;
|
|
int i, j, k, cc, cc2, cc3, cr2, pr;
|
|
int found, nnf;
|
|
|
|
M = mod2sparse_rows(A);
|
|
N = mod2sparse_cols(A);
|
|
|
|
if (mod2sparse_cols(L)!=K || mod2sparse_rows(L)!=M
|
|
|| mod2sparse_cols(U)!=N || mod2sparse_rows(U)!=K)
|
|
{ fprintf (stderr,
|
|
"mod2sparse_decomp: Matrices have incompatible dimensions\n");
|
|
exit(1);
|
|
}
|
|
|
|
if (abandon_number>N-K)
|
|
{ fprintf(stderr,"Trying to abandon more columns than allowed\n");
|
|
exit(1);
|
|
}
|
|
|
|
rinv = chk_alloc (M, sizeof *rinv);
|
|
cinv = chk_alloc (N, sizeof *cinv);
|
|
|
|
if (abandon_number>0)
|
|
{ acnt = chk_alloc (M+1, sizeof *acnt);
|
|
}
|
|
|
|
if (strategy==Mod2sparse_minprod)
|
|
{ rcnt = chk_alloc (M, sizeof *rcnt);
|
|
}
|
|
|
|
mod2sparse_clear(L);
|
|
mod2sparse_clear(U);
|
|
|
|
/* Copy A to B. B will be modified, then discarded. */
|
|
|
|
B = mod2sparse_allocate(M,N);
|
|
mod2sparse_copy(A,B);
|
|
|
|
/* Count 1s in rows of B, if using minprod strategy. */
|
|
|
|
if (strategy==Mod2sparse_minprod)
|
|
{ for (i = 0; i<M; i++)
|
|
{ rcnt[i] = mod2sparse_count_row(B,i);
|
|
}
|
|
}
|
|
|
|
/* Set up initial row and column choices. */
|
|
|
|
for (i = 0; i<M; i++) rows[i] = rinv[i] = i;
|
|
for (j = 0; j<N; j++) cols[j] = cinv[j] = j;
|
|
|
|
/* Find L and U one column at a time. */
|
|
|
|
nnf = 0;
|
|
|
|
for (i = 0; i<K; i++)
|
|
{
|
|
/* Choose the next row and column of B. */
|
|
|
|
switch (strategy)
|
|
{
|
|
case Mod2sparse_first:
|
|
{
|
|
found = 0;
|
|
|
|
for (k = i; k<N; k++)
|
|
{ e = mod2sparse_first_in_col(B,cols[k]);
|
|
while (!mod2sparse_at_end(e))
|
|
{ if (rinv[mod2sparse_row(e)]>=i)
|
|
{ found = 1;
|
|
goto out_first;
|
|
}
|
|
e = mod2sparse_next_in_col(e);
|
|
}
|
|
}
|
|
|
|
out_first:
|
|
break;
|
|
}
|
|
|
|
case Mod2sparse_mincol:
|
|
{
|
|
found = 0;
|
|
|
|
for (j = i; j<N; j++)
|
|
{ cc2 = mod2sparse_count_col(B,cols[j]);
|
|
if (!found || cc2<cc)
|
|
{ e2 = mod2sparse_first_in_col(B,cols[j]);
|
|
while (!mod2sparse_at_end(e2))
|
|
{ if (rinv[mod2sparse_row(e2)]>=i)
|
|
{ found = 1;
|
|
cc = cc2;
|
|
e = e2;
|
|
k = j;
|
|
break;
|
|
}
|
|
e2 = mod2sparse_next_in_col(e2);
|
|
}
|
|
}
|
|
}
|
|
|
|
break;
|
|
}
|
|
|
|
case Mod2sparse_minprod:
|
|
{
|
|
found = 0;
|
|
|
|
for (j = i; j<N; j++)
|
|
{ cc2 = mod2sparse_count_col(B,cols[j]);
|
|
e2 = mod2sparse_first_in_col(B,cols[j]);
|
|
while (!mod2sparse_at_end(e2))
|
|
{ if (rinv[mod2sparse_row(e2)]>=i)
|
|
{ cr2 = rcnt[mod2sparse_row(e2)];
|
|
if (!found || cc2==1 || (cc2-1)*(cr2-1)<pr)
|
|
{ found = 1;
|
|
pr = cc2==1 ? 0 : (cc2-1)*(cr2-1);
|
|
e = e2;
|
|
k = j;
|
|
}
|
|
}
|
|
e2 = mod2sparse_next_in_col(e2);
|
|
}
|
|
}
|
|
|
|
break;
|
|
}
|
|
|
|
default:
|
|
{ fprintf(stderr,"mod2sparse_decomp: Unknown stategy\n");
|
|
exit(1);
|
|
}
|
|
}
|
|
|
|
if (!found)
|
|
{ nnf += 1;
|
|
}
|
|
|
|
/* Update 'rows' and 'cols'. Looks at 'k' and 'e' found above. */
|
|
|
|
if (found)
|
|
{
|
|
if (cinv[mod2sparse_col(e)]!=k) abort();
|
|
|
|
cols[k] = cols[i];
|
|
cols[i] = mod2sparse_col(e);
|
|
|
|
cinv[cols[k]] = k;
|
|
cinv[cols[i]] = i;
|
|
|
|
k = rinv[mod2sparse_row(e)];
|
|
|
|
if (k<i) abort();
|
|
|
|
rows[k] = rows[i];
|
|
rows[i] = mod2sparse_row(e);
|
|
|
|
rinv[rows[k]] = k;
|
|
rinv[rows[i]] = i;
|
|
}
|
|
|
|
/* Update L, U, and B. */
|
|
|
|
f = mod2sparse_first_in_col(B,cols[i]);
|
|
|
|
while (!mod2sparse_at_end(f))
|
|
{
|
|
fn = mod2sparse_next_in_col(f);
|
|
k = mod2sparse_row(f);
|
|
|
|
if (rinv[k]>i)
|
|
{ mod2sparse_add_row(B,k,B,mod2sparse_row(e));
|
|
if (strategy==Mod2sparse_minprod)
|
|
{ rcnt[k] = mod2sparse_count_row(B,k);
|
|
}
|
|
mod2sparse_insert(L,k,i);
|
|
}
|
|
else if (rinv[k]<i)
|
|
{ mod2sparse_insert(U,rinv[k],cols[i]);
|
|
}
|
|
else
|
|
{ mod2sparse_insert(L,k,i);
|
|
mod2sparse_insert(U,i,cols[i]);
|
|
}
|
|
|
|
f = fn;
|
|
}
|
|
|
|
/* Get rid of all entries in the current column of B, just to save space. */
|
|
|
|
for (;;)
|
|
{ f = mod2sparse_first_in_col(B,cols[i]);
|
|
if (mod2sparse_at_end(f)) break;
|
|
mod2sparse_delete(B,f);
|
|
}
|
|
|
|
/* Abandon columns of B with lots of entries if it's time for that. */
|
|
|
|
if (abandon_number>0 && i==abandon_when)
|
|
{
|
|
for (k = 0; k<M+1; k++)
|
|
{ acnt[k] = 0;
|
|
}
|
|
for (j = 0; j<N; j++)
|
|
{ k = mod2sparse_count_col(B,j);
|
|
acnt[k] += 1;
|
|
}
|
|
|
|
cc = abandon_number;
|
|
k = M;
|
|
while (acnt[k]<cc)
|
|
{ cc -= acnt[k];
|
|
k -= 1;
|
|
if (k<0) abort();
|
|
}
|
|
|
|
cc2 = 0;
|
|
for (j = 0; j<N; j++)
|
|
{ cc3 = mod2sparse_count_col(B,j);
|
|
if (cc3>k || (cc3==k && cc>0))
|
|
{ if (cc3==k) cc -= 1;
|
|
for (;;)
|
|
{ f = mod2sparse_first_in_col(B,j);
|
|
if (mod2sparse_at_end(f)) break;
|
|
mod2sparse_delete(B,f);
|
|
}
|
|
cc2 += 1;
|
|
}
|
|
}
|
|
|
|
if (cc2!=abandon_number) abort();
|
|
|
|
if (strategy==Mod2sparse_minprod)
|
|
{ for (j = 0; j<M; j++)
|
|
{ rcnt[j] = mod2sparse_count_row(B,j);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Get rid of all entries in the rows of L past row K, after reordering. */
|
|
|
|
for (i = K; i<M; i++)
|
|
{ for (;;)
|
|
{ f = mod2sparse_first_in_row(L,rows[i]);
|
|
if (mod2sparse_at_end(f)) break;
|
|
mod2sparse_delete(L,f);
|
|
}
|
|
}
|
|
|
|
mod2sparse_free(B);
|
|
free(rinv);
|
|
free(cinv);
|
|
if (strategy==Mod2sparse_minprod) free(rcnt);
|
|
if (abandon_number>0) free(acnt);
|
|
|
|
return nnf;
|
|
}
|
|
|
|
|
|
/* SOLVE A LOWER-TRIANGULAR SYSTEM BY FORWARD SUBSTITUTION. */
|
|
|
|
int mod2sparse_forward_sub
|
|
( mod2sparse *L, /* Matrix that is lower triangular after reordering */
|
|
int *rows, /* Array of indexes (from 0) of rows for new order */
|
|
char *x, /* Vector on right of equation, also reordered */
|
|
char *y /* Place to store solution */
|
|
)
|
|
{
|
|
int K, i, j, ii, b, d;
|
|
mod2entry *e;
|
|
|
|
K = mod2sparse_cols(L);
|
|
|
|
/* Make sure that L is lower-triangular, after row re-ordering. */
|
|
|
|
for (i = 0; i<K; i++)
|
|
{ ii = rows ? rows[i] : i;
|
|
e = mod2sparse_last_in_row(L,ii);
|
|
if (!mod2sparse_at_end(e) && mod2sparse_col(e)>i)
|
|
{ fprintf(stderr,
|
|
"mod2sparse_forward_sub: Matrix is not lower-triangular\n");
|
|
exit(1);
|
|
}
|
|
}
|
|
|
|
/* Solve system by forward substitution. */
|
|
|
|
for (i = 0; i<K; i++)
|
|
{
|
|
ii = rows ? rows[i] : i;
|
|
|
|
/* Look at bits in this row, forming inner product with partial
|
|
solution, and seeing if the diagonal is 1. */
|
|
|
|
d = 0;
|
|
b = 0;
|
|
|
|
for (e = mod2sparse_first_in_row(L,ii);
|
|
!mod2sparse_at_end(e);
|
|
e = mod2sparse_next_in_row(e))
|
|
{
|
|
j = mod2sparse_col(e);
|
|
|
|
if (j==i)
|
|
{ d = 1;
|
|
}
|
|
else
|
|
{ b ^= y[j];
|
|
}
|
|
}
|
|
|
|
/* Check for no solution if the diagonal isn't 1. */
|
|
|
|
if (!d && b!=x[ii])
|
|
{ return 0;
|
|
}
|
|
|
|
/* Set bit of solution, zero if arbitrary. */
|
|
|
|
y[i] = b^x[ii];
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
|
|
|
|
/* SOLVE AN UPPER-TRIANGULAR SYSTEM BY BACKWARD SUBSTITUTION. */
|
|
|
|
int mod2sparse_backward_sub
|
|
( mod2sparse *U, /* Matrix that is upper triangular after reordering */
|
|
int *cols, /* Array of indexes (from 0) of columns for new order */
|
|
char *y, /* Vector on right of equation */
|
|
char *z /* Place to store solution, also reordered */
|
|
)
|
|
{
|
|
int K, i, j, ii, b, d;
|
|
mod2entry *e;
|
|
|
|
K = mod2sparse_rows(U);
|
|
|
|
/* Make sure that U is upper-triangular, after column re-ordering. */
|
|
|
|
for (i = 0; i<K; i++)
|
|
{ ii = cols ? cols[i] : i;
|
|
e = mod2sparse_last_in_col(U,ii);
|
|
if (!mod2sparse_at_end(e) && mod2sparse_row(e)>i)
|
|
{ fprintf(stderr,
|
|
"mod2sparse_backward_sub: Matrix is not upper-triangular\n");
|
|
exit(1);
|
|
}
|
|
}
|
|
|
|
/* Solve system by backward substitution. */
|
|
|
|
for (i = K-1; i>=0; i--)
|
|
{
|
|
ii = cols ? cols[i] : i;
|
|
|
|
/* Look at bits in this row, forming inner product with partial
|
|
solution, and seeing if the diagonal is 1. */
|
|
|
|
d = 0;
|
|
b = 0;
|
|
|
|
for (e = mod2sparse_first_in_row(U,i);
|
|
!mod2sparse_at_end(e);
|
|
e = mod2sparse_next_in_row(e))
|
|
{
|
|
j = mod2sparse_col(e);
|
|
|
|
if (j==ii)
|
|
{ d = 1;
|
|
}
|
|
else
|
|
{ b ^= z[j];
|
|
}
|
|
}
|
|
|
|
/* Check for no solution if the diagonal isn't 1. */
|
|
|
|
if (!d && b!=y[i])
|
|
{ return 0;
|
|
}
|
|
|
|
/* Set bit of solution, zero if arbitrary. */
|
|
|
|
z[ii] = b^y[i];
|
|
}
|
|
|
|
return 1;
|
|
}
|