I have attached the solver part of my code.
you can see its error.
matrix size is 10000x1000.
#include<math.h >
#include<fstrea m.h>
#include<iostre am.h>
#include<conio. h>
#include<stdlib .h>
#include <iomanip.h>
#include <stdio.h>
#include <time.h>
#include <malloc.h>
#include <stdio.h>
#include <math.h>
//************* INPUT Data *************** *************** *************** *************** ********
const int Nx=100; // no. of fine scale grid in x direction
const int Ny=100; // no. of fine scale grid in y direction
//*************** *********** Variables Arrays *************** *************** ***********
double **a,**l,**u,*b, *x;
int i,j,k,n;
//*************** ************** Subroutines and Functions *************** *************** *******/
void linbcg(unsigned long n, double b[], double x[], int itol, double tol,
int itmax, int *iter, double *err);
void solver(double **a,double **u,double **l, double *b,int n,int N); // Solver subroutine
int min(int r,int s); //
int max(int r,int s);
void main()
{
//*************** ***** Matrix Decleration *************** *************** *************** *********
n=Nx*Ny;
b=new double [n+1];
x=new double [n+1];
a=new double *[n+1];
l=new double *[n+1];
u=new double *[n+1];
for(i=0;i<=n;i+ +){
a[i]=new double [n+1];
l[i]=new double [n+1];
u[i]=new double [n+1];
}
//********* a[][]*************** ** adjacency matrix
for(j=1;j<=n;j+ +)
for(i=1;i<=n;i+ +)
a[i][j]=0.0;
j=1;
for(i=1;i<=Nx;i ++){
k=Nx*(j-1)+i;
a[k][k]=4.9;
}
for(j=2;j<Ny;j+ +)
for(i=1;i<=Nx;i ++){
k=Nx*(j-1)+i;
a[k][k]=5;
}
j=Ny;
for(i=1;i<=Nx;i ++){
k=Nx*(j-1)+i;
a[k][k]=4.9;
}
for(j=1;j<=Ny;j ++)
for(i=1;i<Nx;i+ +){
k=Nx*(j-1)+i;
a[k][k+1]=2.2;
}
for(j=1;j<=Ny;j ++)
for(i=2;i<=Nx;i ++){
k=Nx*(j-1)+i;
a[k][k-1]=2.2;
}
for(j=1;j<=Ny-1;j++)
for(i=1;i<=Nx;i ++){
k=Nx*(j-1)+i;
a[k][k+Nx]= 2.4;
}
for(j=2;j<=Ny;j ++)
for(i=1;i<=Nx;i ++){
k=Nx*(j-1)+i;
a[k][k-Nx]=2.4;
}
ofstream matrixA("A.dat" ,ios::in);
for (j=1;j<=n;j++)
{
for (i=1;i<=n;i++)
{
matrixA<<a[i][j]<<" \t";
}
matrixA<<"\n";
}
//*************** *****b[]********** Forcing function matrix *************** *************
for(j=1;j<=Ny;j ++)
for(i=1;i<=Nx;i ++){
k=Nx*(j-1)+i;
b[k]=(0.0);
}
i=1;
for(j=1;j<=Ny;j ++)
{
k=Nx*(j-1)+1;
b[k]= 6;
}
i=Nx;
for(j=1;j<=Ny;j ++)
{
k=Nx*(j-1)+Nx;
b[k]= 1;
}
//********** solver *************** *************** ***************
solver(a,u,l,b, n,Nx);
for(j=1;j<=Ny;j ++)
for(i=1;i<=Nx;i ++)
{
k=Nx*(j-1)+i;
x[k]=b[k];
}
for(j=1;j<=Ny;j ++)
for(i=1;i<=Nx;i ++)
{
k=Nx*(j-1)+i;
cout<<"x("<<k<< ","<<i<<","<<j< <"]="<<x[k]<<endl;
}
}
//*************** *************** *************** *************** *************** *************** *******
//*************** ******** subroutin Matrix solver *************** *************** *********
void solver(double **a,double **u,double **l, double *b,int n,int N)
{
int i,j,k,mi,mx;
double zig,sig;
int p=N+1; int q=N+1;
//******* LU decomposition *************** ******
for(k=1;k<=n-1;k++){
mi=min(k+p,n);
for(i=k+1;i<=mi ;i++)
a[i][k]=a[i][k]/a[k][k];
for(i=k+1;i<=mi ;i++){
mi=min(k+q,n);
for(j=k+1;j<=mi ;j++)
a[i][j]=a[i][j]-a[i][k]*a[k][j];
}
}
for(i=1;i<=n;i+ +)
for(j=1;j<=n;j+ +){
l[i][j]=0.0;
u[i][j]=0.0;
}
/// matrix L[][] ////
for(i=1;i<=n;i+ +){
for(j=1;j<i;j++ )
l[i][j]=a[i][j];
}
//////matrix U[][] //////////////
for(i=1;i<=n;i+ +){
for(j=i;j<=n;j+ +)
u[i][j]=a[i][j];
}
///// overwriting b[] for L[][] matrix
for(i=1;i<=n;i+ +){
mx=max(1,i-p);
zig=0.0;
for(j=mx;j<=i-1;j++)
zig=zig+l[i][j]*b[j];
b[i]=b[i]-zig;
}
/// solution matrix x[] for UX=Y
for(i=n;i>=1;i--){
mi=min(i+q,n);
sig=0.0;
for(j=i+1;j<=mi ;j++)
sig=sig+u[i][j]*b[j];
b[i]=(b[i]-sig)/u[i][i];
}
}
int min(int r,int s){
int mini=0;
mini=r;
if(s<r) mini=s;
return mini ;
}
int max(int r,int s){
int mxi=0;
mxi=r;
if(s>r) mxi=s;
return mxi;
}
if you know how can I debug it just let me know.
Thanks
6
