peg.c

density evolution for LDPC codes design and some related documents

/*
  	FILENAME: peg.c
  	AUTHOR: Tadashi Wadayama

	HOW TO MAKE:
	gcc -O2 -o peg peg.c

*/

#include <stdio.h>
#include <stdlib.h>

void visit_cnode(int c, int v, int l);

int n;
int m;
int wc;
int **H;
int *rweight;			/* row weight of H */
int L;
int *visited;
int len;

void new_column()
{
  int i,j;
  int min, mindex;
  for (i = 0; i <= m-1; i++) {
    H[i][n] = 0;
  }
  for (j = 1; j <= wc; j++) {
    min = 10000;
    for (i = 0; i <= m-1; i++) {
      if (rweight[i] < min) {
	min = rweight[i];
	mindex = i;
      }
    }
    H[mindex][n] = 1;
    rweight[mindex]++;
  }
  n++;
}

void visit_vnode(int c, int v, int l)
{
  int i,j;
  if (l > L) return;
  //for (i = 1; i <=l; i++) printf("\t");
  //printf("visit_vnode: v = %d, l = %d\n",v,l);
  for (i = 0; i <= m-1; i++) {
    if (H[i][v] == 1) visit_cnode(i,v,l+1);
  }
}


void visit_cnode(int c, int v, int l)
{
  int i,j;
  if (l > L) return;
  //for (i = 1; i <=l; i++) printf("\t");
  //printf("visit_cnode: c = %d, l = %d\n",c,l);
  visited[c] = 1;
  for (i = 0; i <= n; i++) {
    if (H[c][i] == 1) visit_vnode(c,i,l+1);
  }
}

int peg_process()
{
  int i,j;
  int min, mindex;
  int sum;

  for (i = 0; i <= m-1; i++) {
    H[i][n] = 0;
  }
  for (i = 0; i <= m-1; i++) visited[i] = 0;

  for (j = 1; j <= wc; j++) {
    //printf("==================== j = %d\n",j);
    min = 10000;
    for (i = 0; i <= m-1; i++) {
      if ((rweight[i] < min)&&(visited[i] == 0)) {
	min = rweight[i];
	mindex = i;
      }
    }
    H[mindex][n] = 1;
    rweight[mindex]++;
    for (i = 0; i <= m-1; i++) visited[i] = 0;
    visit_cnode(mindex,n,0);
    
    //printf("visited:");
    //for (i = 0; i <= m-1; i++) printf("%d ",visited[i]);
    //printf("\n");
    
    sum = 0;
    for (i = 0; i <= m-1; i++) sum += visited[i];
    //printf("number of visited checknodes = %d\n",sum);
    if (sum == m) {
      fprintf(stderr,"PEG process stoped: n = %d\n",n);
      return -1;
    }
  }
  n++;  
}


void print_H()
{
  int i,j;
  for (i = 0; i <= m-1; i++) {
    for (j = 0; j <= n-1; j++) {
      printf("%d ",H[i][j]);
    }
    printf("\n");
  }
}

void init()
{
  int i,j;
  H = malloc(sizeof(int*)*m);
  rweight = malloc(sizeof(int)*m);
  visited = malloc(sizeof(int)*m);
  for (j = 0; j <= m-1; j++) {
    H[j] = malloc(sizeof(int)*5000);
    rweight[j] = 0;
  }
}


void output_spmat()
{
  int *row_w;
  int *col_w;
  int i,j;
  int w;
  int max_rw,max_cw;

  row_w = malloc(sizeof(int*)*m);
  col_w = malloc(sizeof(int*)*n);

  max_rw = 0;
  for (i = 0; i <= m-1; i++) {
    w = 0;
    for (j = 0; j <= n-1; j++) {
      w += H[i][j];
    }
    row_w[i] = w;
    if (w > max_rw) max_rw = w;
  }

  max_cw = 0;
  for (j = 0; j <= n-1; j++) {
    w = 0;
    for (i = 0; i <= m-1; i++) {
      w += H[i][j];
    }
    col_w[j] = w;
    if (w > max_cw) max_cw = w;
  }

  printf("%d %d\n",n,m);
  printf("%d %d\n",max_rw,max_cw);
  for (i = 0; i<=m-1; i++) printf("%d ",row_w[i]);
  printf("\n");
  for (j = 0; j<=n-1; j++) printf("%d ",col_w[j]);
  printf("\n");

  for (i = 0; i <= m-1; i++) {
    for (j = 0; j <= n-1; j++) {
      if (H[i][j] == 1) printf("%d ",j+1);
    }
    printf("\n");
  }
}

int main(int argc,char **argv)
{
  int i,j;
  int r;

  if (argc != 5) {
    fprintf(stderr,"usage : peg m wc girth len\n");
    fprintf(stderr,"m: number of rows\n");
    fprintf(stderr,"wc: column weight\n");
    fprintf(stderr,"girth: girth of the bipartite graph\n");
    fprintf(stderr,"len: code length\n");
    exit(-1);
  }
  m = atoi(argv[1]);
  wc = atoi(argv[2]);
  L = atoi(argv[3])-2;
  len = atoi(argv[4]);
  if  (L % 2 != 0) {
    fprintf(stderr,"girth must be a even number\n");
    exit(-1);
  }
  
  //printf("m = %d\n",m);
  //printf("wc = %d\n",wc);
  
  n = 0;
  init();
  while(1) {
    r = peg_process();
    if (r == -1) break;
    if (n == len) break;
    fprintf(stderr,"n = %d\n",n);
  }
  //print_H();
  output_spmat();
  printf("# girth = %d\n",L+2);
}