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

```/* 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]

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

void mexFunction(
int nlhs,       mxArray *plhs[],
int nrhs, const mxArray *prhs[]
)
{
unsigned char **HH, **GG;
int ii, jj, *ir, *jc, rdep, tmp, d;
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 != 1) {
mexErrMsgTxt("h2g requires one input arguments.");
} else if (nlhs != 2) {
mexErrMsgTxt("h2g requires two output arguments.");
} else if (!mxIsSparse(H_IN)) {
mexErrMsgTxt("h2g requires sparse H matrix.");
}

/* 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] = 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. It is "1" in GF2 */

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

/* eliminate current column using row operations */
for(ii=0 ; ii<M ; ii++)
if(HH[ir[ii]][jc[d]] && (ir[ii] != ir[d]))
for(jj=d ; jj<N ; jj++)
HH[ir[ii]][jc[jj]] = (HH[ir[ii]][jc[jj]]+HH[ir[d]][jc[jj]])%2;
}
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 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;
}

```