mod2dense.c

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

}


/* MULTIPLY TWO DENSE MOD2 MATRICES.  Multiplies the first matrix by the second 
   matrix, storing the result in the third matrix.  Neither of the first 
   two matrices is changed by this operation.  The three matrices must have 
   compatible numbers of rows and columns. 

   The algorithm used runs faster if the second matrix (right operand of the
   multiply) is sparse, but it is also appropriate for dense matrices.  This
   procedure could be speeded up a bit by replacing the call of mod2dense_get
   with in-line code that avoids division, but this doesn't seem worthwhile
   at the moment. 

   The result matrix must not be the same as either operand. */

void mod2dense_multiply 
( mod2dense *m1, 	/* Left operand of multiply */
  mod2dense *m2,	/* Right operand of multiply */
  mod2dense *r		/* Place to store result of multiply */
)
{
  int i, j, k;

  if (mod2dense_cols(m1)!=mod2dense_rows(m2) 
   || mod2dense_rows(m1)!=mod2dense_rows(r) 
   || mod2dense_cols(m2)!=mod2dense_cols(r))
  { fprintf(stderr,
     "mod2dense_multiply: Matrices have incompatible dimensions\n");
    exit(1);
  }

  if (r==m1 || r==m2)
  { fprintf(stderr,
      "mod2dense_multiply: Result matrix is the same as one of the operands\n");
    exit(1);
  }

  mod2dense_clear(r);

  for (j = 0; j<mod2dense_cols(r); j++)
  { for (i = 0; i<mod2dense_rows(m2); i++)
    { if (mod2dense_get(m2,i,j))
      { for (k = 0; k<r->n_words; k++)
        { r->col[j][k] ^= m1->col[i][k];
        }
      }
    }
  }
}


/* SEE WHETHER TWO DENSE MOD2 MATRICES ARE EQUAL.  Returns one if every 
   element of the first matrix is equal to the corresponding element of the
   second matrix.  The two matrices must have the same number of rows and the
   same number of columns. */

int mod2dense_equal
( mod2dense *m1,
  mod2dense *m2
)
{
  int k, j, w;
  mod2word m;

  if (mod2dense_rows(m1)!=mod2dense_rows(m2) 
   || mod2dense_cols(m1)!=mod2dense_cols(m2))
  { fprintf(stderr,"mod2dense_equal: Matrices have different dimensions\n");
    exit(1);
  }

  w = m1->n_words;

  /* Form a mask that has 1s in the lower bit positions corresponding to
     bits that contain information in the last word of a matrix column. */

  m = (1 << (mod2_wordsize - (w*mod2_wordsize-m1->n_rows))) - 1;
  
  for (j = 0; j<mod2dense_cols(m1); j++)
  {
    for (k = 0; k<w-1; k++)
    { if (m1->col[j][k] != m2->col[j][k]) return 0;
    }

    if ((m1->col[j][k]&m) != (m2->col[j][k]&m)) return 0;
  }

  return 1;
}


/* INVERT A DENSE MOD2 MATRIX.  Inverts the first matrix passed, destroying 
   its contents in the process (though it remains a valid matrix for storing
   into later).  The matrix to be inverted must have the same number of
   rows as columns.  The inverse of this matrix is stored in the matrix 
   passed as the second argument (which must already exist, and have the 
   same number of rows and columns as the first).

   The value returned by mod2dense_invert is one if the inversion was 
   successful, zero if the matrix turned out to be singular (in which 
   case the contents of both the original matrix and the result matrix 
   will be garbage). 

   The result matrix must not be the same as the operand matrix. 

   The algorithm used is based on inverting A by transforming the equation
   AI = A to the equation AB = I using column operations, at which point B 
   is the inverse of A.  The representation of matrices used allows easy
   swapping of columns as needed by fiddling pointers.
*/

