I have attached the solver part of my code.
you can see its error.
matrix size is 10000x1000.
#include<math.h>
#include<fstream.h>
#include<iostream.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
