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;
}