mod2sparse.c

ldpc的11个程序 encode decode extract make-gen make-ldpc make-pchk print-gen print-pchk rand-src transm

      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;
}