-
Notifications
You must be signed in to change notification settings - Fork 0
/
LinearSystem.cpp
138 lines (125 loc) · 3.83 KB
/
LinearSystem.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
134
135
136
137
138
#include"Master.hpp"
#include<cassert>
#include<fstream>
#include<iostream>
using namespace std;
// copy matrix and vector so that by gaussian elimination original matrix and vector remains unchanged
LinearSystem::LinearSystem(const Matrix& A, const Vector& b )
{
//Check for size compatibility of matrix and vector
int local_size= A.GetNumberOfRows();
assert(A.GetNumberOfCols() == local_size);
assert(b.GetSize() == local_size);
//Set variables for linear system
mSize=local_size;
mpA = new Matrix(A);
mpb = new Vector(b);
}
// Destructor frees memory
LinearSystem::~LinearSystem()
{
delete mpA;
delete mpb;
}
//Solve the linear system using Gaussian elimination
//This will change the matrix mpA
void LinearSystem::Mumps_input_gen()
{
//Introducing pointer for easy access to Matrix and vector
Matrix rA = *mpA;
Vector rb = *mpb;
//Printing the matrix and rhs in form of irn and jcn form
std::ofstream file_1("irn.txt");
std::ofstream file_2("jcn.txt");
std::ofstream file_3("rhs.txt");
std::ofstream file_4("a.txt");
//std::cout<<"HI\n";
if (!file_1.is_open() && !file_2.is_open() && !file_3.is_open() && !file_4.is_open())
{
std::cout<<"File cannot be opened while creating irn.txt etc.!";
}
else
{
assert(file_1.is_open());
assert(file_2.is_open());
assert(file_3.is_open());
assert(file_4.is_open());
int N = rA.GetNumberOfRows();
// std::cout<<"\n"<<N<<"\n";
for(int i = 1; i < N+1; i++)
{
for(int j = 1; j < N+1; j++)
{
if(rA(i,j) != 0)
{
file_1<<i<<"\n";
file_2<<j<<"\n";
double y = rA(i,j);
file_4<<y<<"\n";
}
}
double z = rb[i-1];
file_3<<z<<"\n";
}
//file_1.flush();
file_1.close();
//file_2.flush();
file_2.close();
//file_3.flush();
file_3.close();
//file_4.flush();
file_4.close();
}
//std::cout<<"Ouptput done!\n";
// //forward sweep of gaussian elimination (conversion to upper triangular form)
// for (int k=0; k < mSize-1; k++)
// {
// //see if pivoting is necessary
// double max=0.0;
// int row = -1;
// for(int i = k; i < mSize; i++)
// {
// if (fabs(rA(i+1,k+1)) > max)
// {
// row = i;
// }
// }
// assert (row > 0);
// // pivot if necessary
// if (row != k)
// {
// //swap matrix k+1 with row+1
// for (int i=0; i < mSize; i++)
// {
// double temp = rA(k+1,i+1);
// rA(k+1,i+1) = rA(row+1,i+1);
// rA(row+1,i+1) = temp;
// }
// //swap vector enteries k+1 with row+1
// double temp = rb(k+1);
// rb(k+1) = rb(row+1);
// rb(row+1) = temp;
// }
// //create 0 in the lower part of column k
// for (int i=k+1; i < mSize; i++)
// {
// m(i+1) = rA(i+1,k+1)/rA(k+1,k+1);
// for (int j=k; j < mSize; j++)
// {
// rA(i+1,j+1) -= rA(k+1,j+1)*m(i+1);
// }
// rb(i+1) -= rb(k+1)*m(i+1);
// }
// }
// //back substitution
// for (int i=mSize-1; i > -1; i--)
// {
// solution(i+1) = rb(i+1);
// for (int j=i+1; j < mSize; j++ )
// {
// solution(i+1) -= rA(i+1,j+1)*solution(j+1);
// }
// solution (i+1) /= rA(i+1,i+1);
// }
// return solution;
}