www.pudn.com > 7894509matlab.rar > ldpc_h2g.c, change:2000-05-09,size:12505b

/* Invert sparse binary H for  LDPC*/
/* Author : Igor Kozintsev   igor@ifp.uiuc.edu
Please let me know if you find bugs in this code (I did test
it but I still have some doubts). All other comments are welcome
too :) !
I use a simple algorithm to invert H.
We convert H to [I | A]
[junk ]
using column reodering and row operations (junk - a few rows of H
which are linearly dependent on the previous ones)
G is then found as G = [A'|I]
G is stored as array of doubles in Matlab which is very inefficient.
Internal representation in this programm is unsigned char. Please modify
the part which writes G if you wish.
*/
#include <math.h>
#include "mex.h"

/* Input Arguments: tentative H matrix*/
#define	H_IN	prhs[0]
#define Q_IN    prhs[1] /* field base */

/* Output Arguments: final matrices*/
#define	H_OUT	plhs[0]
#define	G_OUT	plhs[1]

/************************************ GFq math *******************************/

/* This file contains lookup tables and routines required
* to perform the field operations over GF(q), i.e. addition
* and multiplication.
*
* Addition is easy, we use exclusive-or operation.
* For multiplication we need two tables for each q.
* The first is the logarithm table, and the second
* is the exponential table.  We catch multiplication
* by zero and one separately.
*
* Source for tables: MacWilliams and Sloane, fig 4.5
*
* WARNING: for speed, no check is made that legal values
*          are supplied.
*/
const int log4[4]   = {0,0,1,2}; /* stores i at address \alpha^i */
const int log8[8]   = {0,0,1,3,2,6,4,5};
const int log16[16] = {0,0,1,4,2,8,5,10,3,14,9,7,6,13,11,12};
const int log32[32] = {0,0,1,18,2,5,19,11,3,29,6,27,20,8,12,23,4,\
10,30,17,7,22,28,26,21,25,9,16,13,14,24,15};
const int log64[64] = {0,0,1,6,2,12,7,26,3,32,13,35,8,48,27,18,4,24,\
33,16,14,52,36,54,9,45,49,38,28,41,19,56,5,62,\
25,11,34,31,17,47,15,23,53,51,37,44,55,40,10,\
61,46,30,50,22,39,43,29,60,42,21,20,59,57,58};
const int log128[128] = {0,0,1,31,2,62,32,103,3,7,63,15,33,84,104,\
93, 4,124,8,121,64,79,16,115,34,11,85,38,105,46,94,51,\
5,82,125,60,9,44,122,77,65,67,80,42,17,69,116,23,35,118,\
12,28,86,25,39,57,106,19,47,89,95,71,52,110,6,14,83,92,126,\
30,61,102,10,37,45,50,123,120,78,114,66,41,68,22,81,59,43,76,\
18,88,70,109,117,27,24,56,36,49,119,113,13,91,29,101,87,108,\
26,55,40,21,58,75,107,54,20,74,48,112,90,100,96,97,72,98,53,73,111,99};
const int log256[256] = {0,0,1,25,2,50,26,198,3,223,51,238,27,104,199,75,4,100,\
224,14,52,141,239,129,28,193,105,248,200,8,76,113,5,138,101,47,225,\
36,15,33,53,147,142,218,240,18,130,69,29,181,194,125,106,39,249,185,\
201,154,9,120,77,228,114,166,6,191,139,98,102,221,48,253,226,152,37,\
179,16,145,34,136,54,208,148,206,143,150,219,189,241,210,19,92,131,\
56,70,64,30,66,182,163,195,72,126,110,107,58,40,84,250,133,186,61,202,\
94,155,159,10,21,121,43,78,212,229,172,115,243,167,87,7,112,192,247,\
140,128,99,13,103,74,222,237,49,197,254,24,227,165,153,119,38,184,180,\
124,17,68,146,217,35,32,137,46,55,63,209,91,149,188,207,205,144,135,151,\
178,220,252,190,97,242,86,211,171,20,42,93,158,132,60,57,83,71,109,65,\
162,31,45,67,216,183,123,164,118,196,23,73,236,127,12,111,246,108,161,59,\
82,41,157,85,170,251,96,134,177,187,204,62,90,203,89,95,176,156,169,160,\
81,11,245,22,235,122,117,44,215,79,174,213,233,230,231,173,232,116,214,\
244,234,168,80,88,175};

