mod2sparse.c

用C语言编写的LDPC译码程序

      { 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.  Takes as input a matrix, A, 
   having M rows and N columns, and an integer K.  Finds an LU decomposition
   of a K by K sub-matrix of A.  The decomposition is stored in the matrix L,
   with M rows and K columns, and the matrix U, with K rows and N columns.
   The product of L and U will be equal to the K by K submatrix of A obtained
   by taking only rows and columns that are given in the first K elements of
   the 'rows' and 'cols' arrays, which are set by this procedure, with this
   sub-matrix distributed over the original M rows and N columns.  Furthermore,
   the ordering of the row and column indexes in these arrays will be set so 
   that if the rows of L and the columns of U were rearranged in this order, L 
   would be lower triangular, with zeros in rows past row K, and U would be 
   upper triangular, with zeros in columns past column K.  The 'rows' array
   is M long, and the 'cols' array is N long.  The elements in both arrays
   after the first K contain the indexes of the rows and columns not selected
   to be part of the sub-matrix of A, in arbitrary order.

   If A is not of rank K or more, L will contain fewer than K non-zero columns,
   and U will contain an equal number of non-zero rows.  The entries in the
   'rows' and 'cols' arrays for the extra zero rows or columns will be 
   arbitrary (but valid).  The number of extra zero columns is returned as 
   the value of this procedure (hence a return value of zero indicates that 
   a K by K sub-matrix of full rank was found).

   The matrix A is not altered.  The previous contents of L and U are cleared. 
*/

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, *e3;
  int i, j, k, cc, cc2, cc3, cr2, mxr, pr, fnd;
  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.  Solves Ly=x 
   for y by forward substitution, based on L being lower triangular after 
   its rows are reordered according to the given index array.  The vectors
   x and y are stored unpacked, one bit per character.  If L is M by K, 
   then x should be M long, but only the K bits indexed by 'rows' are looked
   at.  The solution vector, y, must be K long.  Only K rows of L are used,
   as also determined by the K indexes in the 'rows' argument.  If 'rows' is 
   null, the first K rows of L and the first K elements of x are used.  

   If the matrix L does not have 1s on its diagonal (after row rearrangement), 
   there may be no solution, depending on what x is.  If no solution exists, 
   this procedure returns zero, otherwise it returns one.  Any arbitrary
   bits in the solution are set to zero.

   If L is not lower-triangular, a message is displayed on standard error and 
   the program is terminated.
*/

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.  Solves Uz=y
   for z by backward substitution, based on U being upper triangular after 
   its colums are reordered according to the given index array.  The vectors
   y and z are stored unpacked, one bit per character.  If U is K by N,
   then the solution vector, z, should be N long, but only the K bits indexed 
   by 'cols' are set.  The vector y must be K long.  Only K columns of U are 
   used, as also determined by the K indexes in the 'cols' argument.  The other
   columns of U must be zero (this is not checked, but is necessary for the
   method used to work).  If 'cols' is null, the first K columns of U and the 
   first K elements of z are used.  

   If the matrix U does not have 1s on its diagonal (after column rearrangement)
   there may be no solution, depending on what y is.  If no solution exists, 
   this procedure returns zero, otherwise it returns one.  Any arbitrary
   bits in the solution are set to zero.

   If U is not upper-triangular, a message is displayed on standard error and 
   the program is terminated.
*/

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