LU Decomposition not working?

[Miles Hikari]

So for one of my projects I need to write an algorithm for LU Decomposition, the only issue is that part way through the explanation the professor got himself confused and swapped the Lower and Upper matrices by accident mid lecture, no one caught it until afterwards. Anyway after piecing together what the professor was trying to say I got to writing the code. One of the key things about this system is its supposed to take a matrix, and split it into two triangular matrices that are supposed to be able to be multiplied together to once again become the original matrix... But that's quite honestly not happening. Being multiplied together gives me something totally different... Also on a side note, and I know this is poor programming but all of my calculations are being kept in a Matrix header file (as opposed to a header for declarations and cpp) The prof gave us a simple matrix building header and I've been building on it since adding more functions, I just somewhere along the line forgot to split it....

Anyway I even went back and googled this particular code in c++ (I believe it's crout's method only keeping the 1's on the lower) to see what I did wrong and according to what I found everything should be working perfect. I'm just just sure where I messed up. If anyone could offer a hand I'd appreciate it, because my prof is a bit difficult to track down most days. code in next post.

[Miles Hikari]

matrix.h

#include <random>
#include <iomanip>
#include <cmath>
#include <ctime>
#include <iostream>

using namespace std;
class mat
{
public:
    mat()
    {
        A=NULL;
    }
    ~mat() {}
//////////////////////////////////////////////////////////////////////
    mat(int r, int c)
    {
        row=r;
        col=c;         ///Build an empty matrix of row r and column c
        A=new double *[row];
        for(int i=0; i<row; i++)
        {
            A[i]=new double[col];
        }
        for(int i=0; i<row; i++)
        {
            for(int j=0; j<col; j++)
            {
                A[i][j]=0.0;
            }
        }
    }
/////////////////////////////////////////////////////////////////////
    int getrow()
    {
        return row;
    }
    int getcol()
    {
        return col;   ///Retrieve the needed info
    }
    void setValue(int r, int c, double val)
    {
        A[r][c] = val;
    }

    double getValue(int r, int c)
    {
        return A[r][c];
    }
/////////////////////////////////////////////////////////////////////
    void display()
    {
        for(int i=0; i<row; i++)
        {
            for(int j=0; j<col; j++)
            {
                cout<<setw(8)<<A[i][j];
            }
            cout<<endl;
        }
        cout<<endl;
    }
/////////////////////////////////////////////////////////////////////
void random_mat(double m)   ///Fill the matrix with random values
    {
        double sum=0.0;
        double r;
        default_random_engine gen(time(NULL));
        uniform_real_distribution<double> rdist(-1*m, m);
        for(int i=0; i<row; i++)
        {
            for(int j=0; j<col; j++)
            {
                r=rdist(gen);
                A[i][j]=r;
                sum = sum+abs(r);
            }
            A[i][i]=sum;
            sum=0.0;
        }
    }
////////////////////////////////////////////////////////////////////////

mat matMul(mat &a, mat &b)
    {

        int aRow = a.getrow();
        int aCol = a.getcol();
        int bRow = b.getrow();
        int bCol = b.getcol();

        if (aRow == 1 && aCol == 1)             /// If matrix a has only one value treat it like a scaler to b
        {

            mat c(bRow,bCol);
            for(int i=0; i<bRow; i++)
            {
                for(int j=0; j<bCol; j++)
                {
                    c.A[i][j] = b.A[i][j] * a.A[0][0];
                }
            }
            return c;
        }
////////////////////////////////////////////////////////////////////////////////////////////////////////////
        else if (bRow == 1 && bCol == 1)     ///if matrix b only has one value, treat it like a scaler to a
        {
            mat c(aRow,aCol);
            for(int i=0; i<aRow; i++)
            {
                for(int j=0; j<aCol; j++)
                {
                    c.A[i][j] = a.A[i][j] * b.A[0][0];
                }
            }
            return c;
        }
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
        else if (aRow == 1 && bCol == 1 && aCol == bRow)///if a row is multiplied by a column, get a single number
        {
            mat c(1,1);
            double sum = 0;
            for(int i=0; i<aCol; i++)
            {

                sum = sum + a.A[0][i] * b.A[i][0];
            }
            c.A[0][0] = sum;
            return c;
        }
///////////////////////////////////////////////////////////////////////////////////////////////////////
        else if (aCol == 1 && bRow == 1 && bCol == aRow) /// a column multiplied by a row
        {
            mat c(aRow,bCol);
            for(int i=0; i<aRow; i++)
            {
                for(int j=0; j<bCol; j++)
                {
                    c.A[i][j] = b.A[0][j] * a.A[i][0];
                }
            }
            return c;
        }

/////////////////////////////////////////////////////////////////////////////////////////////////////////
        else if (aCol == bRow)  ///Matrix multiplied by another matrix (m*n)*(n*p)
        {
            double sum;
            mat c(aRow,bCol);
            for(int i=0; i<aRow; i++)
            {

                for(int j=0; j<bCol; j++)
                {
                    sum = 0;
                    for(int k=0; k<bRow; k++)
                    {

                        sum = sum + (a.A[i][k] * b.A[k][j]);
                    }
                    c.A[i][j] = sum;
                }
            }
            return c;
        }
    }

    void buildU(mat ori) ///Build the Upper triangle base
    {

        for (int i=0; i<ori.getcol(); i++)
        {
            this->A[0][i] = ori.A[0][i];
        }

        //this->mat() = upper;
    }

    void buildL(mat ori) ///Build the lower triangle base, 1s on diagonal
    {
        this->A[0][0] = 1;
        for (int i=1; i<ori.getrow(); i++)
        {
            this->A[i][0] = (ori.A[i][0] / ori.A[0][0]);
            this->A[i][i] = 1;
        }

    }

    void factorize(mat upper, mat lower) /// Use crouchs method to flesh out the triangular matrices
    {
        for (int i = 1; i < this->getrow(); i++)
        {
            for(int k = i; k < this->getrow(); k++) ///Begin filling out a row of the upper triangle
            {
                double sum = 0;
                for(int j = 0; j < k; j++)
                {
                    sum = sum + (lower.getValue(j,i) * upper.getValue(i,j));
                }
                sum = this->getValue(i,k) - sum;
                upper.setValue(i,k,sum);
                //upper.display();

            }
//////////////////////////////////////////////////////////////////////////////////////
            for(int k = i; k < getrow(); k++) ///Fill in a column of the lower triangle
            {
                double sum = 0;
                for(int j = 0; j < k; j++)
                {
                    sum = sum + (lower.getValue(j,i) * upper.getValue(i,j));
                }
                if (this->getValue(i,i) == 0)
                {
                    cout << "0 on the diagonal" << endl;
                    break;
                }
                sum = (this->getValue(k,i) - sum) / upper.getValue(i,i);
                lower.setValue(k,i,sum);
                //lower.display();
            }

        }

    }
////////////////////////////////////////////////////////////////////////////////////////////////




private:
    int row;
    int col;
    double ** A;
};

main.cpp (the file i'm using to make sure it works

 #include "matrix.h"


int main()
{
mat A(3,3);
mat B(3,3);
mat C(3,3);
///A.random_mat(5);
A.setValue(0,0,1);
A.setValue(0,1,-2);
A.setValue(0,2,3);
A.setValue(1,0,2);
A.setValue(1,1,-5);
A.setValue(1,2,12);
A.setValue(2,0,0);
A.setValue(2,1,2);
A.setValue(2,2,-10);
A.display();

B.buildU(A);

B.display();

C.buildL(A);

C.display();

A.factorize(B, C);

B.display();
C.display();

A = B.matMul(B,C);

A.display();
}

Also yes I know that I could fix matMul to effect the called matrix directly, I had to build that function before i was properly introduced to this->, and havent had the time to go back and fix it.

[lastchance]

*** REQUEST FOR FURTHER FORUM HELP AT THE END OF THIS POST ***


In matrix.h:

Line 199 AND line 212 should have

j < i

and not

j < k

Line 209 should have

int k = i + 1

and not

int k = i

(you have already set the diagonal element)


Line 201 should be
sum = sum + (lower.getValue(i,j) * upper.getValue(j,k));


Line 214 should be
sum = sum + (lower.getValue(k,j) * upper.getValue(j,i));



Line 192 requires references if you want to return some values
void factorize(mat &upper, mat &lower)



In main.cpp the multiplication order should be LU, and you have UL.
Line 34 in main.cpp should be
A = B.matMul(C,B);



OK - that solves your immediate problem (at least, it did on my PC).
However, there are a huge number of problems with your code.

First, on the LU algorithm. It is hugely confusing to have separate buildU and buildL routines to set first row and column and the unit diagonal elements. Just put the appropriate lines of code in the factorize() method.

Second, the matrix multiply. Why do you need all these cases (which won't actually work anyway for a non-square matrix)? You only need a single case: there is nothing wrong with having a loop that only actually runs once.


Then the C++ ...
- all this code in the header file ... ?
- you don't need

this->

all the time: if it's a member function then it has access to all the member data anyway.
- you are using new ... without delete; this is going to cause significant memory leaks; your destructor needs to clear these pointers.


*** FURTHER FORUM HELP REQUIRED ***
I'm going to call for help from the more C++-knowledgeable on this forum as to whether you need to overload the = operator for your (to-be-corrected) line
A = B.matMul(C,B);
since both A and the return matrix from the function will have dynamically-allocated memory.

[48T7M4Gy]

The overloaded = operator etc.

A 'simple' explanation is at http://courses.cms.caltech.edu/cs11/material/cpp/donnie/cpp-ops.html

Ideally operator* overload is a good idea so that multiplying two matrices to give a reusltant C is simply C = A*B