int mod2dense_invert 
( mod2dense *m,		/* The matrix to find the inverse of (destroyed) */
  mod2dense *r		/* Place to store the inverse */
)
{
  mod2word *s, *t;
  int i, j, k, n, w, k0, b0;

  if (mod2dense_rows(m)!=mod2dense_cols(m))
  { fprintf(stderr,"mod2dense_invert: Matrix to invert is not square\n");
    exit(1);
  }

  if (r==m)
  { fprintf(stderr, 
      "mod2dense_invert: Result matrix is the same as the operand\n");
    exit(1);
  }

  n = mod2dense_rows(m);
  w = m->n_words;

  if (mod2dense_rows(r)!=n || mod2dense_cols(r)!=n)
  { fprintf(stderr,
     "mod2dense_invert: Matrix to receive inverse has wrong dimensions\n");
    exit(1);
  }

  mod2dense_clear(r);
  for (i = 0; i<n; i++) 
  { mod2dense_set(r,i,i,1);
  }

  for (i = 0; i<n; i++)
  { 
    k0 = i >> mod2_wordsize_shift;
    b0 = i & mod2_wordsize_mask;

    for (j = i; j<n; j++) 
    { if (mod2_getbit(m->col[j][k0],b0)) break;
    }

    if (j==n) return 0;

    if (j!=i)
    {
      t = m->col[i];
      m->col[i] = m->col[j];
      m->col[j] = t;

      t = r->col[i];
      r->col[i] = r->col[j];
      r->col[j] = t;
    }

    for (j = 0; j<n; j++)
    { if (j!=i && mod2_getbit(m->col[j][k0],b0))
      { s = m->col[j];
        t = m->col[i];
        for (k = k0; k<w; k++) s[k] ^= t[k];
        s = r->col[j];
        t = r->col[i];
        for (k = 0; k<w; k++) s[k] ^= t[k];
      }
    }
  }

  return 1;
}


/* INVERT A DENSE MOD2 MATRIX WITH COLUMNS SELECTED FROM A BIGGER MATRIX. 
   Inverts the matrix obtained by selecting certain columns from the first
   matrix passed (which must have at least as many columns as rows), storing
   the inverse in the second matrix passed.  The second matrix must already 
   exist, and must have the same number of rows and columns as the first 
   matrix.  The result is stored in the selected columns of the result matrix,  
   with the other columns being set to garbage. (Normally, one would extract 
   just the relevant columns afterwards using mod2dense_copycols.)  The 
   contents of the first matrix are destroyed (though it remains a valid 
   matrix for storing into later).

   Indexes of the columns selected are stored, in order, in cols, followed 
   by the columns not selected (in arbitrary order).  The value returned is 
   zero if an invertible matrix was found, and otherwise the number of
   columns left to find when inversion failed.  If inversion fails, the
   partial results are stored, with the selection of columns after that
   point being arbitrary (but valid).

   Note that when the first matrix is square, and non-singular, the result 
   is NOT in general the same as that obtained by calling mod2dense_invert, 
   which orders the columns of the inverse so that it applies to the original 
   ordering of the columns of the first matrix.

   The result matrix must not be the same as the operand matrix. 

   See mod2dense_invert for remarks on the algorithm.
*/

int mod2dense_invert_selected 
( mod2dense *m,		/* Matrix from which to pick a submatrix to invert */
  mod2dense *r,		/* Place to store the inverse */
  int *cols		/* Set to indexes of columns used and not used */
)
{
  mod2word *s, *t;
  int i, j, k, n, n2, w, k0, b0, c, R;

  if (mod2dense_rows(m)>mod2dense_cols(m))
  { fprintf(stderr,"mod2dense_invert_selected: Matrix has too few columns\n");
    exit(1);
  }

  if (r==m)
  { fprintf(stderr, 
      "mod2dense_invert_selected: Result matrix is the same as the operand\n");
    exit(1);
  }

  n = mod2dense_rows(m);
  w = m->n_words;

  n2 = mod2dense_cols(m);

  if (mod2dense_rows(r)!=n || mod2dense_cols(r)!=n2)
  { fprintf(stderr,
 "mod2dense_invert_selected: Matrix to receive inverse has wrong dimensions\n");
    exit(1);
  }

  mod2dense_clear(r);

  for (j = 0; j<n2; j++)
  { cols[j] = j;
  }

  R = 0;

  for (i = 0; i<n; i++)
  { 
    k0 = i >> mod2_wordsize_shift;
    b0 = i & mod2_wordsize_mask;

    for (j = i; j<n2; j++) 
    { if (mod2_getbit(m->col[cols[j]][k0],b0)) break;
    }

    if (j==n2) 
    { R += 1;
      continue;
    }

    c = cols[j];
    cols[j] = cols[i];
    cols[i] = c;

    mod2dense_set(r,i,c,1);

    for (j = 0; j<n2; j++)
    { if (j!=c && mod2_getbit(m->col[j][k0],b0))
      { s = m->col[j];
        t = m->col[c];
        for (k = k0; k<w; k++) s[k] ^= t[k];
        s = r->col[j];
        t = r->col[c];
        for (k = 0; k<w; k++) s[k] ^= t[k];
      }
    }
  }

  return R;
}


