This blog explores the ways to solve Ax=B, where A is nxn matrix, x is nx1 matrix and B is nx1 matrix. On the way, I will take a detour, where I will write a code to convert A into an upper triangular matrix, also I will validate it whether A is a singular matrix.
The code is written in C++. When we have a system of linear equations, we can write it as
sum(Aij*Xj=Bi). In another blog, I will write code to find the inverse of non-singular n*n matrix. For now, We will focus on X's solution and way to convert A into upper triangular matrix.
The main utility function is upT(A, row, col, b), it takes A , its row length, column length and B matrix.
Validate(A, m, n) checks whether A has zero at pivot position. If pivot is found, then it exchanges zero pivot row with non-zero row. First for loop in upT function converts A into upper triangular matrix and then singular(A, m, n-1) checks whether A is a singular matrix. Second for loop converts upper triangular matrix into Identity matrix to find solution X's , where Xj=Ai/Bi (A being diagonal matrix)
Time complexity: O(n^3). (If you can improve it, please let me know)
#include <iostream>
using namespace std;
bool singular(float** A, int m, int n){
for(int i=0;i<m;i++){
int non_zero=0;
for(int j=0;j<n;j++){
if(A[i][j]!=0){
non_zero++;
}
}
if(non_zero<1){
return true;
}
}
return false;
}
void validate(float** a, int m, int n){
for(int i=0;i<m;i++){
for(int j=0;j<n;j++){
if(i==j and a[i][j]==0){
int found=0;
for(int k=i+1;k<m;k++){
if(a[k][j]!=0){
found=1;
for(int u=0;u<n;u++){
swap(a[i][u], a[k][u]);
}
}
if(found==1){
break;
}
}
}
}
}
}
void upT(float** A, int m, int n, float* b){
validate(A, m, n);
for(int i=0;i<m-1;i++){
int j=m-1;
while(j>i){
int k=i;
float r=A[j][i]/A[i][i];
int start=1;
while(k<n){
if(start==1){
A[j][k]=A[j][k]-(A[j][i]/A[i][i])*A[i][k];
start=0;
}else{
A[j][k]=A[j][k]-r*A[i][k];
}
k++;
}
j--;
}
}
if(singular(A, m, n-1)){
cout << "A is a singular matrix." << endl;
return;
}
cout << "Upper triangular matrix: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n-1;j++){
cout << A[i][j] << " ";
}
cout << endl;
}
for(int i=0;i<m-1;i++){
int j=i;
while(j<m-1){
int k=j+1;
int start=1;
int r=A[i][j+1]/A[j+1][j+1];
while(k<n){
if(start==1){
A[i][k]=A[i][k]-A[i][j+1]/A[j+1][j+1]*A[j+1][k];
start=0;
}else{
A[i][k]=A[i][k]-r*A[j+1][k];
}
k++;
}
j++;
}
}
cout << "diagonal matrix A: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n-1;j++){
cout << A[i][j] << " ";
}
cout << endl;
}
cout << "solution for X's: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n;j++){
if(i==j){
cout << b[i]/A[i][i] << endl;
}
}
}
}
int main(){
int m, n;
cout << "row & col of matrix: " << endl;
cin >> m >> n;
float** a;
a=new float*[m];
for(int i=0;i<n;i++){
a[i]=new float[n];
}
cout << "enter matrix elements: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n;j++){
cin >> a[i][j];
}
}
cout << "enter b elements: " << endl;
float b[m];
for(int i=0;i<m;i++){
cin >> b[i];
}
float** A;
int col=n+1;
A=new float*[m];
for(int i=0;i<col;i++){
A[i]=new float[col];
}
for(int i=0;i<m;i++){
for(int j=0;j<col;j++){
A[i][j]=(j==n)?(b[i]):a[i][j];
}
}
upT(A, m, col, b);
return 0;
}
Feel free to comment. If any error is found, please let me know
Reference : Matrix
The code is written in C++. When we have a system of linear equations, we can write it as
sum(Aij*Xj=Bi). In another blog, I will write code to find the inverse of non-singular n*n matrix. For now, We will focus on X's solution and way to convert A into upper triangular matrix.
The main utility function is upT(A, row, col, b), it takes A , its row length, column length and B matrix.
Validate(A, m, n) checks whether A has zero at pivot position. If pivot is found, then it exchanges zero pivot row with non-zero row. First for loop in upT function converts A into upper triangular matrix and then singular(A, m, n-1) checks whether A is a singular matrix. Second for loop converts upper triangular matrix into Identity matrix to find solution X's , where Xj=Ai/Bi (A being diagonal matrix)
Time complexity: O(n^3). (If you can improve it, please let me know)
#include <iostream>
using namespace std;
bool singular(float** A, int m, int n){
for(int i=0;i<m;i++){
int non_zero=0;
for(int j=0;j<n;j++){
if(A[i][j]!=0){
non_zero++;
}
}
if(non_zero<1){
return true;
}
}
return false;
}
void validate(float** a, int m, int n){
for(int i=0;i<m;i++){
for(int j=0;j<n;j++){
if(i==j and a[i][j]==0){
int found=0;
for(int k=i+1;k<m;k++){
if(a[k][j]!=0){
found=1;
for(int u=0;u<n;u++){
swap(a[i][u], a[k][u]);
}
}
if(found==1){
break;
}
}
}
}
}
}
void upT(float** A, int m, int n, float* b){
validate(A, m, n);
for(int i=0;i<m-1;i++){
int j=m-1;
while(j>i){
int k=i;
float r=A[j][i]/A[i][i];
int start=1;
while(k<n){
if(start==1){
A[j][k]=A[j][k]-(A[j][i]/A[i][i])*A[i][k];
start=0;
}else{
A[j][k]=A[j][k]-r*A[i][k];
}
k++;
}
j--;
}
}
if(singular(A, m, n-1)){
cout << "A is a singular matrix." << endl;
return;
}
cout << "Upper triangular matrix: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n-1;j++){
cout << A[i][j] << " ";
}
cout << endl;
}
for(int i=0;i<m-1;i++){
int j=i;
while(j<m-1){
int k=j+1;
int start=1;
int r=A[i][j+1]/A[j+1][j+1];
while(k<n){
if(start==1){
A[i][k]=A[i][k]-A[i][j+1]/A[j+1][j+1]*A[j+1][k];
start=0;
}else{
A[i][k]=A[i][k]-r*A[j+1][k];
}
k++;
}
j++;
}
}
cout << "diagonal matrix A: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n-1;j++){
cout << A[i][j] << " ";
}
cout << endl;
}
cout << "solution for X's: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n;j++){
if(i==j){
cout << b[i]/A[i][i] << endl;
}
}
}
}
int main(){
int m, n;
cout << "row & col of matrix: " << endl;
cin >> m >> n;
float** a;
a=new float*[m];
for(int i=0;i<n;i++){
a[i]=new float[n];
}
cout << "enter matrix elements: " << endl;
for(int i=0;i<m;i++){
for(int j=0;j<n;j++){
cin >> a[i][j];
}
}
cout << "enter b elements: " << endl;
float b[m];
for(int i=0;i<m;i++){
cin >> b[i];
}
float** A;
int col=n+1;
A=new float*[m];
for(int i=0;i<col;i++){
A[i]=new float[col];
}
for(int i=0;i<m;i++){
for(int j=0;j<col;j++){
A[i][j]=(j==n)?(b[i]):a[i][j];
}
}
upT(A, m, col, b);
return 0;
}
Feel free to comment. If any error is found, please let me know
Reference : Matrix
No comments:
Post a Comment