-
Notifications
You must be signed in to change notification settings - Fork 0
/
sequential.cpp
133 lines (112 loc) · 3 KB
/
sequential.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include<bits/stdc++.h>
using namespace std;
void random_initialization(vector<vector<float>> &matrix)
{
cout<<"\nRandomly initializing input matrix A ....."<<endl;
int size = matrix.size();
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
matrix[i][j] = (rand()%10)+1;
}
}
// Making the matrix diagonal dominant, so that it is invertible
int diagonal = 0;
float sum = 0;
for (int i = 0; i < size; i++)
{
for (int j = 0; j < size; j++)
{
sum += abs(matrix[i][j]); // sum of all values in column j
}
sum -= abs(matrix[i][diagonal]); //Removing the diagonal value from the sum
// For diagonal dominance, diagonal > sum of all other elements
matrix[i][diagonal] = sum + (rand()%5)+1;
diagonal++;
sum = 0;
}
}
void print_matrix(vector<vector<float>> matrix)
{
int size = matrix.size();
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
cout <<left<<setw(9)<<setprecision(3)<<matrix[i][j]<<left<<setw(9);
}
cout<<endl;
}
cout<<endl;
}
void LU_Decomposition(vector<vector<float>> &A,vector<vector<float>> &L,vector<vector<float>> &U)
{
int size = A.size();
// Finding lower triangular matrix L
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
if (j < i) // set upper diagonal elements to 0
{
L[j][i] = 0;
continue;
}
L[j][i] = A[j][i];
for (int k = 0; k < i; k++)
{
//deduct from the current L cell the value of these 2 values multiplied
L[j][i] = L[j][i] - L[j][k] * U[k][i];
}
}
}
// Finding upper triangular matrix U
for(int i=0;i<size;i++)
{
for(int j=0;j<size;j++)
{
if (i > j)
{
U[i][j] = 0;
continue;
}
//if they're equal, set u's current index to 1
if (j == i)
{
U[i][j] = 1;
continue;
}
//otherwise, do some math to get the right value
U[i][j] = A[i][j] / L[i][i];
for (int k = 0; k < i; k++)
{
U[i][j] = U[i][j] - ((L[i][k] * U[k][j]) / L[i][i]);
}
}
}
}
int main()
{
cout<<"\n----- SEQUENTIAL IMPLEMENTATION OF LU DECOMPOSITION -----"<<endl;
double runtime = 0.0;
int size;
cout<<"\nEnter the size of the matrix : ";
cin>>size;
vector<vector<float>> A(size,vector<float>(size));
vector<vector<float>> L(size,vector<float>(size));
vector<vector<float>> U(size,vector<float>(size));
random_initialization(A);
cout<<"\nMatrix A : \n"<<endl;
print_matrix(A);
clock_t startTime = clock(); // starting the clock time
LU_Decomposition(A, L, U);
clock_t endTime = clock(); // ending the clock time
cout<<"\nLower triangular matrix (L Matrix) : \n"<<endl;
print_matrix(L);
cout<<"\nUpper triangular matrix (U Matrix) : \n"<<endl;
print_matrix(U);
runtime = (double)(endTime-startTime)/(double)CLOCKS_PER_SEC;
// cout <<"\nTime elapsed : "<<runtime<<" seconds\n"<<endl;
return 0;
}