RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
lu_decomp.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
3  * 2020 Freie Universität Berlin
4  *
5  * This file is subject to the terms and conditions of the GNU Lesser General
6  * Public License v2.1. See the file LICENSE in the top level directory for more
7  * details.
8  */
9 
24 #include <float.h>
25 #include <stdio.h>
26 
27 #include "matrix.h"
28 #include "vector.h"
29 
30 // U will be saved in A
31 uint8_t lu_decomp(uint8_t n, matrix_t A[][n], matrix_t L[][n], matrix_t P[][n])
32 {
33  matrix_t abs_max_col_elem;
34  uint8_t k;
35  uint8_t pivot_index;
36  uint8_t changes;
37  matrix_t multipliers_vec[n - 1];
38 
39  matrix_get_diag_mat(n, n, 1, L);
40  matrix_get_diag_mat(n, n, 1, P);
41  changes = 0;
42  for (uint8_t i = 0; i < n - 1; i++) {
43  abs_max_col_elem =
45  n, n, A, i, i, &k);
46  if (abs_max_col_elem <= FLT_EPSILON) {
47  continue;
48  }
49 
50  pivot_index = k;
51 
52  if (pivot_index != i) {
53  //Swap the rows at the position i and pivot_index of the matrix A(i, i:n)
54  matrix_part_swap_rows(n, A, i, k, i, n - 1);
55 
56  //Swap the matrix P at the rows i and pivot_index
57  matrix_swap_rows(n, P, i, k);
58 
59  //Swap the rows of L in columns 0 through i-1. //FIXED
60  if (i >= 1) {
61 
62  matrix_part_swap_rows(n, L, i, k, 0, i - 1);
63  }
64  changes++;
65  } //if
66 
67  for (uint8_t j = i + 1; j < n; j++) {
68  //compute multipliers
69  multipliers_vec[j - 1] = A[j][i] / A[i][i];
70  //update the matrix A
71  for (uint8_t l = i + 1; l < n; l++) {
72  A[j][l] = A[j][l]
73  - multipliers_vec[j - 1]
74  * A[i][l];
75  }
76  //zero out the elements of column i, from row i+1 to n-1.
77  A[j][i] = 0;
78  // Set the matrix L with the multipliers
79  L[j][i] = multipliers_vec[j - 1];
80  }
81 
82  } //for
83 
84  return changes;
85 }
matrix_swap_rows
void matrix_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j)
Swaps two rows of a matrix.
Definition: matrix.c:725
matrix.h
Matrix computations.
matrix_get_diag_mat
void matrix_get_diag_mat(uint8_t m, uint8_t n, matrix_t value, matrix_t diag_matrix[m][n])
Create a diagonal matrix with a specified value.
Definition: matrix.c:594
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
vector.h
Vector computations.
matrix_part_swap_rows
void matrix_part_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j, uint8_t col_begin, uint8_t col_end)
Swaps two rows of a sub-matrix.
Definition: matrix.c:736
matrix_get_abs_max_elem_and_index_in_part_column
matrix_t matrix_get_abs_max_elem_and_index_in_part_column(uint8_t max_m, uint8_t max_n, matrix_t matrix[max_m][max_n], uint8_t row_num, uint8_t col_num, uint8_t *index)
Get the maximum absolute value and its position in a column vector in a sub-matrix.
Definition: matrix.c:703
lu_decomp
uint8_t lu_decomp(uint8_t n, matrix_t A[][n], matrix_t L[][n], matrix_t P[][n])
Computes the LU decomposition of the matrix.
Definition: lu_decomp.c:31