32 #ifndef GAUSSIAN_ELIMINATION_H
33 #define GAUSSIAN_ELIMINATION_H
44 template<
int N,
typename Precision>
51 for (
int i=0; i<size; ++i) {
53 Precision maxval =
abs(A[i][i]);
55 for (
int ii=i+1; ii<size; ++ii) {
56 double v =
abs(A[ii][i]);
62 Precision pivot = A[argmax][i];
64 Precision inv_pivot =
static_cast<Precision
>(1)/pivot;
66 for (
int j=i; j<size; ++j)
67 swap(A[i][j], A[argmax][j]);
68 swap(b[i], b[argmax]);
71 for (
int j=i+1; j<size; ++j)
75 for (
int u=i+1; u<size; ++u) {
76 double factor = A[u][i];
78 for (
int j=i+1; j<size; ++j)
79 A[u][j] -= factor * A[i][j];
80 b[u] -= factor * b[i];
85 for (
int i=size-1; i>=0; --i) {
87 for (
int j=i+1; j<size; ++j)
88 x[i] -= A[i][j] * x[j];
95 template<
int i,
int j,
int k>
struct Size3
97 static const int s=(i!= -1)?i:(j!=-1?j:k);
106 template<
int R1,
int C1,
int R2,
int C2,
typename Precision>
113 int size=A.num_rows();
115 for (
int i=0; i<size; ++i) {
117 Precision maxval =
abs(A[i][i]);
119 for (
int ii=i+1; ii<size; ++ii) {
120 double v =
abs(A[ii][i]);
126 Precision pivot = A[argmax][i];
128 Precision inv_pivot =
static_cast<Precision
>(1)/pivot;
130 for (
int j=i; j<size; ++j)
131 swap(A[i][j], A[argmax][j]);
133 for(
int j=0; j < b.num_cols(); j++)
134 swap(b[i][j], b[argmax][j]);
137 for (
int j=i+1; j<size; ++j)
138 A[i][j] *= inv_pivot;
141 for (
int u=i+1; u<size; ++u) {
142 double factor = A[u][i];
144 for (
int j=i+1; j<size; ++j)
145 A[u][j] -= factor * A[i][j];
146 b[u] -= factor * b[i];
151 for (
int i=size-1; i>=0; --i) {
152 for(
int k=0; k <b.num_cols(); k++)
155 for (
int j=i+1; j<size; ++j)
156 x[i][k] -= A[i][j] * x[j][k];
Everything lives inside this namespace.
static void test(int s1, int s2)
Vector< N, Precision > gaussian_elimination(Matrix< N, N, Precision > A, Vector< N, Precision > b)