Working on creating a matrix with given paremeters and then solving the matrix. The matrices A and C come out just fine, but there is a seg fault when I try to run the lud files. Can anyone help?
Thank you very much!

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#define NR_END 1
#define FREE_ARG char*

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <cmath>
#include <fstream>
#define TINY 1.0e-20;

using namespace std;

void free_vector(float *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}


float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
	float *v;

	v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
	return v-nl+NR_END;
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
	int *v;

	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
	return v-nl+NR_END;
}



void ludcmp(float **a, int dim, int *indx, float *d)
{
int i,imax,j,k;
float big,dum,sum,temp;
float *vv;

vv=vector(1,dim);
*d=1.0;

for (i=1;i<=dim;i++)
{
big=0.0;
for (j=1;j<=dim;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
if (big == 0.0)
vv[i]=1.0/big;
}

for (j=1;j<=dim;j++)
{
for (i=1;i<j;i++)
{
sum=a[i][j];
for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<=dim;i++)
{
sum=a[i][j];
for (k=1;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big)
{
big=dum;
imax=i;
}
}
if (j != imax)
{
for (k=1;k<=dim;k++)
{
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
*d = -(*d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0) a[j][j]=TINY;
if (j != dim)
{
dum=1.0/(a[j][j]);
for (i=j+1;i<=dim;i++) a[i][j] *= dum;
}
}
free_vector(vv,1,dim);
}

void lubksb(float **a, int dim, int *indx, float b[])
{
int i,ii=0,ip,j;
float sum;
for (i=1;i<=dim;i++)
{
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
if (ii)
for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
else if (sum) ii=i;
b[i]=sum;
}
for (i=dim;i>=1;i--)
{
sum=b[i];
for (j=i+1;j<=dim;j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i];

for (i=1; i<=dim; i++)
{
cout<<b[i]<<"\n";
}
}
}



int main ()
{
double lf, xmax, ybase, ymax, tf, xfin, Ab, To, Tb, P, q, ka, ha, x, y, dx, dy;
int N, n, r, c, cf;
float d, *C, **A;
int *indx;
  
indx = (int *)malloc(816 * sizeof(int));
int k;
    for (k = 0; k < 816; k++)
       indx[k] = 0;

  // Allocate memory
  A = new float*[816];
  for (int i = 0; i < 816; ++i)
    A[i] = new float[816];


C = (float *)malloc(816 * sizeof(float));
    for (k = 0; k < 816; k++)
       C[k] = 0;

cf=0;
lf=.02;
xmax=.0075;
ybase=.005;
ymax=ybase+lf;
tf=.005;
xfin=xmax-tf/2;
Ab=.2*.2;
To=25+273;
Tb=105;
P=80;
q=P/(Ab);
ka=10;
ha=100;
n=0;
x=0.0;
y=0.0;
dx=.0005;
dy=.0005;
N=xmax/dx+1;

while ((y<=ymax+dx)&& (x<=xmax+dx)){
    while((x<=xmax+dx/5)){
//         %Section I
        if ((x==0)&&(y==0)){
            A[n][n+1]=1.0;
            A[n][n+N]=1.0;
            A[n][n]=-2.0;
            C[n]=-q*dx/ka;
            n=n+1;
		}
//         %Section 2
        if ((x>0)&&(x<xmax)&&(y==0)){
            A[n][n-1]=1.0;
            A[n][n+1]=1.0;
            A[n][n+N]=2.0;
            A[n][n]=-4.0;
            C[n]=-2*q*dx/ka;
            n=n+1;
		}
//         %Section 3
        if ((x>=xmax)&&(y==0.0)){
            A[n][n-1]=1.0;
            A[n][n+N]=1.0;
            A[n][n]=-2.0;
            C[n]=-q*dx/ka;
            n=n+1;
		}
//         %Section 4
        if ((x==0)&&(y>0)&&(y<ybase)){
            A[n][n-N]=1.0;
            A[n][n+N]=1.0;
            A[n][n+1]=2.0;
            A[n][n]=-4.0;
            n=n+1;
		}
//         %Section 5
        if (((x>0)&& (x<xmax)&&(y>0)&&(y<ybase))||((x>xfin+dx/2)&&(x<xmax)&&(y>=ybase-dy/2)&&(y<ymax))){
            A[n][n-N]=1.0;
            A[n][n+N]=1.0;
            A[n][n+1]=1.0;
            A[n][n-1]=1.0;
            A[n][n]=-4.0;
            n=n+1;
		}
//         %Section 6
        if ((x>=xmax)&&(y>0)&&(y<ymax)){
            A[n][n-N]=1.0;
            A[n][n+N]=1.0;
            A[n][n-1]=2.0;
            A[n][n]=-4.0;
            n=n+1;
		}
//         %Section 7
        if ((x==0)&&(y>ybase-dy/5)&&(y<ybase+dy/5)){
            A[n][n-N]=1.0;
            A[n][n+1]=1.0;
            A[n][n]=-2.0*(ha*dx/ka+1);
            C[n]=-2.0*ha*To*dx/ka;
            n=n+1;
		}
//         %Section 8
        if ((x>0+dx/2)&&(x<xfin)&&(y>ybase-dy/5)&&(y<ybase+dy/5)){
            A[n][n-N]=2.0;
            A[n][n+1]=1.0;
            A[n][n-1]=1.0;
            A[n][n]=-2.0*(2.0+ha*dx/ka);
            C[n]=-2.0*ha*To*dx/ka;
            n=n+1;
		}
//         %Section 9
        if ((x>xfin-dx/5)&&(x<xfin+dx/5)&&(y>ybase-dy/5)&&(y<ybase+dy/5)){
            A[n][n-N]=2.0;
            A[n][n+1]=2.0;
            A[n][n-1]=1.0;
            A[n][n+N]=1.0;
            A[n][n]=-2.0*(3.0+ha*dx/ka);
            C[n]=-2.0*ha*To*dx/ka;
            n=n+1;
		}
//         %Section 10
        if ((x>xfin-dx/5)&&(x<xfin+dx/5)&&(y>ybase+dx/2)&&(y<ymax)){
            A[n][n-N]=1.0;
            A[n][n+1]=2.0;
            A[n][n+N]=1.0;
            A[n][n]=-2.0*(2+ha*dx/ka);
            C[n]=-2.0*ha*To*dx/ka;
            n=n+1;
		}
//         %Section 11
        if ((x>xfin-dx/5)&&(x<xfin+dx/5)&&(y>=ymax)){
            A[n][n-N]=1.0;
            A[n][n+N]=1.0;
            A[n][n]=-2.0*(1.0+ha*dx/ka);
            C[n]=-2.0*ha*To*dx/ka;
            n=n+1;
		}
//         %Section 12
        if ((x>xfin+dx/2)&&(x<xmax)&&(y>=ymax)){
            A[n][n-N]=2.0;
            A[n][n+N]=1.0;
            A[n][n-1]=1.0;
            A[n][n]=-2.0*(2.0+ha*dx/ka);
            C[n]=-2*ha*To*dx/ka;
            n=n+1;
		}
//         %Section 13
        if ((x>=xmax)&&(y>=ymax)){
            A[n][n-N]=1.0;
            A[n][n-1]=1.0;
            A[n][n]=-(2.0+ha*dx/ka);
            C[n]=-ha*To*dx/ka;
            n=n+1;
		}
//         %Section 14
        if ((x>=0)&&(x<xfin)&&(y>ybase+dy/2)&&(y<ymax+dy/5)){
            A[n][n]=1.0;
            C[n]=To;
            n=n+1;
		}
        x=x+dx;
	}
    x=0.0;
    y=y+dy;
    cf=cf+1;


}

ludcmp(A,n,indx,&d);
lubksb(A,n,indx,C);
cout<<"The A matrix is: \n";
// while (c<n){
// 	while(r<n){
// 		cout<<A[r][c]<<" ";
// 		r=r+1;
// 	}
// 	cout<<"\n";
// 	r=0;
// 	c=c+1;
// }
c=0;
cout<<"The C matrix is: \n";
while (c<n){
	cout<<indx[c]<<" ";
	c++;
}
return 0;
}

Recommended Answers

All 2 Replies

> v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
> return v-nl+NR_END;
Do you want to explain these?

Problems include
- meaningless variable names
- poorly to non-existent indentation
- a chaotic mix of C and C++ code.

First thing do is work out where the problem is.

For a Windows build run it in debugging mode in Visual studio.

For a unix-like or cygwin GNU build generate the core file by executing:

ulimit -c unlimited

and then running a debug build of your program.

Then open it in the gnu debugger:

gdb <executable_name> <core.file>

Type 'wh' to give you your source and type 'up' a few times (as you're likely in the standard library) until you get to your code.

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.