I have this code which do some matrix operation.
It can be run on Microsoft Visual c++ 6.0 but i got problem when run it using Dev C++ 4.9.9.2..
it says, main must return and int
Can somebody help me?
// FDduct.cpp , finite difference. method. for duct
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <time.h>
#define NR_END 1 //the parameter NR_END is used as a number of extra storage
// locations allocated at the beginning of every vector or matrix block, simply for the purpose
// of making offset pointer references guaranteed-representable.
#define FREE_ARG char*
double p1,p2; // function for sin(px)
int ij_k(int i,int j,int Nr); // mapping from (i,j) to (k)
void k_ij(int k,int *i,int *j,int Nr); // mapping from (k) to (i,j);
double theta,dtheta,r_i[100],th_i[100],ra,r120,r12pi,pi;
double angle,radius_r12, grad_r12,radiusO;
double x[1000][1000],y[1000][1000];
double xy[1000][2]; //[I][1=x,2=y];
double H[200][200]; // convection matrix [row][column]
int N_theta,N_r;
int element[1000][3],N_e; //element[no.][1st,2nd,3rd]
int BoundSTemp[100][3+1], BoundEss[100],BoundNat[100],BoundEssN,BoundNatN,BoundConvN; // [line index][1:start, 2: stop]
int i,j,ei,ej,i_,j_,k_; //matrix[row][column]
double length(int e,int node1,int node2); //length element
double area[1000],a[100][100],a_reduce[100][100]; // area[element no.]
double x1_,x2_,x3_,y1_,y2_,y3_,temp;
double conv_coef_h,thermaCond_k,heatflux_coef_q,fix_heat;
double Bmat[200][2+1][3+1],coordinate[200][2+1][3+1]; // matrix Bmat[element][row][col] ,Bmat[2+1][3+1];//
double Kmat[1000][1000]; // Kmat[G node][G node]
double r_conv_inf[1000],r_heat_q[1000],r_Q[1000],t_inf,Q; // r_conv_inf[G node],r_Q[G node]
double *dvector(long nl, long nh);
void free_dvector(double *v, long nl, long nh);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);
void MatrixInv(double **a, int N,double **y); //matrix a will be destroyed.
void nrerror(char error_text[]);
void LinearSys(double **a, int N,double *b); // The answer x will be given back in b. original matrix A will havebeen destroyed.
double funcA(int i,int j);
double funcB(int i,int j);
double funcC(int i,int j);
double funcD(int i,int j);
double w[5],point[5][3+1];
void main(void)
{
double u[6000],u1[6000],sigma,Normal,t,alpha,dt,sqdt,p_bar,eta,log_pbar,mu;
double dx;
double p1,p2,ans;
int i,j,k,n,nn,c,counter,iteration,type,MaxLoop,i1,i2,i3,ii,jj;
double temp,temp1,temp2,temp3;
double int_r,discount,dR,epsilon,error,dif;
long iseed,iseed1;
float NormalDev(long *idum);
double func_f(double x),x_start,x_end;
struct tm *date_time;
time_t timer;
time(&timer);
date_time=gmtime(&timer);
FILE *fpo,*fpo1,*fpo2,*fpo3;
double *Vecx,*Vecb,**Amat;
Vecb=dvector(1,200);
Vecx=dvector(1,200);
Amat=dmatrix(1,200,1,200);
//initialize...
fpo=fopen("mesh.txt","w+t");
fpo1=fopen("output1.txt","w+t");
pi=3.1415926535897932384626433832795;
fprintf(fpo,"%.19s\n",asctime(date_time));
fprintf(fpo,"%i\n",date_time->tm_sec);
// discretization......
fix_heat=100.00;
heatflux_coef_q=1.0;
Q=1.0;
t_inf=80.0;
conv_coef_h=25.0;
thermaCond_k=1.0;
N_theta=7;
N_r=3;
dtheta=(double) pi/((N_theta)-1);
angle=0.0; //radian....
r120=2.0; // b
r12pi=4.0; //c
grad_r12=(r12pi-r120)/pi;
radius_r12=r120;
angle=0.0; //radian....
ra=1.5; //radius from origin to r1;a
// triangular gauss quadrature integration
w[1]=-9./32.0; w[2]=25./96.0;
w[3]=w[2];w[4]=w[2];
point[1][1]=1./3.;point[1][2]=1./3.; point[1][3]=1./3.;
k_ij(3,&i,&j,N_r); // mapping from (k) to (i,j)
printf("%s %d %d %d \n","k_ij, k,i,j",3,i,j);
k_ij(4,&i,&j,N_r); // mapping from (k) to (i,j)
printf("%s %d %d %d \n","k_ij, k,i,j",4,i,j);
k_ij(5,&i,&j,N_r); // mapping from (k) to (i,j)
printf("%s %d %d %d \n","k_ij, k,i,j",5,i,j);
for (j=1;j<=N_theta;j++) { // theta discretize.
radius_r12=r120+angle*grad_r12;
for (i=1;i<=N_r;i++) { // r12 discretize.
radiusO=ra+(i-1)*radius_r12/double(N_r-1);
//create location (x,y) using polar system.....
x[i][j]=radiusO*cos(angle);
y[i][j]=radiusO*sin(angle);
xy[ij_k(i,j,N_r)][1]=x[i][j];
xy[ij_k(i,j,N_r)][2]=y[i][j];
// x[j]=(r120+radius_r12)*cos(pi);
// y[i]=(r120+radius_r12)*sin(pi);
}//i
//.........
angle=angle+dtheta;
} //jj
fprintf(fpo1,"%s \n","the discretize coord...");
for (j=1;j<=N_r*N_theta;j++) {
fprintf(fpo,"%f %f \n",xy[j][1],xy[j][2]);
} //j
//BoundEss[100],BoundNat[100],BoundEssN,BoundNatN
//create essential boundary conditions:
for (i=1;i<=N_r;i++) {
j=1;
BoundEss[i]=ij_k(i,j,N_r);
}
for (i=1;i<=N_r;i++) {
j=N_theta;
BoundEss[i+N_r]=ij_k(i,j,N_r);
}
BoundEssN=2*N_r;
fprintf(fpo1,"%s \n","the essential bondary conditions...");
for (k=1;k<=BoundEssN;k++) {
k_ij(BoundEss[k],&i,&j,N_r);
fprintf(fpo,"%d %d %d %d \n",k,BoundEss[k],i,j);
} //j
for (i=1;i<=N_r;i++) {
for (j=1;j<=N_theta;j++) {
ii=ij_k(i,j,N_r);
if(i>1) a[ii][ij_k(i-1,j,N_r)]=a[ii][ij_k(i-1,j,N_r)]+funcA(i,j);
else a[ii][ij_k(i+1,j,N_r)]=a[ii][ij_k(i+1,j,N_r)]+funcA(i,j); //du/dr=0 -> u[0][j]=u[2][j]
a[ii][ij_k(i,j,N_r)]=a[ii][ij_k(i,j,N_r)]-2.0*funcA(i,j);
if(i>1) a[ii][ij_k(i-1,j,N_r)]=a[ii][ij_k(i-1,j,N_r)]-funcB(i,j);
else a[ii][ij_k(i+1,j,N_r)]=a[ii][ij_k(i+1,j,N_r)]-funcB(i,j); //du/dr=0 -> u[0][j]=u[2][j]
if(j>1) a[ii][ij_k(i,j-1,N_r)]=a[ii][ij_k(i,j-1,N_r)]+funcC(i,j);
else
a[ii][ij_k(i,j,N_r)]=a[ii][ij_k(i,j,N_r)]-2.0*funcC(i,j);
a[ii][ij_k(i,j+1,N_r)]=a[ii][ij_k(i,j+1,N_r)]+funcC(i,j); //du/dr=0 -> u[i][0]=u[i][2]
if(j>1) a[ii][ij_k(i,j-1,N_r)]=a[ii][ij_k(i,j-1,N_r)]-funcD(i,j);
else
a[ii][ij_k(i,j+1,N_r)]=a[ii][ij_k(i,j+1,N_r)]+funcD(i,j); //du/dr=0 -> u[i][0]=u[i][2]
if(i>2) a[ii][ij_k(i-1,j,N_r)]=a[ii][ij_k(i-1,j,N_r)]+funcA(i,j);
else a[ii][ij_k(i+1,j,N_r)]=a[ii][ij_k(i+1,j,N_r)]+funcA(i,j); //du/dr=0 -> u[1][j]=u[3][j]
a[ii][ij_k(i,j,N_r)]=a[ii][ij_k(i,j,N_r)]-2.0*funcA(i,j);
if(i>2) a[ii][ij_k(i-1,j,N_r)]=a[ii][ij_k(i-1,j,N_r)]-funcB(i,j);
else a[ii][ij_k(i+1,j,N_r)]=a[ii][ij_k(i+1,j,N_r)]-funcB(i,j); //du/dr=0 -> u[1][j]=u[3][j]
if(j>2) a[ii][ij_k(i,j-1,N_r)]=a[ii][ij_k(i,j-1,N_r)]+funcC(i,j);
else
a[ii][ij_k(i,j,N_r)]=a[ii][ij_k(i,j,N_r)]-2.0*funcC(i,j);
a[ii][ij_k(i,j+1,N_r)]=a[ii][ij_k(i,j+1,N_r)]+funcC(i,j); //du/dr=0 -> u[i][1]=u[i][3]
if(j>1) a[ii][ij_k(i,j-1,N_r)]=a[ii][ij_k(i,j-1,N_r)]-funcD(i,j);
else
a[ii][ij_k(i,j+1,N_r)]=a[ii][ij_k(i,j+1,N_r)]+funcD(i,j); //du/dr=0 -> u[i][1]=u[i][3]
}
return 0;
}
//checking for double funcA(int i,int j)
fprintf(fpo,"%s %f \n","funcA(1,1)",funcA(1,1));
//checking for double funcA(int i,int j)
fprintf(fpo,"%s %f \n","funcA(3,2)",funcA(3,2));
////checking for double funcB(int i,int j)
fprintf(fpo,"%s %f \n","funcB(1,1)",funcB(1,1));
////checking for double funcB(int i,int j)
fprintf(fpo,"%s %f \n","funcB(2,2)",funcB(2,2));
////checking for double funcC(int i,int j)
fprintf(fpo,"%s %f \n","funcC(3,2)",funcC(3,2));
////checking for double funcD(int i,int j)
fprintf(fpo,"%s %f \n","funcD(1,1)",funcD(1,1));
////checking for double funcD(int i,int j)
fprintf(fpo,"%s %f \n","funcD(3,2)",funcD(3,2));
// print Vecb[i] (new)
fprintf(fpo,"\n Vecb[i](new): \n");
for (i=1; i<=N_theta*N_r; i++)
{
fprintf(fpo,"%.4f ",Vecb[i]);
fprintf(fpo,"\n");
}//i
exit(0);
fprintf(fpo," %s \n","We are using Gauss Elimination.. ");
//LinearSys(Amat,N_theta*N_r,Vecb); //Amat will be destroyed. ans return in Vecb.
fprintf(fpo,"\n Final answer: \n");
for (k=1; k<=N_theta*N_r; k++) //global node.
{
// k_ij(k,&i,&j,N_r);
xy[ k][1]; xy[ k][2];
fprintf(fpo,"%.4f %.4f %.4f ",xy[ k][1], xy[ k][2], Vecb[k]);
fprintf(fpo,"\n");
}//i
free_dvector(Vecb,1, 200);
free_dvector(Vecx,1, 200);
free_dmatrix(Amat,1,200,1,200);
fclose(fpo);
fclose(fpo1);
} // end. .....................main..................
double funcA(int i,int j)
{
double temp,thetaj,ri,m,b,dr;
thetaj=(double) (j-1.0)*dtheta;
ri=(double) (i-1.0)/(N_r-1.0);
b=r120;
m=grad_r12;
dr=(double) 1.0/(N_r-1.0);
temp=(m*thetaj+b)*pow(dr,2);
return 1.0/temp;
}// ....................................
double funcB(int i,int j)
{
double temp,temp2,temp3,thetaj,ri,m,b,dr,R;
thetaj=(double) (j-1.0)*dtheta;
ri=(double) (i-1.0)/(N_r-1.0);
b=r120;
m=grad_r12;
R=(double) ra+ri*(m*thetaj+b);
dr=(double) 1.0/(N_r-1.0);
temp=R*(m*thetaj+b)+(ri*pow(m,2));
temp2=pow(m*thetaj+b,2)*(pow(R,2)+pow(ri*m,2));
temp3= (temp/temp2)*(1/(2.0*dr));
return temp3;
} //....................
double funcC(int i,int j)
{
double temp,thetaj,ri,m,b,R,dr;
thetaj=(double) (j-1.0)*dtheta;
ri=(double) (i-1.0)/(N_r-1.0);
b=r120;
m=grad_r12;
R=(double) ra+ri*(m*thetaj+b);
dr=(double) 1.0/(N_r-1.0);
temp=(pow(R,2)+pow(ri*m,2))*pow(dtheta,2);
return 1.0/temp;
}// ................
double funcD(int i,int j)
{
double temp,temp2,temp3,thetaj,ri,m,b,R,dr;
thetaj=(double) (j-1.0)*dtheta;
ri=(double) (i-1.0)/(N_r-1.0);
b=r120;
m=grad_r12;
R=(double) ra+ri*(m*thetaj+b);
dr=(double) 1.0/(N_r-1.0);
temp=m-((R*ri*m)*(m*thetaj+b));
temp2=(m*thetaj+b)*(pow(R,2)+pow(ri*m,2));
temp3=(temp/temp2)*(1/2*dtheta);
return temp3;
}// ................
double length(int e,int node1,int node2)
{
double temp,x1,y1,x2,y2;
x1=xy[element[e][node1]][1];
y1=xy[element[e][node1]][2];
x2=xy[element[e][node2]][1];
y2=xy[element[e][node2]][2];
temp=pow(x2-x1,2)+pow(y2-y1,2);
return sqrt(temp);
}// ....................................
int ij_k(int i,int j,int Nr) // mapping from (i,j) to (k), correct...
{
return (j-1)*Nr+i;
}
void k_ij(int k,int *i,int *j,int Nr) // mapping from (k) to (i,j), correct....
{
int div;
div=int(k/Nr);
*i=k%Nr; // modulo Nr
*j=div+1;
if(*i==0)
{
*i=Nr;
*j=*j-1;
}
}//...........................
double func_f(double x)
{
//double p1;
//p1=2.0;
//return sin(p1*x);
return sin(p1*x)+sin(p2*x);
// return -x;
} // func_f
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
} // ...............................
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (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)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}//...................
void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}//.....................
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
} //..........................................
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
printf("%s\n",error_text);
printf("...now exiting to system...\n");
exit(1);
}
#define TINY 1.0e-25 // 1.0e-20 // A small number.
void ludcmp(double **a, int n, int *indx, double *d) // matrix a will be destroyed...
/* Given a matrix a[1..n][1..n], this routine replaces it by the LU decomposition of a rowwise
permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above;
indx[1..n] is an output vector that records the row permutation effected by the partial
pivoting; d is output as ±1 depending on whether the number of row interchanges was even
or odd, respectively. This routine is used in combination with lubksb to solve linear equations
or invert a matrix.
*/
{
int i,imax,j,k;
double big,dum,sum,temp;
double *vv; // vv stores the implicit scaling of each row.
vv=dvector(1,n);
*d=1.0; // No row interchanges yet.
for (i=1;i<=n;i++) { // Loop over rows to get the implicit scaling information.
big=0.0;
for (j=1;j<=n;j++)
if ((temp=fabs(a[i][j])) > big) big=temp;
if (big == 0.0) nrerror("Singular matrix in routine ludcmp");
// No nonzero largest element.
vv[i]=1.0/big; //Save the scaling.
} // for i
for (j=1;j<=n;j++) { //This is the loop over columns of Crout’s method.
for (i=1;i<j;i++) { //This is equation (2.3.12) except for i = j.
sum=a[i][j];
for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
a[i][j]=sum;
} // for i
big=0.0; //Initialize for the search for largest pivot element.
for (i=j;i<=n;i++) { //This is i = j of equation (2.3.12) and i = j+1. . .N of equation (2.3.13).
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) {
//Is the figure of merit for the pivot better than the best so far?
big=dum;
imax=i;
} // if
}// for i
if (j != imax) { //Do we need to interchange rows?
for (k=1;k<=n;k++) { //Yes, do so...
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}// for k
*d = -(*d); //...and change the parity of d.
vv[imax]=vv[j]; //Also interchange the scale factor.
} // if j
indx[j]=imax;
// if (a[j][j] == 0.0) a[j][j]=TINY; //original
if (fabs(a[j][j]) < TINY) a[j][j]=TINY; //change.
/* If the pivot element is zero the matrix is singular (at least to the precision of the
algorithm). For some applications on singular matrices, it is desirable to substitute
TINY for zero. */
if (j != n) { //Now, finally, divide by the pivot element.
dum=1.0/(a[j][j]);
for (i=j+1;i<=n;i++) a[i][j] *= dum;
} // if j
}//for j //Go back for the next column in the reduction.
free_dvector(vv,1,n);
} // ...........................ludcmp .........................
void lubksb(double **a, int n, int *indx, double b[])
/*
Solves the set of n linear equations A·X = B. Here a[1..n][1..n] is input, not as the matrix
A but rather as its LU decomposition, determined by the routine ludcmp. indx[1..n] is input
as the permutation vector returned by ludcmp. b[1..n] is input as the right-hand side vector
B, and returns with the solution vector X. a, n, and indx are not modified by this routine
and can be left in place for successive calls with different right-hand sides b. This routine takes
into account the possibility that b will begin with many zero elements, so it is efficient for use
in matrix inversion. */
{
int i,ii=0,ip,j;
double sum;
for (i=1;i<=n;i++) { //When ii is set to a positive value, it will become the
// index of the first nonvanishing element of b. We now do the forward substitution, equation (2.3.6). The
// only new wrinkle is to unscramble the permutation as we go.
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; //A nonzero element was encountered, so from now on we
// will have to do the sums in the loop above.
b[i]=sum;
}// for i
for (i=n;i>=1;i--) { // Now we do the backsubstitution, equation (2.3.7).
sum=b[i];
for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
b[i]=sum/a[i][i]; //Store a component of the solution vector X.
} // for i, // All done!
} // ..................lubksb ........................
void MatrixInv(double **a, int N,double **y) // matrix a will be destroyed. correct..
{
double d,*col;
int i,j,*indx;
col=dvector(1,N);
indx=ivector(1,N);
ludcmp(a,N,indx,&d); //Decompose the matrix just once.
for(j=1;j<=N;j++) { //Find inverse by columns.
for(i=1;i<=N;i++) col[i]=0.0;
col[j]=1.0;
lubksb(a,N,indx,col);
for(i=1;i<=N;i++) y[i][j]=col[i];
} // for j
free_dvector(col,1,N);
free_ivector(indx,1,N);
}// ............................MatrixInv............................
void LinearSys(double **A, int N,double *b) // correct.
{// The answer x will be given back in b. original matrix A will havebeen destroyed.
double d;
int *indx;
indx=ivector(1,N);
ludcmp(A,N,indx,&d);
lubksb(A,N,indx,b);
free_ivector(indx,1,N);
}// ............................MatrixInv............................