Peano 4
Loading...
Searching...
No Matches
LUDecomposition.cpph
Go to the documentation of this file.
1template<int Rows, typename Scalar>
3 const Matrix<Rows,Rows,Scalar>& R,
4 const Vector<Rows,Scalar>& f
5) {
7
8 for(int k = Rows-1; k >= 0; k--){
9 x(k) = f(k);
10 for(int i = k+1; i < Rows; i++){
11 x(k) -= R(k, i)*x(i);
12 }
13 x(k) = x(k)/R(k, k);
14 }
15
16 return x;
17}
18
19
20template<int Rows, typename Scalar>
23 Vector<Rows,int>& pivots
24) {
25 for (int k=0; k < Rows; k++){
26 int maxIndex = k;
27 Scalar max = std::abs(A(k,k));
28 for (int i=k+1; i < Rows; i++){
29 Scalar current = std::abs(A(i,k));
30 if (current > max){
31 maxIndex = i;
32 max = current;
33 }
34 }
35 pivots(k) = maxIndex;
36
37 // Exchange lines
38 for (int j=0; j < Rows; j++){
39 Scalar temp = A(k,j);
40 A(k,j) = A(maxIndex,j);
41 A(maxIndex,j) = temp;
42 }
43
44 // Compute scaling elements
45 for (int i=k+1; i < Rows; i++){
46 A(i,k) /= A(k,k);
47 }
48
49 // Subtract contributions from each line
50 for (int i=k+1; i < Rows; i++){
51 for (int j=k+1; j < Rows; j++){
52 A(i,j) -= A(i,k) * A(k,j);
53 }
54 }
55 }
56}
57
58
59template<int Rows, typename Scalar>
62) {
63 for (int k=0; k < Rows; k++){
64 // Compute scaling elements
65 for (int i=k+1; i < Rows; i++){
66 A(i,k) /= A(k,k);
67 }
68
69 // Subtract contributions from each line
70 for (int i=k+1; i < Rows; i++){
71 for (int j=k+1; j < Rows; j++){
72 A(i,j) -= A(i,k) * A(k,j);
73 }
74 }
75 }
76}
77
78
79template<typename Scalar>
81 const Matrix<2,2,Scalar>& R
82) {
83 return R(0,0)*R(1,1)-R(0,1)*R(1,0);
84}
85
86#ifdef UseLapack
87template<int Rows, typename Scalar>
90){
91 // make a copy
92 [[maybe_unused]] tarch::la::Matrix<Rows,Rows,Scalar> M = R;
93
94 // aux variables
95 int returnCode;
96 int ipiv[Rows];
97
98 // step 1 - compute LU factorisation. Modify M in-place.
99 returnCode = LAPACKE_dgetrf( LAPACK_ROW_MAJOR, Rows, Rows, M.data(), Rows, ipiv );
100 // check return code.
101 assertion(returnCode==0);
102
103 // step 2 - use the result of the above to compute the inverse. Modify M in-place.
104 returnCode = LAPACKE_dgetri( LAPACK_ROW_MAJOR, Rows, M.data(), Rows, ipiv );
105 // check return code
106 assertion(returnCode==0);
107
108 return M;
109}
110
111#else
112
113template<typename Scalar>
115 const Matrix<2,2,Scalar>& R
116) {
117 Matrix<2,2,Scalar> result;
118
119 assertion(det(R)!=0.0);
120 const double scaling = 1.0/det(R);
121
122 result(0,0) = R(1,1) * scaling;
123 result(0,1) = -R(0,1) * scaling;
124 result(1,0) = -R(1,0) * scaling;
125 result(1,1) = R(0,0) * scaling;
126
127 return result;
128}
129
130
131
132template<typename Scalar>
134 const Matrix<3,3,Scalar>& R
135) {
136 return R(0, 0) * (R(1, 1) * R(2, 2) - R(2, 1) * R(1, 2)) -
137 R(0, 1) * (R(1, 0) * R(2, 2) - R(1, 2) * R(2, 0)) +
138 R(0, 2) * (R(1, 0) * R(2, 1) - R(1, 1) * R(2, 0));
139}
140
141
142template<typename Scalar>
144 const Matrix<3,3,Scalar>& R
145) {
146 Matrix<3,3,Scalar> result;
147
148 assertion(det(R)!=0.0);
149 double invdet = 1 / det(R);
150
151 result(0, 0) = (R(1, 1) * R(2, 2) - R(2, 1) * R(1, 2)) * invdet;
152 result(0, 1) = (R(0, 2) * R(2, 1) - R(0, 1) * R(2, 2)) * invdet;
153 result(0, 2) = (R(0, 1) * R(1, 2) - R(0, 2) * R(1, 1)) * invdet;
154 result(1, 0) = (R(1, 2) * R(2, 0) - R(1, 0) * R(2, 2)) * invdet;
155 result(1, 1) = (R(0, 0) * R(2, 2) - R(0, 2) * R(2, 0)) * invdet;
156 result(1, 2) = (R(1, 0) * R(0, 2) - R(0, 0) * R(1, 2)) * invdet;
157 result(2, 0) = (R(1, 0) * R(2, 1) - R(2, 0) * R(1, 1)) * invdet;
158 result(2, 1) = (R(2, 0) * R(0, 1) - R(0, 0) * R(2, 1)) * invdet;
159 result(2, 2) = (R(0, 0) * R(1, 1) - R(1, 0) * R(0, 1)) * invdet;
160
161 return result;
162}
163
164#endif
#define assertion(expr)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
Static (i.e.
Definition Matrix.h:31
double f(const tarch::la::Vector< Dimensions, double > &x)
int i
Definition makeIC.py:51
Vector< Rows, Scalar > backSubstitution(const Matrix< Rows, Rows, Scalar > &R, const Vector< Rows, Scalar > &f)
Accepts an upper triangular matrix and a rhs.
double det(const Matrix< 2, 2, Scalar > &R)
double max(double a, double b, double c)
I need the maximum of three values all the time, to I decided to write a function for this.
Definition Scalar.cpp:8
void lu(Matrix< Rows, Rows, Scalar > &A, Vector< Rows, int > &pivots)
Performs an in-situ LU-decomposition of the square matrix A.
Matrix< 2, 2, Scalar > invert(const Matrix< 2, 2, Scalar > &R)
Simple vector class.
Definition Vector.h:134