0
#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.

  • 1
    The compiler is giving 5 warnings from this code. Two of them are about functions which do not return a value, but should. The next step is to enable full warnings (if you aren't getting them) and fix them all. – Weather Vane Jan 31 '22 at 07:27
  • In fact the value from those functions is not used. Two more are about conversion of `double` to `int` (which might or might not matter), another is about `exit;` which should be `exit(something)`. Also there is a mistake in `#define TINY 0.000000001;` which should be `#define TINY 0.000000001` but hasn't caused an error. – Weather Vane Jan 31 '22 at 07:31
  • @WeatherVane I changed the function type to void, it still doesn't print the matrix – Shubhendu Das Jan 31 '22 at 07:34
  • I see - my compilation crashes, for a matrix size 2. – Weather Vane Jan 31 '22 at 07:35
  • Your loops are incorrect: `for(i=0;i<=N;i++)` should be `for(i=0;i=0;i--)` which should be `for(i=N-1;i>=0;i--)` – Weather Vane Jan 31 '22 at 07:38
  • If you `#include `, why don't you use `M_PI`? Instead, you `#define` it, with less precision than it's `double` value representation. – Luis Colorado Feb 02 '22 at 11:57
  • Can you edit your code and use meaningful variable names? starting with `a`, `b`, ... `z`, `aa`, ... is not the way. This way it is impossible to follow your algorithm, even for you. – Luis Colorado Feb 02 '22 at 12:00
  • @LuisColorado to be fair, [`M_PI` isn't standard](https://stackoverflow.com/q/26065359/995714). Only since C++20 you can use `std::numbers::pi` [Does C++11, 14, 17 or 20 introduce a standard constant for pi?](https://stackoverflow.com/q/49778240/995714) but there's still no pi in C – phuclv Feb 03 '22 at 07:36
  • `M_PI` is Included in X/Open Systems Interface standard, which is an extension of ANSI-C. It's also part of POSIX.1-2017 as you can see [here](https://pubs.opengroup.org/onlinepubs/9699919799/basedefs/math.h.html). Of course you decide which standard you want to follow, but today almost no unix is not POSIX compliant. – Luis Colorado Feb 03 '22 at 08:52
  • I made all the changes which were suggested. I still can't get indexing right – Shubhendu Das Feb 03 '22 at 23:19
  • @WeatherVane Can you tell me which one of my indexing is still wrong please? – Shubhendu Das Feb 03 '22 at 23:23

1 Answers1

0

There are a number of problems with your code:

  1. In your for-loops, i <= N should be i < N and i = N should be i = N - 1.
  2. The absolute value of a double is returned by fabs, not abs.
  3. The statement exit should be exit(1) or exit(EXIT_FALILURE).
  4. Two of your functions lack a return statement.

You should also free the memory you have allocated with the function free. When you compile a C program you should also enable all warnings.

August Karlstrom
  • 10,773
  • 7
  • 38
  • 60