const int exp4[3]   = {1,2,3}; /* stores \alpha^i at address i */
const int exp8[7]   = {1,2,4,3,6,7,5};
const int exp16[15] = {1,2,4,8,3,6,12,11,5,10,7,14,15,13,9};
const int exp32[31] = {1,2,4,8,16,5,10,20,13,26,17,7,14,28,29,31,\
27,19,3,6,12,24,21,15,30,25,23,11,22,9,18};
const int exp64[63] = {1,2,4,8,16,32,3,6,12,24,48,35,5,10,20,40,19,\
38,15,30,60,59,53,41,17,34,7,14,28,56,51,37,\
9,18,36,11,22,44,27,54,47,29,58,55,45,25,50,\
39,13,26,52,43,21,42,23,46,31,62,63,61,57,49,33};
const int exp128[127] = {1,2,4,8,16,32,64,9,18,36,72,25,50,100,65,11,\
22,44,88,57,114,109,83,47,94,53,106,93,51,102,69,3,6,12,24,\
48,96,73,27,54,108,81,43,86,37,74,29,58,116,97,75,31,62,124,\
113,107,95,55,110,85,35,70,5,10,20,40,80,41,82,45,90,61,122,\
125,115,111,87,39,78,21,42,84,33,66,13,26,52,104,89,59,118,101,\
67,15,30,60,120,121,123,127,119,103,71,7,14,28,56,112,105,91,63,\
126,117,99,79,23,46,92,49,98,77,19,38,76,17,34,68};

const int exp256[255] = {1,2,4,8,16,32,64,128,29,58,116,232,205,135,19,38,76,\
152,45,90,180,117,234,201,143,3,6,12,24,48,96,192,157,39,78,156,\
37,74,148,53,106,212,181,119,238,193,159,35,70,140,5,10,20,40,80,\
160,93,186,105,210,185,111,222,161,95,190,97,194,153,47,94,188,101,\
202,137,15,30,60,120,240,253,231,211,187,107,214,177,127,254,\
225,223,163,91,182,113,226,217,175,67,134,17,34,68,136,13,26,52,104,\
208,189,103,206,129,31,62,124,248,237,199,147,59,118,236,197,151,51,\
102,204,133,23,46,92,184,109,218,169,79,158,33,66,132,21,42,84,168,\
77,154,41,82,164,85,170,73,146,57,114,228,213,183,115,230,209,191,99,\
198,145,63,126,252,229,215,179,123,246,241,255,227,219,171,75,150,49,\
98,196,149,55,110,220,165,87,174,65,130,25,50,100,200,141,7,14,28,56,\
112,224,221,167,83,166,81,162,89,178,121,242,249,239,195,155,43,86,172,\
69,138,9,18,36,72,144,61,122,244,245,247,243,251,235,203,139,11,22,44,\
88,176,125,250,233,207,131,27,54,108,216,173,71,142};

/* For testing:
main(){

u_int i,j,q;

while(1){
printf("please enter a b q:");
scanf("%u %u %u",&i,&j,&q);
printf("\n a * b in GF(q) = %u\n",GFq_m(i,j,q));
}
}
*/

int GFq_m(int a, int b, int q)
{

if ( a == 0 || b == 0 )  return 0 ;
if ( a == 1 )  return b ;
if ( b == 1 )  return a ;
switch (q){
case 256:
return exp256[(log256[a]+log256[b])%255];
case 128:
return exp128[(log128[a]+log128[b])%127];
case 64:
return exp64[(log64[a]+log64[b])%63];
case 32:
return exp32[(log32[a]+log32[b])%31];
case 16:
return exp16[(log16[a]+log16[b])%15];
case 8:
return exp8[(log8[a]+log8[b])%7];
case 4:
return exp4[(log4[a]+log4[b])%3];
}
mexErrMsgTxt(1,"GFq_m: I'm afraid I don't know how to multiply in GFq\n");
return 0 ;
}

int GFq_inv(int a, int q)
{

if ( a == 0) mexErrMsgTxt(1,"GFq_inv: no inverse for 0!\n");;
if ( a == 1 )  return 1 ;
switch (q){
case 256:
return exp256[(255-log256[a])];
case 128:
return exp128[(127-log128[a])];
case 64:
return exp64[(63-log64[a])];
case 32:
return exp32[(31-log32[a])];
case 16:
return exp16[(15-log16[a])];
case 8:
return exp8[(7-log8[a])];
case 4:
return exp4[(3-log4[a])];
}
mexErrMsgTxt(1,"GFq_inv: not defined inverse for  GFq\n");
return 0 ;
}

int GFq_a(int a, int b)
{
return a^b;
}

/************************************ end GFq math *******************************/

