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;
}