PandA-2024.02
LUdecomposition.c
Go to the documentation of this file.
1 // int Upper_Triangular_Solve(float *U, float *B, float x[], int n) //
3 // //
4 // Description: //
5 // This routine solves the linear equation Ux = B, where U is an n x n //
6 // upper triangular matrix. (The subdiagonal part of the matrix is //
7 // not addressed.) //
8 // The algorithm follows: //
9 // x[n-1] = B[n-1]/U[n-1][n-1], and //
10 // x[i] = [B[i] - (U[i][i+1] * x[i+1] + ... + U[i][n-1] * x[n-1])] //
11 // / U[i][i], //
12 // for i = n-2, ..., 0. //
13 // //
14 // Arguments: //
15 // float *U Pointer to the first element of the upper triangular //
16 // matrix. //
17 // float *B Pointer to the column vector, (n x 1) matrix, B. //
18 // float *x Pointer to the column vector, (n x 1) matrix, x. //
19 // int n The number of rows or columns of the matrix U. //
20 // //
21 // Return Values: //
22 // 0 Success //
23 // -1 Failure - The matrix U is singular. //
24 // //
25 // Example: //
26 // #define N //
27 // float A[N][N], B[N], x[N]; //
28 // //
29 // (your code to create matrix A and column vector B) //
30 // err = Upper_Triangular_Solve(&A[0][0], B, x, n); //
31 // if (err < 0) printf(" Matrix A is singular\n"); //
32 // else printf(" The solution is \n"); //
33 // ... //
35 // //
36 int Upper_Triangular_Solve(float *U, float B[], float x[], int n)
37 {
38  int i, k;
39 
40 // Solve the linear equation Ux = B for x, where U is an upper
41 // triangular matrix.
42 
43  for (k = n-1, U += n * (n - 1); k >= 0; U -= n, k--) {
44  if (*(U + k) == 0.0) return -1; // The matrix U is singular
45  x[k] = B[k];
46  for (i = k + 1; i < n; i++) x[k] -= x[i] * *(U + i);
47  x[k] /= *(U + k);
48  }
49 
50  return 0;
51 }
52 
54 // void Unit_Lower_Triangular_Solve(float *L, float *B, float x[], int n) //
55 // //
56 // Description: //
57 // This routine solves the linear equation Lx = B, where L is an n x n //
58 // unit lower triangular matrix. (Only the subdiagonal part of the matrix//
59 // is addressed.) The diagonal is assumed to consist of 1's and is not //
60 // addressed. //
61 // The algorithm follows: //
62 // x[0] = B[0], and //
63 // x[i] = B[i] - (L[i][0] * x[0] + ... + L[i][i-1] * x[i-1]), //
64 // for i = 1, ..., n-1. //
65 // //
66 // Arguments: //
67 // float *L Pointer to the first element of the unit lower triangular //
68 // matrix. //
69 // float *B Pointer to the column vector, (n x 1) matrix, B. //
70 // float *x Pointer to the column vector, (n x 1) matrix, x. //
71 // int n The number of rows or columns of the matrix L. //
72 // //
73 // Return Values: //
74 // void //
75 // //
76 // Example: //
77 // #define N //
78 // float A[N][N], B[N], x[N]; //
79 // //
80 // (your code to create matrix A and column vector B) //
81 // Unit_Lower_Triangular_Solve(&A[0][0], B, x, n); //
82 // printf(" The solution is \n"); //
83 // ... //
85 // //
86 void Unit_Lower_Triangular_Solve(float *L, float B[], float x[], int n)
87 {
88  int i, k;
89 
90 // Solve the linear equation Lx = B for x, where L is a unit lower
91 // triangular matrix.
92 
93  x[0] = B[0];
94  for (k = 1, L += n; k < n; L += n, k++)
95  for (i = 0, x[k] = B[k]; i < k; i++) x[k] -= x[i] * *(L + i);
96 }
97 
99 // int Doolittle_LU_Decomposition(float *A, int n) //
100 // //
101 // Description: //
102 // This routine uses Doolittle's method to decompose the n x n matrix A //
103 // into a unit lower triangular matrix L and an upper triangular matrix U //
104 // such that A = LU. //
105 // The matrices L and U replace the matrix A so that the original matrix //
106 // A is destroyed. //
107 // Note! In Doolittle's method the diagonal elements of L are 1 and are //
108 // not stored. //
109 // Note! The determinant of A is the product of the diagonal elements //
110 // of U. (det A = det L * det U = det U). //
111 // This routine is suitable for those classes of matrices which when //
112 // performing Gaussian elimination do not need to undergo partial //
113 // pivoting, e.g. positive definite symmetric matrices, diagonally //
114 // dominant band matrices, etc. //
115 // For the more general case in which partial pivoting is needed use //
116 // Doolittle_LU_Decomposition_with_Pivoting. //
117 // The LU decomposition is convenient when one needs to solve the linear //
118 // equation Ax = B for the vector x while the matrix A is fixed and the //
119 // vector B is varied. The routine for solving the linear system Ax = B //
120 // after performing the LU decomposition for A is Doolittle_LU_Solve //
121 // (see below). //
122 // //
123 // The Doolittle method is given by evaluating, in order, the following //
124 // pair of expressions for k = 0, ... , n-1: //
125 // U[k][j] = A[k][j] - (L[k][0]*U[0][j] + ... + L[k][k-1]*U[k-1][j]) //
126 // for j = k, k+1, ... , n-1 //
127 // L[i][k] = (A[i][k] - (L[i][0]*U[0][k] + . + L[i][k-1]*U[k-1][k])) //
128 // / U[k][k] //
129 // for i = k+1, ... , n-1. //
130 // The matrix U forms the upper triangular matrix, and the matrix L //
131 // forms the lower triangular matrix. //
132 // //
133 // Arguments: //
134 // float *A Pointer to the first element of the matrix A[n][n]. //
135 // int n The number of rows or columns of the matrix A. //
136 // //
137 // Return Values: //
138 // 0 Success //
139 // -1 Failure - The matrix A is singular. //
140 // //
141 // Example: //
142 // #define N //
143 // float A[N][N]; //
144 // //
145 // (your code to intialize the matrix A) //
146 // //
147 // err = Doolittle_LU_Decomposition(&A[0][0], N); //
148 // if (err < 0) printf(" Matrix A is singular\n"); //
149 // else { printf(" The LU decomposition of A is \n"); //
150 // ... //
152 // //
153 int Doolittle_LU_Decomposition(float *A, int n)
154 {
155  int i, j, k, p;
156  float *p_k, *p_row, *p_col;
157 
158 // For each row and column, k = 0, ..., n-1,
159 // find the upper triangular matrix elements for row k
160 // and if the matrix is non-singular (nonzero diagonal element).
161 // find the lower triangular matrix elements for column k.
162 
163  for (k = 0, p_k = A; k < n; p_k += n, k++) {
164  for (j = k; j < n; j++) {
165  for (p = 0, p_col = A; p < k; p_col += n, p++)
166  *(p_k + j) -= *(p_k + p) * *(p_col + j);
167  }
168  if ( *(p_k + k) == 0.0 ) return -1;
169  for (i = k+1, p_row = p_k + n; i < n; p_row += n, i++) {
170  for (p = 0, p_col = A; p < k; p_col += n, p++)
171  *(p_row + k) -= *(p_row + p) * *(p_col + k);
172  *(p_row + k) /= *(p_k + k);
173  }
174  }
175  return 0;
176 }
177 
178 
180 // int Doolittle_LU_Solve(float *LU, float *B, float *x, int n) //
181 // //
182 // Description: //
183 // This routine uses Doolittle's method to solve the linear equation //
184 // Ax = B. This routine is called after the matrix A has been decomposed //
185 // into a product of a unit lower triangular matrix L and an upper //
186 // triangular matrix U without pivoting. The argument LU is a pointer to //
187 // the matrix the subdiagonal part of which is L and the superdiagonal //
188 // together with the diagonal part is U. (The diagonal part of L is 1 and //
189 // is not stored.) The matrix A = LU. //
190 // The solution proceeds by solving the linear equation Ly = B for y and //
191 // subsequently solving the linear equation Ux = y for x. //
192 // //
193 // Arguments: //
194 // float *LU Pointer to the first element of the matrix whose elements //
195 // form the lower and upper triangular matrix factors of A. //
196 // float *B Pointer to the column vector, (n x 1) matrix, B //
197 // float *x Solution to the equation Ax = B. //
198 // int n The number of rows or columns of the matrix LU. //
199 // //
200 // Return Values: //
201 // 0 Success //
202 // -1 Failure - The matrix A is singular. //
203 // //
204 // Example: //
205 // #define N //
206 // float A[N][N], B[N], x[N]; //
207 // //
208 // (your code to create matrix A and column vector B) //
209 // err = Doolittle_LU_Decomposition(&A[0][0], N); //
210 // if (err < 0) printf(" Matrix A is singular\n"); //
211 // else { //
212 // err = Doolittle_LU_Solve(&A[0][0], B, x, n); //
213 // if (err < 0) printf(" Matrix A is singular\n"); //
214 // else printf(" The solution is \n"); //
215 // ... //
216 // } //
218 // //
219 int Doolittle_LU_Solve(float *LU, float B[], float x[], int n)
220 {
221 
222 // Solve the linear equation Lx = B for x, where L is a lower
223 // triangular matrix with an implied 1 along the diagonal.
224 
225  Unit_Lower_Triangular_Solve(LU, B, x, n);
226 
227 // Solve the linear equation Ux = y, where y is the solution
228 // obtained above of Lx = B and U is an upper triangular matrix.
229 
230  return Upper_Triangular_Solve(LU, x, x, n);
231 }
232 
233 int invertMatrix(float *LU, float *invA, float *I)
234 {
235  int i, j;
236  // float I[4][4] = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
237  float resultColumn[4];
238 
239  for (i = 0; i < 4; ++i)
240  {
241  int res = Doolittle_LU_Solve(LU, I + i*4, resultColumn, 4);
242 
243  if (res != 0) return res;
244  for (j = 0; j < 4; ++j)
245  *(invA + i + j * 4) = resultColumn[j];
246  }
247 
248  return 0;
249 }
250 
251 //float A[4][4] = {{1, 1, 1, 1}, {1, 4, 2, 3}, {1, 2, 1, 2}, {1, 1, 1, 0}};
252 //float invA[4][4]= {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}};
253 
254 int fun(float *A, float *invA, float *b, float *x, float *I)
255 {
256  int res = Doolittle_LU_Decomposition((float *)A, 4);
257 
258  if (res != 0) return res;
259 
260  // float b[4] = {63, 105, 48, 186};
261  // float x[4];
262 
263  res = Doolittle_LU_Solve((float *)A, b, x, 4);
264 
265  if (res != 0) return res;
266 
267  res = invertMatrix((float *)A, (float *)invA, I);
268 
269  return res;
270 }
Upper_Triangular_Solve
int Upper_Triangular_Solve(float *U, float B[], float x[], int n)
Definition: LUdecomposition.c:36
A
#define A
Definition: generate.c:13
invertMatrix
int invertMatrix(float *LU, float *invA, float *I)
Definition: LUdecomposition.c:233
k
static const uint32_t k[]
Definition: sha-256.c:22
L
#define L
Definition: spmv.h:13
Doolittle_LU_Solve
int Doolittle_LU_Solve(float *LU, float B[], float x[], int n)
Definition: LUdecomposition.c:219
Doolittle_LU_Decomposition
int Doolittle_LU_Decomposition(float *A, int n)
Definition: LUdecomposition.c:153
x
x
Return the smallest n such that 2^n >= _x.
Definition: math_function.hpp:170
fun
int fun(float *A, float *invA, float *b, float *x, float *I)
Definition: LUdecomposition.c:254
B
#define B
Definition: generate.c:14
Unit_Lower_Triangular_Solve
void Unit_Lower_Triangular_Solve(float *L, float B[], float x[], int n)
Definition: LUdecomposition.c:86
Generated on Mon Feb 12 2024 13:02:47 for PandA-2024.02 by doxygen 1.8.13