I have few very large arrays which I want to use to store some calculations. The arrays are of size 10000003,25003*3, 1000000,1000000.I have used malloc to dynamically allocate memory to all the arrays before using any of them.
The code tries to implement an algorithm for a discritized space. I have tried using smaller sizes for my arrays and everything seems to be working fine.
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
double mu_a=0.0;
double mu_b=0.0;
int n_b_d1=50.0;
int n_b_d2=50.0;
double dot(int m,int n,int o,int p,double * array[]){
m=m-1;
n=n-1;
o=o-1;
p=p-1;
return ((array[m][0]-array[n][0])*(array[o][0]-array[p][0]))+((array[m][1]-array[n][1])*(array[o][1]-array[p][1]))+((array[m][2]-array[n][2])*(array[o][2]-array[p][2]));
}
double * linspace(double x1, double x2,int n,double array[]){
double * x=calloc(n,sizeof(double));
double step=(x2-x1)/(double)(n-1);
for (int i=0;i<n;i++){
array[i]=x1+((double)i*step);
x[i]=x1+((double)i*step);
}
return x;
}
/*array length=*(&array+1)-array)*/
void points(double P_1[],double P_2[],double P_3[],double P_4[],double temp[]){
double dst=400*(((P_3[0]-P_1[0])*((P_2[1]-P_1[1])-(P_4[1]-P_3[1])))-((P_3[1]-P_1[1])*((P_2[0]-P_1[0])-(P_4[0]-P_3[0]))));
if (dst>0.001){
temp[3]=999999; return;
}
double * array[4]={P_1,P_2,P_3,P_4};
mu_a=((dot(1,3,4,3,array)*dot(4,3,2,1,array))-(dot(1,3,2,1,array)*dot(4,3,4,3,array)))/((dot(2,1,2,1,array)*dot(4,3,4,3,array))-(dot(4,3,2,1,array)*dot(4,3,2,1,array)));
mu_b=(dot(1,3,4,3,array)+(mu_a*dot(4,3,2,1,array)))/dot(4,3,4,3,array);
double P_a[3]={0,0,0};
double P_b[3]={0,0,0};
for (int i=0;i<3;i++){
P_a[i]=P_1[i]+mu_a*(P_2[i]-P_1[i]);
P_b[i]=P_3[i]+mu_b*(P_4[i]-P_3[i]);
}
//printf("%lf,%lf,%lf,n",P_a[0],P_a[1],P_a[2]);
dst=pow((pow((P_a[0]-P_b[0]),(double)(2))+pow((P_a[1]-P_b[1]),(double)(2))+pow(((P_a[2]-P_b[2])/10),(double)(2))),(0.5));
//printf("%lf",dst);
temp[0]=(P_a[0]+P_b[0])/2;
temp[1]=(P_a[1]+P_b[1])/2;
temp[2]=(P_a[2]+P_b[2])/2;
temp[3]=dst;
}
double omega(double a,double b,double d){
double alpha=a/(double)2/d;
double beta=b/(double)2/d;
return (a*b/pow(d,(double)2));
}
double omega_offaxis(double A,double B,double a, double b,double d){
return (omega(2*(A+a),2*(B+b),d)-omega(2*A,2*(B+b),d)-omega(2*(A+a),2*B,d)+omega(2*A,2*B,d))/(double)4;
}
double calc_prob(double b[],double d1[][3],double d2[][3]){
double d[2][3];
double dist;
if (b[2]<=200){
/*#tot=omega_4quad(b[0],b[1],50,50,b[2])*/
dist=400.0-b[2];
for(int i=0;i<2;i++){
for (int j=0;j<3;j++){
d[i][j]=d2[i][j];
}
}
}
else{
/*#tot=omega_4quad(b[0],b[1],50,50,400-b[2])*/
//double d[2][3];
for(int i=0;i<2;i++){
for (int j=0;j<3;j++){
d[i][j]=d1[i][j];
}
}
dist=b[2];
}
double x1=d1[0][0],y1=d1[0][1],z1=d1[0][2];
double x2=d1[1][0],y2=d1[1][1],z2=d1[1][2];
double x3=d2[0][0],y3=d2[0][1],z3=d2[0][2];
double x4=d2[1][0],y4=d2[1][1],z4=d2[1][2];
double x_c_1=((x3-x1)*b[2]/(double)(400))+x1;
double x_c_2=((x4-x2)*b[2]/(double)(400))+x2;
double y_c_1=((y3-y1)*b[2]/(double)(400))+y1;
double y_c_2=((y4-y2)*b[2]/(double)(400))+y2;
if ((b[0]<x_c_1) || (b[0]>x_c_2) || (b[1]<y_c_1) || (b[1]>y_c_2)){
return 0;
}
double x_a=d[0][0],y_a=d[0][1],z_a=d[0][2];
double x_b=d[1][0],y_b=d[1][1],z_b=d[1][2];
double a1=10;
double b1=10;
return (omega_offaxis(x_a-b[0],y_a-b[1],a1, b1, dist)/(double)4*M_PI);
}
//double blocks[1000000][3];
void main(){
double ** blocks=(double **)malloc(1000000*sizeof(double *));
for(int i=0;i<1000000;i++)
{
blocks[i]=(double*)malloc(3*sizeof(double));
}
double ** intersect=(double **)malloc(100000*sizeof(double *));
for(int i=0;i<100000;i++)
{
intersect[i]=(double*)malloc(3*sizeof(double));
}
//double detectors[2500][3][3]
double *** detectors=(double ***)malloc(2500*sizeof(double **));
for(int i=0;i<2500;i++)
{
detectors[i]=(double**)malloc(3*sizeof(double*));
for(int j=0;j<3;j++){detectors[i][j]=(double*)malloc(3*sizeof(double));
}}
printf("here1");
printf("%lf",M_PI);
//double blocking[6][6][41];
double x_blocking[51];
double y_blocking[51];
double z_blocking[401];
linspace(0.0,50.0,51,x_blocking);
linspace(0.0,50.0,51,y_blocking);
linspace(0.0,400.0,401,z_blocking);
//printf("%f,%f,%f",blocking[1][4][5],blocking[3][4][5],blocking[1][2][3]);
FILE * filo;
filo=fopen("input.txt","r");
double data[52822][6];
if (filo==NULL){
printf("%s","noice");
}
printf("%f",filo);
for(int i=0;i<52822;i++){
//double y=200;
fscanf(filo,"%lf %lf %lf %lfn",&data[i][0],&data[i][1],&data[i][3],&data[i][4]);
data[i][2]=0;
data[i][5]=400;
}
printf("here2");
FILE * gilo;
gilo=fopen("output.txt","w");
int kl =0;
for (int i=0;i<52822;i++){
for (int j=i;j<52822;j++){
double temp[4];
double P_1[3]={data[i][0],data[i][1],data[i][2]};
double P_2[3]={data[i][3],data[i][4],data[i][5]};
double P_3[3]={data[j][0],data[j][1],data[j][2]};
double P_4[3]={data[j][3],data[j][4],data[j][5]};
// """ if project_xy(P_1,P_2,P_3,P_4) && project_xz(P_1,P_2,P_3,P_4) && ////////project_yz(P_1,P_2,P_3,P_4):"""
points(P_1,P_2,P_3,P_4,temp);
//printf("%lfn",temp[3]);
if (temp[3]<=0.0001){
//printf("nice||||");
//#xs.append((temp[0][0]+temp[1][0])/2)
//#ys.append((temp[0][1]+temp[1][1])/2)
//#zs.append((temp[0][2]+temp[1][2])/2)
intersect[kl][0]=temp[0];
intersect[kl][1]=temp[1];
intersect[kl][2]=temp[2];
kl++;
//fprintf(gilo,"%lf,%lf,%lf",intersect[kl][0],intersect[kl][1],intersect[kl][2],intersect[kl][3]);
}}}
fprintf(gilo,"//////////////////////////////START\\\\\\\\\\\\\\\\\\\");
//printf("%f,%f,%f",data[4][0],data[4][1],data[4][5]);
printf("here3");;
int h=0;
int h1=0;
for(int i=1;i<51;i++){
for (int j=1;j<51;j++){
for (int k=1;k<401;k++){
blocks[h][0]=(x_blocking[i-1]+x_blocking[i])/2;
blocks[h][1]=(y_blocking[j-1]+y_blocking[j])/2;
blocks[h][2]=(z_blocking[k-1]+z_blocking[k])/2;
h++;
if (k==1 || k==400){
if (k==400) {k+=1;}
for(int index=-1;index<=0;index++){
detectors[h1][index+1][0]=x_blocking[i+index];
detectors[h1][index+1][1]=y_blocking[j+index];
detectors[h1][index+1][2]=z_blocking[k-1];
}
h1++;
}
}
}}
/*
for (int i=0;i<h;i++){
printf("i,%lfn",calc_prob(blocks[i],detectors[1],detectors[49]));
}
*/
/*fprintf(gilo,"||||||detectors||||||");
fprintf(gilo,"n");
for (int i=0;i<h1;i++){
fprintf(gilo,"[%lf,%lf,%lf],[%lf,%lf,%lf]",detectors[i][0][0],detectors[i][0][1],detectors[i][0][2],detectors[i][1][0],detectors[i][1][1],detectors[i][1][2]);
fprintf(gilo,"n");
}
fprintf(gilo,"||||||blocks||||||");
fprintf(gilo,"n");
for (int i=0;i<h;i++){
fprintf(gilo,"[%lf,%lf,%lf]",blocks[i][0],blocks[i][1],blocks[i][2]);
fprintf(gilo,"n");
}
*/
///////////////////////////////////////////////////
printf("here4");
int nstar(int i, int j){
int tempu=0;
for(int k=0;k<h;k++){
if( (data[k][0]>detectors[i][0][0]) &&( data[k][0]<detectors[i][1][0]) && (data[k][3]>detectors[i][0][1]) && (data[k][3]<detectors[i][1][1])){
if( (data[k][1]>detectors[j][0][0]) &&(data[k][1]<detectors[j][1][0]) && (data[k][4]>detectors[j][0][1] )&&( data[k][4]<detectors[j][1][1])){
tempu++;
}}
}
return tempu;
}
//double lambda_0ld[h];
double * lambda_0ld=(double *)malloc(h*sizeof(double ));
//double Lambda[h];
double * Lambda=(double *)malloc(h*sizeof(double ));
void classify(){
for (int j=0;j<h;j++){
lambda_0ld[j]=0;
for (int i=0;i<kl;i++){
if((intersect[i][0]<(blocks[j][0]+0.5))&&(intersect[i][0]>(blocks[j][0]-0.5))&&(intersect[i][1]<(blocks[j][1]+0.5))&&(intersect[i][1]>(blocks[j][1]-0.5))&&(intersect[i][2]<(blocks[j][2]+0.5))&&(intersect[i][2]>(blocks[j][2]-0.5))){
lambda_0ld[j]++;};
}}
}
classify();
printf("here5");
double lambda(int b){double sum=0;
classify();
double var;
for(int i=0;i<h1;i++){
for(int j=i;j<h1;j++){
if (detectors[i][2]==detectors[j][2]) continue;
double temp=0.0;
for (int p=0;p<h;p++){var=0.0;
//if (detectors[i][2]!=detectors[j][2]){
var=(lambda_0ld[p]*calc_prob(blocks[p],detectors[i],detectors[j]));
//printf("n%lf,%lf,%lfn",calc_prob(blocks[p],detectors[i],detectors[j]),lambda_0ld[p],var);
temp=temp+var;
}// if (temp!=temp) printf("popdsfoopds");
sum=sum+((double)(nstar(i,j)*calc_prob(blocks[b],detectors[i],detectors[j]))/(double)temp);}}
return (double)lambda_0ld[b]*(sum);
}
//printf("%lf",lambda_0ld[j];)
fprintf(gilo,"//////////////////////////////MLEM\\\\\\\\\\\\\\\\\\\n");
//printf("%lf,%lfn",lambda(20),lambda_0ld[20]);
printf("here6");
for (int u=0;u<h;u++){
Lambda[u]=lambda(u);
fprintf(gilo,"%d-%lf,%lf,%lf-%lfn",u,blocks[u][0],blocks[u][1],blocks[u][2],Lambda[u]);
//printf("%lfn",lambda_0ld[u]);
}
//printf();
fclose(filo);
/*
for (int i=0;i<6;i++){
for (int j=0;j<6;j++){
for(int k=0;k<41;k++){
fprintf(gilo,"%lf,",blocking[i][j][k]);
}fprintf(gilo,"n");}
fprintf(gilo,"new");}*/
fclose(gilo);}
All the help I can find on the internet and stack overflow says dynamic allocation should fix the segmentation fault. But it does not in my case. I am returning back to C after a lot of years. Please point out if I am making any logical error while handling pointers.
PS Please ignore all of the comments. Most of them were either some unimplemented stuff or used for debugging.