void mexFunction(
int nlhs,       mxArray *plhs[],
int nrhs, const mxArray *prhs[]
)
{
unsigned char **HH, **GG;
int ii, jj, *ir, *jc, rdep, tmp, d, q, scale;
double *sr1, *sr2, *g;
int N,M,K,i,j,k,kk,nz,*irs1,*jcs1, *irs2, *jcs2;

/* Check for proper number of arguments */
if (nrhs != 2) {
mexErrMsgTxt("h2g requires two input arguments.");
} else if (nlhs != 2) {
mexErrMsgTxt("h2g requires two output arguments.");
} else if (!mxIsSparse(H_IN)) {
mexErrMsgTxt("h2g requires sparse H matrix.");
}

/* get the field base */
q = (int)mxGetScalar(Q_IN);

/* read sparse matrix H */
sr1  = mxGetPr(H_IN);
irs1 = mxGetIr(H_IN);  /* row */
jcs1 = mxGetJc(H_IN);  /* column */
nz = mxGetNzmax(H_IN); /* number of nonzero elements (they are ones)*/
M = mxGetM(H_IN);
N = mxGetN(H_IN);

/* create working array HH[row][column]*/
HH = (unsigned char **)mxMalloc(M*sizeof(unsigned char *));
for(i=0 ; i<M ; i++){
HH[i] = (unsigned char *)mxMalloc(N*sizeof(unsigned char));
}
for(i=0 ; i<M ; i++)
for(j=0 ; j<N ; j++)
HH[i][j] = 0; /* initialize all to zero */

k=0;
for(j=0 ; j<N ; j++) {
for(i=0 ; i<(jcs1[j+1]-jcs1[j]) ; i++) {
ii = irs1[k]; /* index in column j*/
HH[ii][j] = (unsigned char)sr1[k]; /* put  nonzeros */
k++;
}
}

/* invert HH matrix here */
/* row and column indices */
ir = (int *)mxMalloc(M*sizeof(int));
jc = (int *)mxMalloc(N*sizeof(int));
for( i=0 ; i<M ; i++)
ir[i] = i;
for( j=0 ; j<N ; j++)
jc[j] = j;

/* perform Gaussian elimination on H, store reodering operations */
rdep = 0; /* number of dependent rows in H*/
d = 0;    /* current diagonal element */

while( (d+rdep) < M) { /* cycle through independent rows of H */

j = d; /* current column index along row ir[d] */
while( (HH[ir[d]][jc[j]] == 0) && (j<(N-1)) )
j++;            /* find first nonzero element in row i */
if( HH[ir[d]][jc[j]] ) { /* found nonzero element */

/* swap columns */
tmp = jc[d]; jc[d] = jc[j]; jc[j] = tmp;

if(q==2) { /* GF2 */
/* eliminate current column using row operations */
for(ii=0 ; ii<M ; ii++)
if(HH[ir[ii]][jc[d]] && (ii != d)) /* nonzero and non-diagonal */
for(jj=d ; jj<N ; jj++)
HH[ir[ii]][jc[jj]] = (HH[ir[ii]][jc[jj]]+HH[ir[d]][jc[jj]])%2;
}
else { /* GFq */

scale = GFq_inv(HH[ir[d]][jc[d]],q); /* inverse of the diag. element */
/* scale the current row to make the first element 1 */
for(jj=0 ; jj<N ; jj++)
HH[ir[d]][jc[jj]] = GFq_m(HH[ir[d]][jc[jj]],scale,q);

/* eliminate current column using row operations */
for(ii=0 ; ii<M ; ii++) {
if(HH[ir[ii]][jc[d]] && (ii != d)) {
scale = HH[ir[ii]][jc[d]];
for(jj=d ; jj<N ; jj++) {
tmp = GFq_m(HH[ir[d]][jc[jj]],scale,q);
HH[ir[ii]][jc[jj]] = GFq_a(HH[ir[ii]][jc[jj]],tmp);
}
}
}
}
}
else { /* all zeros -  need to delete this row and update indices */
rdep++; /* increase number of dependent rows */
tmp = ir[d];
ir[d] = ir[M-rdep];
ir[M-rdep] = tmp;
d--; /* no diagonal element is found */
}
d++; /* increase the number of diagonal elements */

}/*while i+rdep*/
/* done inverting HH */

K = N-M+rdep; /* true K */

/* create G matrix  G = [A'| I] if H = [I|A]*/
GG = (unsigned char **)mxMalloc(K*sizeof(unsigned char *));
for(i=0 ; i<K ; i++){
GG[i] = (unsigned char *)mxMalloc(N*sizeof(unsigned char));
}
for(i=0 ; i<K ; i++)
for(j=0 ; j<(N-K) ; j++) {
tmp = (N-K+i);
GG[i][j] = HH[ir[j]][jc[tmp]];
}

for(i=0 ; i<K ; i++)
for(j=(N-K); j<N ; j++)
if(i == (j-N+K) ) /* diagonal */
GG[i][j] = 1;
else
GG[i][j] = 0;

/* NOTE, it is a very inefficient way to store G. Change to taste!*/
G_OUT = mxCreateDoubleMatrix(K, N, mxREAL);
/* Assign pointers to the output matrix */
g = mxGetPr(G_OUT);
for(i=0 ; i<K ; i++)
for(j=0 ; j<N; j++)
g[i+j*K] = GG[i][j];

H_OUT = mxCreateSparse(M,N,nz,mxREAL);
sr2  = mxGetPr(H_OUT);
irs2 = mxGetIr(H_OUT);  /* row */
jcs2 = mxGetJc(H_OUT);  /* column */
/* Write H_OUT swapping columns according to jc */
k = 0;
for (j=0; (j<N ); j++) {
jcs2[j] = k;
tmp = jcs1[jc[j]+1]-jcs1[jc[j]];
for (i=0; i<tmp ; i++) {
kk = jcs1[jc[j]]+i;
sr2[k] = sr1[kk];
irs2[k] = irs1[kk];
k++;
}
}
jcs2[N] = k;

/* free the memory */
for( j=0 ; j<M ; j++) {
mxFree(HH[j]);
}
mxFree(HH);
mxFree(ir);
mxFree(jc);
for(i=0;i<K;i++){
mxFree(GG[i]);
}
mxFree(GG);
return;
}