/* FORCIBLY INVERT A DENSE MOD2 MATRIX.  Inverts the first matrix passed, 
   destroying its contents in the process (though it remains a valid matrix 
   for storing into later).  The matrix to be inverted must have the same 
   number of rows as columns.  The inverse of this matrix is stored in the 
   matrix passed as the second argument (which must already exist, and have 
   the same number of rows and columns as the first).

   If the matrix to be inverted is singular, the inversion proceeds anyway,
   with bits of the matrix being changed as needed to create a non-singular
   matrix.  The value returned by mod2dense_invert is the number of elements 
   that had to be changed to make inversion possible (zero, if the original 
   matrix was non-singular). 

   The row and column indexes of the elements of the original matrix 
   that were changed are stored in the arrays passed as the last two
   elements.  These arrays must have as many elements as the dimension
   of the matrix.  (This is so even if it is known that fewer elements
   than this will be changed, as these arrays are also used as temporary
   storage by this routine.) 

   The result matrix must not be the same as the operand matrix. */

int mod2dense_forcibly_invert 
( mod2dense *m, 	/* The matrix to find the inverse of (destroyed) */
  mod2dense *r,		/* Place to store the inverse */
  int *a_row,		/* Place to store row indexes of altered elements */
  int *a_col		/* Place to store column indexes of altered elements */
)
{
  mod2word *s, *t;
  int i, j, k, n, w, k0, b0;
  int u, c;

  if (mod2dense_rows(m)!=mod2dense_cols(m))
  { fprintf(stderr,
      "mod2dense_forcibly_invert: Matrix to invert is not square\n");
    exit(1);
  }

  if (r==m)
  { fprintf(stderr, 
      "mod2dense_forcibly_invert: Result matrix is the same as the operand\n");
    exit(1);
  }

  n = mod2dense_rows(m);
  w = m->n_words;

  if (mod2dense_rows(r)!=n || mod2dense_cols(r)!=n)
  { fprintf(stderr,
 "mod2dense_forcibly_invert: Matrix to receive inverse has wrong dimensions\n");
    exit(1);
  }

  mod2dense_clear(r);
  for (i = 0; i<n; i++) 
  { mod2dense_set(r,i,i,1);
  }

  for (i = 0; i<n; i++)
  { a_row[i] = -1;
    a_col[i] = i;
  }

  for (i = 0; i<n; i++)
  { 
    k0 = i >> mod2_wordsize_shift;
    b0 = i & mod2_wordsize_mask;

    for (j = i; j<n; j++) 
    { if (mod2_getbit(m->col[j][k0],b0)) break;
    }

    if (j==n)
    { j = i;
      mod2dense_set(m,i,j,1);
      a_row[i] = i;
    }

    if (j!=i)
    { 
      t = m->col[i];
      m->col[i] = m->col[j];
      m->col[j] = t;

      t = r->col[i];
      r->col[i] = r->col[j];
      r->col[j] = t;

      u = a_col[i];
      a_col[i] = a_col[j];
      a_col[j] = u;
    }

    for (j = 0; j<n; j++)
    { if (j!=i && mod2_getbit(m->col[j][k0],b0))
      { s = m->col[j];
        t = m->col[i];
        for (k = k0; k<w; k++) s[k] ^= t[k];
        s = r->col[j];
        t = r->col[i];
        for (k = 0; k<w; k++) s[k] ^= t[k];
      }
    }
  }

  c = 0;
  for (i = 0; i<n; i++)
  { if (a_row[i]!=-1)
    { a_row[c] = a_row[i];
      a_col[c] = a_col[i];
      c += 1;
    }
  }

  return c;
}