root/mod2sparse.c
/* [<][>][^][v][top][bottom][index][help] */
/* 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;
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);
}
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;
}