Hey,

I am writing a program that takes the size of an nxn matrix [A], randomly creates that matrix, as well as an nx1 matrix [S], multiplies them together to create [A]*[S]=.
Then, using Gaussian Elimination, I use matrix [A] and to find [x] such that [A]*[x]=.

I can see that I am not going about my algorithm the right way for this, and if anyone has some insight on what may be going wrong, I would appreciate it. I can tell based on the fact that [x] and [S] are not even a smidge close to each other.

I'm going to copy the entire code below so that it can be run.

Thanks :)

#include <stdlib.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <new>
#include <time.h>
using namespace std;

//aCtr and mCtr are used for analysis assignment
int aCtr = 0;
int mCtr = 0;


//double** make2DMatrix(int);        
//double* make1DMatrix(int);
//double randomDouble();
//void fixDiagonals(double**, int);
//void print2DMatrix(double**, int);
//void print1DMatrix(double*, int);

double** make2DMatrix(int mat_size)             
{
    double** m; 
   
    m = new double* [mat_size];
   
    for (int i = 0; i < mat_size; i++ )
    {
        m[i] = new double[mat_size];
    }
    return m;
} 

double* make1DMatrix(int mat_size)
{
    double* m;
    m = new double[mat_size];
    return m;
}

double randomDouble()
{
    int randInt = rand() % 10001;
    double randDouble = ((double)randInt) / 100.00;
    return randDouble;    
}

void fixDiagonals(double** &mat, int mat_size)
{
    double temp;
    for (int i = 0; i < mat_size; i++)
    {
        temp = 0;
        for (int j = 0; j < mat_size; j++)
        {
            if (i == j)
            {
                mat[i][j] = 0;
            }
           
            temp += 2 * mat[i][j];
        }
        mat[i][i] = temp;
    }
}
void print2DMatrix(double** &mat, int mat_size)
{
    for(int i = 0; i  < mat_size; i++)
    {
        for(int j = 0; j < mat_size; j++)
        {
            cout.precision(2);
            cout.setf(ios::fixed);
            cout << setw(8);
            cout << mat[i][j];
        }
        cout << "\n";
    }
    cout << endl;
    
}
void print1DMatrix(double* &mat, int mat_size)
{
    for(int i = 0; i < mat_size; i++)
    {
        cout.precision(2);
        cout.setf(ios::fixed);
        cout << setw(8);
        cout << mat[i] << "\n";        
    }
    cout << endl;
}

int main()
{
    double** matrixA;
    double* matrixB;
    double x, sum;
    int size;
    
    /* initialize random seed: */
    srand(time(NULL));
    
    int n;
    double diff;
    clock_t start, stop;
    
    cout << "Input the size of Matrix A, which will randomly generate an n x n matrix: ";
    cin >> size;      
    
    //Starts the clock           
    start = clock()+CLK_TCK;  
                          
    //Creating matrix A, given the inputted size
    matrixA = make2DMatrix(size);        
    for(int i = 0; i  < size; i++)
    {
        for(int j = 0; j < size; j++)
        {
            matrixA[i][j]= randomDouble();
        }
    }
    fixDiagonals(matrixA, size);
    print2DMatrix(matrixA, size);
            
    //Randomly generates the solution matrix [s]
    cout<<"\n\nRandomly generated solution matrix [s]:\n";
    double * matrixS = make1DMatrix(size); 
    for (int i = 0; i < size; i++)
    {
        matrixS[i] = randomDouble();
    }
    print1DMatrix(matrixS, size);
    
    /*
    Creates matrix [b] by multiplying [A] and [s]
    such that [A]*[s]=[b]    
    */
    matrixB = make1DMatrix(size);
    for (int i = 0; i < size; i++)
    {
        for (int j = 0; j < size; j++)
        {
            matrixB[i] += matrixA[i][j] * matrixS[i];
        } 
    }
    
    
   //Create a lower-triangul
    for (int j = 0; j < size; j++)
    {
        if (matrixA[j][j] == 0)
        {
            //The pivot point is in the right place.
            break;
        }
        else
        {
            double scale = 1.00 / matrixA[j][j];
            for(int m = j; m < size; m++)
            {
                matrixA[j][m] *= scale;
            }
            matrixB[j] *= scale;
            
            for (int i = j + 1; i < size; i++)
            {
                scale = matrixA[j][j] / matrixA[i][j];
                for (int k = 0; k < size; k++)
                {
                    matrixA[i][k] *= scale;
                    mCtr++;
                    matrixA[i][k] -= matrixA[j][k];
                    aCtr++;
                }
                matrixS[i] *= scale;
                mCtr++;
                matrixS[i] -= matrixS[j];
                aCtr++;
            } 
        }
    }
    print2DMatrix(matrixA, size);
    
    
    //*******************BACK SUBSTITUTION TO SOLVE FOR B****************** 
    matrixB[size - 1] = matrixB[size - 1] / matrixA[size -1 ][size - 1];
    for (int i = size - 1; i >= 0; i--)
    {
        double sum = 0.00;
        for (int j = i + 1; j < size; j++)
        {
            sum += matrixA[i][j] * matrixB[j];
        }
        matrixB[i] = (matrixB[i] - sum) / matrixA[i][i];
    }
    
    
    cout <<"\n\n\nSolution matrix [x]:\n\n";
    print1DMatrix(matrixB, size);
    system("pause");
   
    //Stop time and calculate total running time.
    stop = clock()+CLK_TCK;                     
    double totalTime = ((double)(stop - start)) / 1000;           
    cout << "\n\nGaussian Elimination total time: "<< totalTime <<" seconds." << endl;
  
    int complexity = aCtr + mCtr;
    cout << "Computational Complexity: " << complexity << endl; 

    
    system("pause");
    return 0;
}

Your creation of the matrix is not correct. First, when making a new matrix, the elements of it are not automatically set to zero, you have to do that. Second, at line 145, I think you have the wrong index for S (it should be j, not i).

At line 153, you cannot really compare a floating-point value to 0 (because they rarely will be exactly zero). You should check that the absolute value is less than some small number (epsilon, like 1E-10 or something similar).

At line 169, at this point, the value of MatrixA[j][j] should be 1.0 (it was just normalized) so there is no point in using the value in the matrix, just use 1.0.

At line 177 and 179, you should be using matrix B, not matrix S (remember, the matrix S is just for testing the result, it should not be part of your Gaussian elimination).

From line 169 to 180, I think you should avoid using a division by MatrixM[j] for the "scale" variable. This is because that number could very easily be zero, in which case, it should still work. So, instead of dividing row i by that number and subtracting row j from it, you should subtract "row j multiplied by the inverse scale" from row i.

For the rest, it seems alright.

Thank you :)

The only problem I'm having now is more of a cosmetic issue. When I output my matrices, even though it outputs 0.00, I think in some cases it's not exactly 0. Because, some of them output as -0.00. Is there a way to fix this issue? Should I worry about actually making those values equal to 0?

>>Is there a way to fix this issue?
You could just check if the absolute value is under epsilon and then output zero instead (as an integer). That's just for cosmetics, but I don't think anyone would mind (seeing -0.00 is very typical, nobody cares).

>>Should I worry about actually making those values equal to 0?
You should not worry about making the values equal to zero. You can display them as exactly zero as I said above, but setting the actual value to zero is meaningless. One thing that you can do, as part of the Gaussian elimination, is set the values that are cancelled to zero instead of the result of a subtraction with itself (that is also slightly more efficient... not that you should ever worry about efficiency with a Gaussian elimination since it's already (almost) the worst algorithm for solving a linear system).

This article has been dead for over six months. Start a new discussion instead.