#include<stdio.h>
#include<math.h>
#include<stdlib.h>
const int N = 3;
void LUBKSB(double b[], double a[N][N], int N, int *indx)
{
int i, ii, ip, j;
double sum;
ii = 0;
for(i=0;i<N;i++)
{
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if (ii)
{
for(j = ii;j<i-1;j++)
{
sum = sum - a[i][j] * b[j];
}
}
else if(sum)
{
ii = i;
}
b[i] = sum;
}
for(i=N-1;i>=0;i--)
{
sum = b[i];
for (j = i; j<N;j++)
{
sum = sum - a[i][j] * b[j];
}
b[i] = sum/a[i][i];
}
for (i=0;i<N;i++)
{
printf("b[%d]: %lf \n",i,b[i]);
}
}
void ludecmp(double a[][3], int N)
{
int i, imax, j, k;
double big, dum, sum, temp, d;
double *vv = (double *) malloc(N * sizeof(double));
int *indx = (int *) malloc(N * sizeof(double));
double TINY = 0.000000001;
double b[3] = {2*M_PI,5*M_PI,-8*M_PI};
d = 1.0;
for(i=0;i<N;i++)
{
big = 0.0;
for(j=0;j<N;j++)
{
temp = fabs(a[i][j]);
if (temp > big)
{
big = temp;
}
}
if (big == 0.0)
{
printf("Singular matrix\n");
exit(1);
}
vv[i] = 1.0/big;
}
for(j=0;j<N;j++)
{
for(i=0;i<j-1;i++)
{
sum = a[i][j];
for(int k=0;k<i-1;k++)
{
sum = sum - (a[i][k] * a[k][j]);
}
a[i][j] = sum;
}
big = 0.0;
for(i=j;i<N;i++)
{
sum = a[i][j];
for(k=0;k<j-1;k++)
{
sum = sum - a[i][k] * a[k][j];
}
a[i][j] =sum;
dum = vv[i] * fabs(a[i][j]);
if(dum >= big)
{
big = dum;
imax = i;
}
}
if(j != imax)
{
for(k=0;k<N;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)
{
a[j][j] = TINY;
}
if (j != N)
{
dum = 1.0/a[j][j];
for(i = j; i<N; i++)
{
a[i][j] = a[i][j] * dum;
}
}
}
LUBKSB(b,a,N,indx);
free(vv);
free(indx);
}
int main()
{
int N, i, j;
N = 3;
double a[3][3] = { 1, 2, -1, 6, -5, 4, -9, 8, -7};
ludecmp(a,N);
}
I am using these algorithms to find LU decomposition of matrix and trying to find solution A.x = b
Given a N ×N matrix A denoted as {a}N,Ni,j=1, the routine replaces it by the LU decomposition of a rowwise permutation of itself. “a” and “N” are input. “a” is also output, modified to apply the LU decomposition; {indxi}N i=1 is an output vector that records the row permutation effected by the partial pivoting; “d” is output and adopts ±1 depending on whether the number of row interchanges was even or odd. This routine is used in combination with algorithm 2 to solve linear equations or invert a matrix.
Solves the set of N linear equations A . x = b. Matrix {a} N,N i,j=1 is actually the LU decomposition of the original matrix A, obtained from algorithm 1. Vector {indxi} N i=1 is input as the permutation vector returned by algorithm 1. Vector {bi} N i=1 is input as the righthand side vector B but returns with the solution vector X. Inputs {a} N,N i,j=1, N, and {indxi} N i=1 are not modified in this algorithm.