Created
March 28, 2011 17:23
-
-
Save tomhermann/890877 to your computer and use it in GitHub Desktop.
Found some of my earlier college C++ work on an old hard drive... man this was bad.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
template <class T> | |
Matrix<T> gaussian_solver<T>::operator()(const Matrix<T>& A, const Vector<T>& b) { | |
if(!A.is_square()) { | |
throw CSignal("Gaussian solver can only be applied to square matrices"); | |
} | |
Matrix<T> system = A; // make a copy because we're going to mess it up | |
Vector<T> row_temp; | |
system = system.augment(b); // Augment vector b, now we have a system | |
const size_t m = system.getRows(); | |
const size_t n = system.getCols(); | |
size_t i=0, j=0, max_ind=0, temp=0; | |
T val, max_val; | |
while (i < m && j < n) | |
{ | |
//Find pivot in column j, starting in row i | |
max_val = system(i,j); | |
max_ind = i; | |
for(size_t k=i; k<m; k++) | |
{ | |
val = system(k, j); | |
if(tabs(val) > tabs(max_val)) | |
{ | |
max_val = val; | |
max_ind = k; | |
} | |
} | |
if(tabs(max_val) > epsilon) | |
{ | |
system.switchRow(i,max_ind); | |
// Now system(i, j) will contain the same value as max_val | |
// divide row i by max_val | |
for(size_t k=0; k < n; k++) | |
{ | |
system(i,k) = system(i,k) / max_val; | |
} | |
// Now sytem(i, j) will have the value 1 | |
for(size_t u=0; u < m; u++) | |
{ | |
if (u != i) | |
{ | |
// subtract A[u, j] * row i from row u | |
for(size_t k=i+1; k<n; k++) | |
{ | |
system(u,k) -= (system(u,j) * system(i,k)); | |
} | |
} | |
} | |
i++; | |
} | |
j++; | |
} | |
return system; | |
} | |
template <class T> | |
void gaussian_solver<T>::setprecision(double eps) | |
{ | |
if(eps < 1 && eps > 0) | |
epsilon = eps; | |
} | |
template <class T> | |
T gaussian_solver<T>::tabs(const T number) | |
{ | |
T temp = number; | |
if(temp < 0) { | |
temp = -temp; | |
} | |
return temp; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment