[Rd] Hello! I got error in C - R

[kentjin at stat.tamu.edu]

Hello,
My name is Ick Hoon Jin and I am Ph. D. student in Texas A & M Univ..
When I run the C embedded in R in the Linux system, I confront the
following error after 6,000 iteration. By googling I found this error is
from the problem in C.

*** caught segfault ***
address (nil), cause 'memory not mapped'

My C code is following:

---------------------------------------------------------------------------
#include<stdio.h>
#include<stdlib.h>
#include<R.h>
#include<Rmath.h>
#include<malloc.h>

#define Min(X,Y) ((X) < (Y) ? (X) : (Y))

void SecondMH(int *n, int *X, int *length_X, double *theta_new, int
*n_para, double *out);
int MH_Result(double *theta_new, int n_para, int *M1, int *M2, int n);
double Count_Edges(int *X, int n);
double Cal_GWD(int *X, int n);
int Cal_Degree(int *M3, int n);
double Inner_Product( double *vector1, double *vector2, int length);
int rdtsc();

void SecondMH(int *n, int *X, int *length_X, double *theta_new, int
*n_para, double *out)
{
	int i, j, a, k;
	int *M1, *M2;
	M1 = (int *)calloc(*length_X, sizeof(int));
	M2 = (int *)calloc(*length_X, sizeof(int));

	for(i = 1; i < *n; i++)
	{
		for(j = 0; j < i; j++)
		{
			for(k = 0; k < *length_X; k++ )
			{
				M1[k] = X[k];
				M2[k] = X[k];
			}

			if(X[i * *n + j] == 1)
			{
				M1[i * *n + j] = 0;
				a = MH_Result(theta_new, *n_para, M1, M2, *n);
				if(a == 1)
				{
					X[i * *n + j] = 0;
				}
			}
			else
			{
				M1[i * *n + j] = 1;
				a = MH_Result(theta_new, *n_para, M1, M2, *n);
				if(a == 1)
				{
					X[i * *n + j] = 1;
				}
			}
		}
	}

	for(i = 1; i < *n; i++)
	{
		for(j = 0; j < i; j++)
		{
			X[j * *n + i] = 0;
		}
	}

	for(i = 0; i < *length_X; i++)
	{
		out[i] = (double)X[i];
	}

	free(M1);
	free(M2);

	return;
}

int MH_Result(double *theta_new, int n_para, int *M1, int *M2, int n)
{
	double *M1_STAT, *M2_STAT;
	double pi_Num, pi_Denom, MH_Ratio, v;

	M1_STAT = (double *)calloc( n_para, sizeof( double ) );
	M2_STAT = (double *)calloc( n_para, sizeof( double ) );

	M1_STAT[0] = Count_Edges(M1, n);
	M2_STAT[0] = Count_Edges(M2, n);

	M1_STAT[1] = Cal_GWD(M1, n);
	M2_STAT[1] = Cal_GWD(M2, n);

	pi_Num = Inner_Product(theta_new, M1_STAT, n_para);
	pi_Denom = Inner_Product(theta_new, M2_STAT, n_para);
	MH_Ratio = pi_Num - pi_Denom;

	srand(rdtsc());
	v = (double)rand() / ( (double)RAND_MAX + (double)1 );

	if( log( v ) < Min( 0, MH_Ratio ) )
		return 1;
	else
		return 0;

	free(M1_STAT);
	free(M2_STAT);
}

double Count_Edges(int *X, int n)
{
	double temp;
	int i, j;

	temp = 0;

	for(i = 1; i < n; i++)
	{
		for(j = 0; j < i; j++)
		{
			temp += X[i * n + j];
		}
	}

	return temp;
}

double Cal_GWD(int *X, int n)
{

	int *M3, *Degree, *D_Y;
	int i, j;
	double theta = 0.25;
	double GWD = 0.00;

	M3 = (int *)calloc(n, sizeof(int));
	Degree = (int *)calloc(n, sizeof(int));
	D_Y = (int *)calloc((n - 1), sizeof(int));

	for(i = 0; i < n; i++)
	{
		for(j = 0; j < n; j++)
		{
			M3[j] = X[i * n + j];
		}

		Degree[i] = Cal_Degree(M3, n);
	}

	for(i = 0; i < n - 1; i++)
	{
		D_Y[i] = 0;

		for(j = 0; j < n; j++)
		{
			if(Degree[j] = i + 1)
			{
				D_Y[i] += 1;
			}
		}
	}

	for(i = 0; i < n-1; i++)
	{
			GWD += ( 1 - pow( (1 - exp(theta*(-1))), (i + 1) ) ) * (double)D_Y[i];
	}

	GWD = GWD * exp(theta);

	free(M3);
	free(Degree);
	free(D_Y);

	return GWD;
}

int Cal_Degree(int *M3, int n)
{
	int i, result;

	result = 0;

	for(i = 0; i < n; i++)
	{
		result += M3[i];
	}

	return result;
}

double Inner_Product(double *vector1, double *vector2, int length)
{
	int i;
	double result;

	result = 0;
	for(i = 0; i < length; i++)
	{
		result += vector1[i] * vector2[i];
	}

	return result;
}

int rdtsc()
{
    __asm__ __volatile__("rdtsc");
}
-------------------------------------------------------------------------

and my R-code to call C function is following;

#-------------------------------------------------------------------------------

# Function for Second Gibbs with Metropolis-Hasting Step

#-------------------------------------------------------------------------------

dyn.load("SecondMH.so")

second.MH<-function(Adj.Matrix,theta.new,nodes){

  n<-as.integer(nodes)

  Y<-as.matrix.network(X,matrix.type="adjacency",directed=FALSE)

  Z<-as.integer(as.vector(Y))

  length.Z<-as.integer(length(Z))

  theta.star<-as.double(as.vector(theta.new))

  n.para<-as.integer(length(theta.star))

  result<-.C("SecondMH",n,Z,length.Z,theta.star,n.para,out=as.double(rep(0.0,length.Z)))

  X.new<-matrix(result$out,n,n)
  Adj.Matrix<-network(X.new, directed=FALSE)

  return(Adj.Matrix)}


--------------------------------------------------------------------------

Until this time, I have used only R and I am really a beginner for C
programming. How can I solve the segfault error? Please answer it.
Thank you very much.
Sincerely,



Jin, Ick Hoon
Ph. D. Student, Department of Statistics
Texas A&M University, College Station, TX