RcdMathLib_doc
Open Source Library for Linear and Non-linear Algebra
solve.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 
25 #include <stdio.h>
26 
27 #include "solve.h"
28 #include "matrix.h"
29 #include "qr_givens.h"
30 #include "qr_common.h"
31 #include "lu_decomp.h"
32 #include "qr_householder.h"
34 
35 int8_t solve(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m],
36  matrix_t x_sol[n], enum ALGORITHM algo)
37 {
38  int8_t status;
39 
40  switch (algo) {
41  case Moore_Penrose:
42  {
43  //puts("Moore Penrose !!!!");
44  matrix_t pinv_A[n][m];
45  status = moore_penrose_get_pinv(m, n, A, pinv_A);
46  matrix_mul_vec(n, m, pinv_A, b, x_sol);
47  break;
48  }
49  case Householder:
50  //puts("Householder !!!!");
51  status = solve_householder(m, n, A, b, x_sol);
52  break;
53 
54  case Givens:
55  //puts("Givens !!!");
56  status = solve_givens(m, n, A, b, x_sol);
57  break;
58 
59  case Gauss:
60  //puts("Gauss !!!");
61  status = solve_lu_decomp(m, n, A, b, x_sol);
62  break;
63 
64  default:
65  {
66  //puts("Default: Moore Penrose");
67  matrix_t pinv_A[n][m];
68  status = moore_penrose_get_pinv(m, n, A, pinv_A);
69  matrix_mul_vec(n, m, pinv_A, b, x_sol);
70  }
71 
72  }
73 
74  return status;
75 }
76 
77 int8_t solve_householder(uint8_t m, uint8_t n, matrix_t A[][n],
78  matrix_t b[m],
79  matrix_t x_sol[n])
80 {
81 
82  matrix_t Q[m][n];
83  matrix_t qt_b[n];
84 
85  if (m >= n) {
86  int8_t status;
87  status = qr_householder_decomp(m, n, A, n, Q, true);
88 
89  /* qt_b = Q'*b --> Rx = Q'b (R is stored in A) */
90  matrix_trans_mul_vec(m, n, Q, m, b, qt_b);
91  qr_common_backward_subst(n, n, A, qt_b, x_sol);
92 
93  return status;
94  }
95  else {
96  puts("[solve_householder]: The equation is not solvable !!!");
97 
98  return -2;
99  }
100 }
101 
102 int8_t solve_givens(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m],
103  matrix_t x_sol[n])
104 {
105  matrix_t Q[m][n];
106  matrix_t qt_b[n];
107 
108  if (m >= n) {
109  int8_t status;
110  status = qr_givens_decomp(m, n, A, n, Q, true);
111 
112  /* qt_b = Q'*b --> Rx = Q'b (R is stored in A) */
113  matrix_trans_mul_vec(m, n, Q, m, b, qt_b);
114  qr_common_backward_subst(n, n, A, qt_b, x_sol);
115 
116  return status;
117  }
118  else {
119  puts("[solve_givens]: The equation is not solvable !!!");
120  return -2;
121  }
122 }
123 
124 int8_t solve_lu_decomp(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m],
125  matrix_t x_sol[n])
126 {
127 
128  if (m == n) {
129  //puts("solve_lu_decomp !!!");
130  matrix_t L[m][m];
131  matrix_t P[m][m];
132  lu_decomp(m, A, L, P);
133 
134  /* lu_b = inv(L)*P*b */
135  matrix_t inv_L[m][m];
136  matrix_get_inv_low_triang(m, m, L, inv_L);
137  matrix_t tmp[m][m];
138  matrix_mul(m, m, inv_L, m, m, P, tmp);
139  matrix_t lu_b[m];
140  matrix_mul_vec(m, m, tmp, b, lu_b);
141 
142  /* U*x = inv(L)*P*b = lu_b*/
143  qr_common_backward_subst(n, n, A, lu_b, x_sol);
144  return 1;
145  }
146 
147  else {
148  puts(
149  "[solve_lu_decomp]: Only quadratic linear equation systems are supported !!!");
150  return -1;
151  }
152 }
moore_penrose_get_pinv
int8_t moore_penrose_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t pinv_A[n][m])
Calculate the Moore–Penrose inverse of a rectangular matrix.
Definition: moore_penrose_pseudo_inverse.c:38
matrix_mul_vec
void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t vec[n], matrix_t dst_arr[m])
Compute the multiplication of a matrix with a column vector.
Definition: matrix.c:434
Givens
Givens algorithm.
Definition: pseudo_inverse.h:33
qr_householder_decomp
int8_t qr_householder_decomp(uint8_t m, uint8_t n, matrix_t A[][n], uint8_t q_col_num, matrix_t Q[][q_col_num], bool reduced)
Computes the QR decomposition of the matrix A by using the Householder algorithm.
Definition: qr_householder.c:33
qr_householder.h
Householder algorithm for the QR-decomposition.
Moore_Penrose
Moore–Penrose algorithm.
Definition: pseudo_inverse.h:31
lu_decomp.h
Computes the LU decomposition of the matrix.
qr_common_backward_subst
void qr_common_backward_subst(uint8_t m, uint8_t n, matrix_t U[][n], matrix_t b[m], matrix_t x_sol[m])
Implements the backward substitution algorithm.
Definition: qr_common.c:28
solve
int8_t solve(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m], matrix_t x_sol[n], enum ALGORITHM algo)
Solve an (m n) linear system Ax = b by using the Moore–Penrose, Householder, or the Givens algorithm...
Definition: solve.c:35
matrix_trans_mul_vec
void matrix_trans_mul_vec(uint8_t m, uint8_t n, matrix_t A[m][n], uint8_t b_size, matrix_t b_vec[m], matrix_t c_vec[n])
Compute the multiplication of transposed matrix with column vector.
Definition: matrix.c:511
solve.h
Enables to solve systems of linear equations Ax = b for x.
matrix.h
Matrix computations.
moore_penrose_pseudo_inverse.h
Moore–Penrose algorithm to compute the pseudo-inverse of a matrix.
Householder
Householder algorithm.
Definition: pseudo_inverse.h:32
solve_lu_decomp
int8_t solve_lu_decomp(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m], matrix_t x_sol[n])
Solve an (m n) linear system Ax = b, using the Gaussian Elimination with pivoting algorithm.
Definition: solve.c:124
matrix_get_inv_low_triang
void matrix_get_inv_low_triang(uint8_t m, uint8_t n, matrix_t L[][n], matrix_t inv_L[][m])
Computes the inverse a lower triangular matrix.
Definition: matrix.c:817
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
ALGORITHM
ALGORITHM
Possible algorithms to compute the pseudo-inverse matrix.
Definition: pseudo_inverse.h:30
qr_common.h
Common definitions and implementations for the QR-decomposition. Provide necessary methods to constru...
Gauss
Gaussian elimination with pivoting algorithm.
Definition: pseudo_inverse.h:34
solve_givens
int8_t solve_givens(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m], matrix_t x_sol[n])
Solve an (m n) linear system Ax = b, using the Givens algorithm.
Definition: solve.c:102
qr_givens_decomp
int8_t qr_givens_decomp(uint8_t m, uint8_t n, matrix_t A[][n], uint8_t q_col_num, matrix_t Q[][q_col_num], bool reduced)
Computes the QR decomposition of the matrix A by using the Givens algorithm.
Definition: qr_givens.c:33
matrix_t
#define matrix_t
Define the data type of the matrix elements.
Definition: matrix.h:38
solve_householder
int8_t solve_householder(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m], matrix_t x_sol[n])
Solve an (m n) linear system Ax = b, using the Householder algorithm.
Definition: solve.c:77
qr_givens.h
Givens algorithm for the QR-decomposition. Provide necessary methods to construct Q- and R- matrices ...
matrix_mul
void matrix_mul(uint8_t a_line_num, uint8_t a_col_num, matrix_t a_matrix[a_line_num][a_col_num], uint8_t b_line_num, uint8_t b_col_num, matrix_t b_matrix[b_line_num][b_col_num], matrix_t dest_matrix[a_line_num][b_col_num])
Compute the multiplication of two matrices.
Definition: matrix.c:363