diff --git a/linear_algebra/basic_operations/doc.txt b/linear_algebra/basic_operations/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f585bf653f76a9550102a1b6181cc99a649f4cab
--- /dev/null
+++ b/linear_algebra/basic_operations/doc.txt
@@ -0,0 +1,21 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup basic_operations BASIC_OPERATIONS
+ * @ingroup linear_algebra
+ *
+ * @brief Matrix basic operations
+ *
+ * @details This module provides functions to perform basic matrix operations such as
+ * matrix addition, multiplication, or transposition.
+ *
+ * @author Zakaria Kasmi
+ *
+ */
diff --git a/linear_algebra/basic_operations/include/matrix.h b/linear_algebra/basic_operations/include/matrix.h
new file mode 100644
index 0000000000000000000000000000000000000000..fb2ae7cf33ba46a59cdd6c498f7863fae9d8bb79
--- /dev/null
+++ b/linear_algebra/basic_operations/include/matrix.h
@@ -0,0 +1,746 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup basic_operations
+ * @{
+ *
+ * @file
+ * @brief Matrix computations.
+ *
+ * @details Matrix computations include operations such as addition,
+ * subtraction, and transposition.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef MATRIX_H_
+#define MATRIX_H_
+
+#include <inttypes.h>
+#include <stdbool.h>
+
+/*!
+ \def matrix_t
+ @brief Define the data type of the matrix elements
+ @details The user can choose between a double or floating point arithmetic.
+ */
+//#define matrix_t float
+#ifndef matrix_t
+#define matrix_t double
+#endif
+
+/*!
+ \def MACHEPS
+ The upper bound on the relative error due to rounding in floating point arithmetic.
+ */
+#ifndef MACHEPS
+#define MACHEPS 2E-16
+#endif
+
+/*!
+ \def M_PI
+ Pi, the ratio of a circle's circumference to its diameter.
+ */
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+/**
+ * A structure to define the row and column number of a matrix.
+ */
+typedef struct {
+ uint8_t row_num; /**< the row number */
+ uint8_t col_num; /**< the column number */
+
+} matrix_dim_t;
+
+/**
+ * @brief Initialize all the elements of the matrix with a specified value.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[out] matrix[][] pointer to the matrix.
+ * @param[in] value value to be set.
+ *
+ */
+void matrix_init(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t value);
+
+/**
+ * @brief Clear all the elements of the vector.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[out] matrix[][] pointer to the matrix.
+ *
+ */
+void matrix_clear(uint8_t m, uint8_t n, matrix_t matrix[m][n]);
+
+/**
+ * @brief Computes the transpose of a matrix.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] src_matrix[][] pointer to the matrix to transpose.
+ * @param[out] dest_matrix[][] pointer to the destination matrix.
+ *
+ */
+void matrix_transpose(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
+ matrix_t dest_matrix[n][m]);
+
+/**
+ * @brief Computes the in-place transpose of a matrix.
+ * Transpose the matrix without auxiliary memory.
+ *
+ * @note This function is limited to square matrices.
+ *
+ * @param[in] m row and column number of the matrix.
+ * @param[in,out] matrix[][] pointer to the matrix to transpose.
+ *
+ */
+void matrix_in_place_transpose(uint8_t m, matrix_t matrix[][m]);
+
+/**
+ * @brief Copy the elements of a matrix to another matrix.
+ *
+ * @param[in] m row number of the matrix to copy.
+ * @param[in] n column number of the matrix to copy.
+ * @param[in] src_matrix[][] pointer to the source matrix
+ * @param[out] dest_matrix[][] pointer to the destination matrix.
+ *
+ */
+void matrix_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
+ matrix_t dest_matrix[m][n]);
+/**
+ * @brief Copy a part of a matrix to another matrix or sub-matrix.
+ * A part of the source matrix can be copied in a sub-part of
+ * the destination matrix (sub-matrix). The source and destination
+ * sub-matrices are limited by the row and column indices.
+ *
+ * @param[in] m row number of the source matrix.
+ * @param[in] n column number of the source matrix.
+ * @param[in] src_matrix[][] pointer to the source matrix.
+ * @param[in] start_row_ind the start index of the rows of the source sub-matrix.
+ * @param[in] end_row_ind the end index of the rows of the source sub-matrix.
+ * @param[in] start_col_ind the start index of the columns of the source sub-matrix.
+ * @param[in] end_col_ind the end index of the columns of the source sub-matrix.
+ * @param[in] dest_row_num the row number of the destination sub-matrix.
+ * @param[in] dest_col_num the column number of the destination sub-matrix.
+ * @param[out] dest_matrix[][] pointer to the destination (sub)-matrix.
+ *
+ */
+void matrix_part_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
+ uint8_t start_row_ind, uint8_t end_row_ind,
+ uint8_t start_col_ind, uint8_t end_col_ind,
+ uint8_t dest_row_num, uint8_t dest_col_num,
+ matrix_t dest_matrix[][dest_col_num]);
+
+/**
+ * @brief Compute the rank of a matrix.
+ * The SVD must be previously invoked to get the singular values of
+ * the matrix.
+ *
+ * @note This function should be invoked after the call of the @ref svd method.
+ *
+ * @param[in] m row number of the source matrix.
+ * @param[in] n column number of the source matrix.
+ * @param[in] singl_values_arr[] array containing the singular values of the
+ * matrix.
+ * @param[in] length length of the singular values array.
+ *
+ * @return the rank of the matrix.
+ *
+ */
+uint8_t matrix_get_rank(uint8_t m, uint8_t n, matrix_t singl_values_arr[],
+ uint8_t length);
+
+/**
+ * @brief Compute the addition of two matrices.
+ * Add matrix B to matrix A and return the result in A_plus_B matrix.
+ *
+ * @param[in] m row number of the matrix to add.
+ * @param[in] n column number of the matrix to add.
+ * @param[in] A[][] pointer to the first matrix.
+ * @param[in] B[][] pointer to the second matrix.
+ * @param[out] A_plus_B[][] pointer to the destination matrix.
+ *
+ */
+void matrix_add(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n],
+ matrix_t A_plus_B[m][n]);
+
+/**
+ * @brief Add a number to diagonal elements of a matrix.
+ *
+ * @param[in] n column number of the matrix.
+ * @param[in,out] A[][] pointer to the matrix.
+ * @param[in] diag_el_num number of diagonal elements to overwrite.
+ * @param[in] value the value to add to the diagonal elements.
+ *
+ */
+void matrix_add_to_diag(uint8_t n, matrix_t A[][n], uint8_t diag_el_num,
+ matrix_t value);
+
+/**
+ * @brief Compute the subtraction of two matrices.
+ * Subtract matrix B from matrix A and return the result in
+ * A_minus_B matrix.
+ *
+ * @param[in] m row number of the matrix to add.
+ * @param[in] n column number of the matrix to add.
+ * @param[in] A[][] pointer to the first matrix.
+ * @param[in] B[][] pointer to the second matrix.
+ * @param[out] A_minus_B[][] pointer to the destination matrix.
+ *
+ */
+void matrix_sub(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n],
+ matrix_t A_minus_B[m][n]);
+
+/**
+ * @brief Compute the multiplication of two matrices.
+ *
+ * @param[in] a_line_num row number of the first matrix.
+ * @param[in] a_col_num column number of the first matrix.
+ * @param[in] a_matrix[][] pointer to the first matrix.
+ * @param[in] b_line_num row number of the second matrix.
+ * @param[in] b_col_num column number of the second matrix.
+ * @param[in] b_matrix[][] pointer to the second matrix.
+ * @param[out] dest_matrix[][] pointer to the destination matrix.
+ *
+ */
+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]);
+
+/**
+ * @brief Compute the partial multiplication of two matrices.
+ *
+ * @details Enables the calculation of matrix product of parts of two matrices.
+ *
+ * @param[in] a_col_num_max column number of the first matrix.
+ * @param[in] a_matrix[][] pointer to the first matrix.
+ * @param[in] b_col_num_max column number of the second matrix.
+ * @param[in] b_matrix[][] pointer to the second matrix.
+ * @param[in] a_start_row_ind row begin of the first, partial matrix.
+ * @param[in] a_end_row_ind row end of the first, partial matrix.
+ * @param[in] a_start_col_ind column begin of the first, partial matrix.
+ * @param[in] a_end_col_ind column end of the first, partial matrix.
+ * @param[in] b_start_row_ind row begin of the second, partial matrix.
+ * @param[in] b_end_row_ind row end of the second, partial matrix.
+ * @param[in] b_start_col_ind column begin of the second, partial matrix.
+ * @param[in] b_end_col_ind column end of the second, partial matrix.
+ * @param[in] dest_col_size column size of the destination matrix.
+ * @param[out] dest_matrix[][] pointer to the destination matrix.
+ */
+void matrix_part_mul(uint8_t a_col_num_max, matrix_t a_matrix[][a_col_num_max],
+ uint8_t b_col_num_max, matrix_t b_matrix[][b_col_num_max],
+ uint8_t a_start_row_ind, uint8_t a_end_row_ind,
+ uint8_t a_start_col_ind, uint8_t a_end_col_ind,
+ uint8_t b_start_row_ind, uint8_t b_end_row_ind,
+ uint8_t b_start_col_ind, uint8_t b_end_col_ind,
+ uint8_t dest_col_size,
+ matrix_t dest_matrix[][dest_col_size]);
+
+/**
+ * @brief Compute the multiplication of a matrix with a column vector.
+ * Return Am,n * Vn,1 = Bm,1, where the A is a (mxn)-matrix, V is a
+ * n-dimensional column vector, and the result is a m-dimensional
+ * column vector.
+ *
+ * @param[in] m row number of the matrix to multiply.
+ * @param[in] n column number of the matrix to multiply.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] vec pointer to the n-dimensional column vector.
+ * @param[out] dst_arr pointer to the destination m-dimensional column
+ * vector.
+ *
+ */
+void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ matrix_t vec[n], matrix_t dst_arr[m]);
+
+//A'(n,m)*b(m,1) = c(n,1)
+/**
+ * @brief Compute the multiplication of transposed matrix with column vector.
+ * @details Transpose A and return A'n,m * Vm,1 = Bn,1, where A is a
+ * (mxn)-matrix, V is a m-dimensional column vector, and the result is
+ * a n-dimensional column vector.
+ *
+ * @param[in] m row number of the matrix to transpose and multiply.
+ * @param[in] n column number of the matrix to transpose and multiply.
+ * @param[in] A[][] pointer to the matrix.
+ * @param[in] b_size size of the column vector.
+ * @param[in] b_vec[] pointer to the column vector of m-dimension.
+ * @param[out] c_vec[] pointer to the destination, column vector of n-dimension.
+ *
+ */
+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]);
+
+
+/**
+ * @brief Compute the multiplication of a column and row vector.
+ * Return Cm,1 * R1,n = Mm,n, where C is a m-dimensional column vector,
+ * R is a n-dimensional row vector, and the result is a (mxn)-matrix.
+ *
+ * @param[in] m row number of the column vector.
+ * @param[in] col_vec[] pointer to the column vector.
+ * @param[in] n column number of the row vector.
+ * @param[in] row_vec[] pointer to the row vector.
+ * @param[in] max_n column number of the result matrix.
+ * @param[out] res_mat[][] pointer to the (mxn) result matrix.
+ *
+ */
+void matrix_mul_col_vec_row_vec(uint8_t m, matrix_t col_vec[m], uint8_t n,
+ matrix_t row_vec[n], uint8_t max_n,
+ matrix_t res_mat[][max_n]);
+
+
+/**
+ * @brief Compute the multiplication of row vector and a matrix.
+ * Return R1,m * Am,n = R1,n, where R is a m-dimensional row vector,
+ * A is a (mxn)-matrix, and the result is a row vector of n-dimension.
+ *
+ * @param[in] m size of the row vector and row number of the
+ * matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] vec[] pointer to the row vector.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[out] dst_arr pointer to the destination row vector of
+ * n-dimension.
+ *
+ */
+void matrix_vec_mul_matr(uint8_t m, uint8_t n, matrix_t vec[m],
+ matrix_t matrix[m][n], matrix_t dst_arr[n]);
+
+
+/**
+ * @brief Compute the multiplication of a scalar with row vector and a matrix.
+ * Return scal*R1,m * Am,n = R1,n, where scal is a scalar, R is a
+ * m-dimensional row vector, A is a (mxn)-matrix, and the result is a
+ * row vector of n-dimension.
+ *
+ * @param[in] m size of the row vector and row number of the
+ * matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] scalar scalar value.
+ * @param[in] vec[] pointer to the row vector.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[out] dst_arr pointer to the destination row vector of
+ * n-dimension.
+ *
+ */
+void matrix_mul_scalar_vec_matr(uint8_t m, uint8_t n, matrix_t scalar,
+ matrix_t vec[m], matrix_t matrix[m][n],
+ matrix_t dst_arr[n]);
+
+
+/**
+ * @brief Compute the multiplication of the transpose of a matrix with itself.
+ * Transpose the matrix A and multiply it with the matrix A: A'*A.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] A[][] pointer to the matrix.
+ * @param[out] AT_mul_A[][] pointer to the destination matrix (A'*A).
+ *
+ */
+void matrix_trans_mul_itself(uint8_t m, uint8_t n, matrix_t A[m][n],
+ matrix_t AT_mul_A[n][n]);
+
+/**
+ * @brief Set all the diagonal elements of a matrix with a specified value.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] value value to be set.
+ * @param[in,out] diag_matrix[][] pointer to the matrix.
+ *
+ */
+void matrix_set_diag_elements(uint8_t m, uint8_t n, matrix_t value,
+ matrix_t diag_matrix[m][n]);
+
+/**
+ * @brief Set all the diagonal elements of a matrix with values that are saved
+ * in a vector.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in,out] diag_matrix[][] pointer to the matrix.
+ * @param[in] length size of the vector.
+ * @param[in] vec pointer to the vector containing diagonal
+ * elements.
+ *
+ */
+void matrix_get_diag_mat_new(uint8_t m, uint8_t n, matrix_t diag_matrix[m][n],
+ uint8_t length, matrix_t vec[]);
+
+/**
+ * @brief Create a diagonal matrix with a specified value.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] value value of the diagonal elements.
+ * @param[in,out] diag_matrix[][] pointer to the diagonal matrix.
+ *
+ */
+void matrix_get_diag_mat(uint8_t m, uint8_t n, matrix_t value,
+ matrix_t diag_matrix[m][n]);
+
+/**
+ * @brief Multiply all elements of a matrix with a specified value.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] mat_src pointer to the source matrix.
+ * @param[in] value multiplication factor.
+ * @param[out] mat_dest[][] pointer to the destination matrix.
+ *
+ */
+void matrix_mul_scalar(uint8_t m, uint8_t n, matrix_t mat_src[m][n],
+ matrix_t value, matrix_t mat_dest[m][n]);
+
+/**
+ * @brief Get a column of a matrix.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] col_num number of the requested column.
+ * @param[out] col_vec[] pointer to the column vector.
+ *
+ */
+void matrix_get_column_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t col_num, matrix_t col_vec[m]);
+
+/**
+ * @brief Get a part of a column of a matrix.
+ *
+ * @param[in] max_m row number of the matrix.
+ * @param[in] max_n column number of the matrix.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] col_num number of the requested column.
+ * @param[in] offset points to the start position of the column vector.
+ * @param[out] col_vec[] pointer to the column vector.
+ *
+ */
+void matrix_get_part_column_vec(uint8_t max_m, uint8_t max_n,
+ matrix_t matrix[max_m][max_n], uint8_t col_num,
+ uint8_t offset,
+ matrix_t col_vec[max_m - offset]);
+
+/**
+ * @brief Get the largest element of a column vector in a matrix.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] col_num column number.
+ *
+ * @return the largest element of a column.
+ *
+ */
+matrix_t matrix_get_max_elem_in_column(uint8_t m, uint8_t n,
+ matrix_t matrix[m][n], uint8_t col_num);
+
+/**
+ * @brief Get the maximum absolute value of a column vector in a matrix.
+ *
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] col_num column number.
+ *
+ * @return the maximum absolute value of a column.
+ *
+ */
+matrix_t matrix_get_abs_max_elem_in_column(uint8_t m, uint8_t n,
+ matrix_t matrix[m][n],
+ uint8_t col_num);
+
+/**
+ * @brief Get the largest element of a column vector in a sub-matrix.
+ *
+ * @param[in] max_m total row number of the matrix.
+ * @param[in] max_n total column number of the matrix.
+ * @param[in] matrix[][] pointer to the entire matrix.
+ * @param[in] row_num row number of the sub-matrix.
+ * @param[in] col_num column number of the sub-matrix.
+ *
+ * @return the largest element of a partial column.
+ *
+ */
+matrix_t matrix_get_max_elem_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);
+
+/**
+ * @brief Get the maximum absolute value of a column vector in a sub-matrix.
+ *
+ * @param[in] max_m total row number of the matrix.
+ * @param[in] max_n total column number of the matrix.
+ * @param[in] matrix[][] pointer to the entire matrix.
+ * @param[in] row_num row number of the sub-matrix.
+ * @param[in] col_num column number of the sub-matrix.
+ *
+ * @return the maximum absolute value of a partial column.
+ *
+ */
+matrix_t matrix_get_abs_max_elem_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);
+
+/**
+ * @brief Get the maximum absolute value and its position in a column vector in a sub-matrix.
+ *
+ * @param[in] max_m total row number of the matrix.
+ * @param[in] max_n total column number of the matrix.
+ * @param[in] matrix[][] pointer to the entire matrix.
+ * @param[in] row_num row number of the sub-matrix.
+ * @param[in] col_num column number of the sub-matrix.
+ * @param[out] index pointer to the variable holding the position of the maximum absolute
+ * value in the column vector of the sub-matrix.
+ *
+ * @return the maximum absolute value of a partial 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);
+/**
+ * @brief Swaps two rows of a matrix.
+ * @details Swaps the rows i and j of a matrix.
+ *
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] i the i-the row of the matrix.
+ * @param[in] j the j-the row of the matrix.
+ *
+ */
+void matrix_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j);
+
+/**
+ * @brief Swaps two rows of a sub-matrix.
+ * @details Swaps the rows i and j of a part of a matrix.
+ *
+ * @param[in] n column number of the entire matrix.
+ * @param[in, out] matrix[][] pointer to the entire matrix.
+ * @param[in] i the i-the row of the sub-matrix.
+ * @param[in] j the j-the row of the sub-matrix.
+ * @param[in] col_begin the column begin of the sub-matrix.
+ * @param[in] col_end the column end of the sub-matrix.
+ *
+ */
+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);
+/**
+ * @brief Get the 2-norm of a matrix that is equal to the largest singular value.
+ *
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] A[][] pointer to the matrix.
+ *
+ * @return the 2-norm of a matrix.
+ *
+ */
+double matrix_get_two_norm(uint8_t m, uint8_t n, matrix_t A[][n]);
+
+/**
+ * @brief Get the Frobenius norm of a matrix.
+ *
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] A[][] pointer to the matrix.
+ *
+ * @return the Frobenius norm a matrix.
+ *
+ */
+double matrix_get_frob_norm(uint8_t m, uint8_t n, matrix_t A[][n]);
+
+/**
+ * @brief Computes the inverse an upper triangular matrix.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] U[][] pointer to the upper triangular matrix.
+ * @param[out] inv_U[][] pointer to the inverse matrix.
+ *
+ */
+void matrix_get_inv_upp_triang(uint8_t m, uint8_t n, matrix_t U[][n],
+ matrix_t inv_U[][m]);
+/**
+ * @brief Computes the inverse a lower triangular matrix.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] L[][] pointer to the matrix.
+ * @param[out] inv_L[][] pointer to the inverse matrix.
+ *
+ */
+void matrix_get_inv_low_triang(uint8_t m, uint8_t n, matrix_t L[][n],
+ matrix_t inv_L[][m]);
+/**
+ * @brief Gets the upper triangular part of a matrix.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] A[][] pointer to the matrix.
+ * @param[out] tr_up_A[][] pointer to the upper triangular part of the matrix.
+ *
+ */
+void matrix_get_upp_triang(uint8_t m, uint8_t n, matrix_t A[][n],
+ matrix_t tr_up_A[][n]);
+/**
+ * @brief Gets the lower triangular part of a matrix.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] A[][] pointer to the matrix.
+ * @param[out] tr_low_A[][] pointer to the lower triangular part of the matrix.
+ *
+ */
+void matrix_get_low_triang(uint8_t m, uint8_t n, matrix_t A[][n],
+ matrix_t tr_low_A[][n]);
+
+/**
+ * @brief Display the values of the matrix elements.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the entire matrix.
+ *
+ */
+void matrix_print(uint8_t m, uint8_t n, matrix_t matrix[m][n]);
+
+/**
+ * @brief Display the values of the sub-matrix elements.
+ *
+ * @param[in] m total row number of the matrix.
+ * @param[in] n total column number of the matrix.
+ * @param[in] matrix[][] pointer to the entire matrix.
+ * @param[in] start_row_ind start row number of the sub-matrix.
+ * @param[in] end_row_ind end row number of the sub-matrix.
+ * @param[in] start_col_ind start column number of the sub-matrix.
+ * @param[in] end_col_ind end column number of the sub-matrix.
+ *
+ */
+void matrix_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t start_row_ind, uint8_t end_row_ind,
+ uint8_t start_col_ind, uint8_t end_col_ind);
+
+/**
+ * @brief Display the values of the matrix elements.
+ * This function allows the user to determine the precision as well as
+ * the with of the numbers to display.
+ *
+ * @note This function is more memory-consuming than @ref matrix_print.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the entire matrix.
+ * @param[in] before_dec the number of digits to be printed before the
+ * decimal point.
+ * @param[in] after_dec the number of digits to be printed after the
+ * decimal point.
+ *
+ */
+void matrix_flex_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t before_dec, uint8_t after_dec);
+
+
+/**
+ * @brief Display the values of the sub-matrix elements.
+ * This function allows the user to determine the precision as well as
+ * the with of the numbers to display.
+ *
+ * @note This function is more memory-consuming than @ref matrix_part_print.
+ *
+ * @param[in] m total row number of the matrix.
+ * @param[in] n total column number of the matrix.
+ * @param[in] matrix[][] pointer to the entire matrix.
+ * @param[in] start_row_ind start row number of the sub-matrix.
+ * @param[in] end_row_ind end row number of the sub-matrix.
+ * @param[in] start_col_ind start column number of the sub-matrix.
+ * @param[in] end_col_ind end column number of the sub-matrix.
+ * @param[in] before_dot the number of digits to be printed before the
+ * decimal point.
+ * @param[in] after_dot the number of digits to be printed after the
+ * decimal point.
+ *
+ */
+void matrix_flex_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t start_row_ind, uint8_t end_row_ind,
+ uint8_t start_col_ind, uint8_t end_col_ind,
+ uint8_t before_dot, uint8_t after_dot);
+
+
+/**
+ * @brief Compute the multiplication of a scalar with row vector and a
+ * sub-matrix.
+ * @details Return scal*R1,m * Am,n = R1,n, where scal is a scalar, R is a
+ * m-dimensional row vector, A is a (mxn)-matrix, and the result is a
+ * row vector of n-dimension.
+ *
+ * @param[in] max_m size of the row vector and row number of the
+ * matrix.
+ * @param[in] max_n column number of the matrix.
+ * @param[in] scalar scalar value.
+ * @param[in] vec[] pointer to the row vector.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] begin_row start row number of the sub-matrix.
+ * @param[in] begin_column start column number of the sub-matrix.
+ * @param[out] dst_arr pointer to the destination row vector of n-dimension.
+ *
+ */
+void matrix_part_mul_scalar_vec_matr(uint8_t max_m, uint8_t max_n,
+ matrix_t scalar, matrix_t vec[max_m],
+ matrix_t matrix[max_m][max_n],
+ uint8_t begin_row, uint8_t begin_column,
+ matrix_t dst_arr[max_n - begin_row]);
+
+/**
+ * @brief Get the value of a matrix at the position (i,j).
+ *
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] i row number.
+ * @param[in] j column number.
+ *
+ * @return the value of the matrix at the row i and column j.
+ *
+ */
+matrix_t matrix_read(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i,
+ uint8_t j);
+/**
+ * @brief Write a value in a matrix at the position (i,j).
+ *
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[out] matrix[][] pointer to the matrix.
+ * @param[in] i row number.
+ * @param[in] j column number.
+ * @param[in] val value to write in the matrix.
+ *
+ */
+void matrix_write(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i,
+ uint8_t j, matrix_t val);
+
+#endif /* MATRIX_H_ */
diff --git a/linear_algebra/basic_operations/include/vector.h b/linear_algebra/basic_operations/include/vector.h
new file mode 100644
index 0000000000000000000000000000000000000000..64b1295cf574b99d22cddda3991b51b14f468c13
--- /dev/null
+++ b/linear_algebra/basic_operations/include/vector.h
@@ -0,0 +1,340 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup basic_operations
+ * @{
+ *
+ * @file
+ * @brief Vector computations.
+ * @details Vector computations include operations such as addition,
+ * subtraction, and inner product (dot product).
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef VECTOR_H_
+#define VECTOR_H_
+
+#include <inttypes.h>
+#include <stdbool.h>
+
+/** @brief Define the data type of the vector elements */
+//#define vector_t float
+#ifndef vector_t
+#define vector_t double
+#endif
+
+/**
+ * @brief Clear all the elements of the vector.
+ *
+ * @param[in] size size of the vector.
+ * @param[in] arr[] pointer to the vector.
+ *
+ */
+void vector_clear(uint8_t size, vector_t arr[]);
+
+/**
+ * @brief Copy the elements of the source vector to the destination vector.
+ *
+ * @param[in] size number of elements to copy.
+ * @param[in] src_arr[] pointer to the source vector.
+ * @param[in] dest_arr[] pointer to the destination vector.
+ *
+ */
+void vector_copy(uint8_t size, vector_t src_arr[], vector_t dest_arr[]);
+
+/**
+ * @brief Clear all the elements of the vector.
+ *
+ * @param[in] size size of the vector.
+ * @param[in] arr[] pointer to the vector.
+ *
+ */
+void vector_clear(uint8_t size, vector_t arr[]);
+
+/**
+ * @brief Compute the 2-norm norm of a vector.
+ *
+ * @param[in] length size of the vector.
+ * @param[in] arr[] pointer to the vector.
+ *
+ * @return the 2-norm of the vector.
+ *
+ */
+vector_t vector_get_norm2(uint8_t length, vector_t arr[]);
+
+/**
+ * @brief Compute the squared 2-norm norm of a vector .
+ *
+ * @param[in] length size of the vector.
+ * @param[in] arr[] pointer to the vector.
+ *
+ * @return the squared 2-norm of the vector.
+ *
+ */
+vector_t vector_get_square_norm2(uint8_t length, vector_t arr[]);
+
+/**
+ * @brief Compute the sum of the elements of a vector.
+ *
+ * @param[in] length size of the vector.
+ * @param[in] arr[] pointer to the vector.
+ *
+ * @return the sum of the elements of the vector.
+ *
+ */
+vector_t vector_get_sum(uint8_t length, vector_t arr[]);
+
+/**
+ * @brief Compute the average or mean value of a vector.
+ *
+ * @param[in] length size of the vector.
+ * @param[in] arr[] pointer to the vector.
+ *
+ * @return the mean value of the vector.
+ *
+ */
+vector_t vector_get_mean_value(uint8_t length, vector_t arr[]);
+
+/**
+ * @brief Compute the subtraction of two vectors.
+ * Substact b_vec from a_vec and return the result in a_minus_b.
+ *
+ * @param[in] size number of elements to subtract.
+ * @param[in] a_vec[] pointer to the first vector.
+ * @param[in] b_vec[] pointer to the second vector.
+ * @param[out] a_minus_b[] pointer to the destination vector.
+ *
+ */
+void vector_sub(uint8_t size, vector_t a_vec[], vector_t b_vec[],
+ vector_t a_minus_b[]);
+
+/**
+ * @brief Compute the addition of two vectors.
+ * Add b_vec to a_vec and return the result in a_plus_b.
+ *
+ * @param[in] size number of elements to subtract.
+ * @param[in] a_vec[] pointer to the first vector.
+ * @param[in] b_vec[] pointer to the second vector.
+ * @param[out] a_plus_b_vec[] pointer to the destination vector.
+ *
+ */
+void vector_add(uint8_t size, vector_t a_vec[size], vector_t b_vec[size],
+ vector_t a_plus_b_vec[size]);
+/**
+ * @brief Compute the multiplication of two vectors.
+ * Multiple vectors a_vec and b_vec element by element and return
+ * the result in a_mul_b.
+ *
+ * @param[in] size number of elements to multiply.
+ * @param[in] a_vec[] pointer to the first vector.
+ * @param[in] b_vec[] pointer to the second vector.
+ * @param[out] a_mul_b_vec[] pointer to the destination vector.
+ *
+ */
+void vector_mul(uint8_t size, vector_t a_vec[size], vector_t b_vec[size],
+ vector_t a_mul_b_vec[size]);
+
+/**
+ * @brief Compute the square of a vector.
+ * Square the elements of vector vec and return the result in
+ * square_vec.
+ *
+ * @param[in] n number of elements to square.
+ * @param[in] vec[] pointer to the source vector.
+ * @param[out] square_vec[] pointer to the destination vector.
+ *
+ */
+void vector_square(uint8_t n, vector_t vec[n], vector_t square_vec[n]);
+
+/**
+ * @brief Compute the product of a vector with a real number.
+ * Multiple the elements of a vector with a scalar and return the
+ * result in the vector itself.
+ *
+ * @param[in] size number of elements to multiply with a scalar.
+ * @param[in, out] a_vec[] pointer to the source/destination vector.
+ * @param[in] scl a scalar.
+ *
+ */
+void vector_in_place_scalar_mul(uint8_t size, vector_t a_vec[size],
+ vector_t scl);
+
+/**
+ * @brief Compute the product of a vector with a real number.
+ * Multiple the elements of a vector with a scalar and return the
+ * result in other vector.
+ *
+ * @param[in] size number of elements to multiply with a scalar.
+ * @param[in] src_vec[] pointer to the source vector.
+ * @param[in] scl a scalar.
+ * @param[out] dest_vec pointer to the destination vector.
+ *
+ */
+void vector_scalar_mul(uint8_t size, vector_t src_vec[size], vector_t scl,
+ vector_t dest_vec[]);
+
+/**
+ * @brief Compute the division of a vector with a real number.
+ * Divide the elements of a vector with a scalar and return the
+ * result in the vector itself.
+ *
+ * @param[in] size number of elements to divide with a scalar.
+ * @param[in, out] a_vec[] pointer to the source/destination vector.
+ * @param[in] scl a scalar.
+ *
+ */
+void vector_scalar_div(uint8_t size, vector_t a_vec[size], vector_t scl);
+
+//Euclidean distance: d = sum((x-y).^2).^0.5
+
+/**
+ * @brief Compute the Euclidean distance between two vectors.
+ *
+ * @param[in] length size of the vector.
+ * @param[in] vec1[] pointer to the first vector.
+ * @param[in] vec2[] pointer to the second vector.
+ *
+ * @return the Euclidean distance.
+ *
+ */
+vector_t vector_get_euclidean_distance(uint8_t length, vector_t vec1[],
+ vector_t vec2[]);
+/**
+ * @brief Compute the dot product of two vectors.
+ *
+ * @param[in] n size of the vectors.
+ * @param[in] vec1[] pointer to the first vector.
+ * @param[in] vec2[] pointer to the second vector.
+ *
+ * @return the scalar product of two vectors.
+ *
+ */
+vector_t vector_get_scalar_product(uint8_t n, vector_t vec1[n],
+ vector_t vec2[n]);
+
+/**
+ * @brief Determine the equality of two vectors.
+ *
+ * @param[in] length size of the vector.
+ * @param[in] vec_1[] pointer to the first vector.
+ * @param[in] vec_2[] pointer to the second vector.
+ *
+ * @return true, if the two vectors are equal.
+ * @return false, if not.
+ *
+ */
+bool vector_is_equal(uint16_t length, vector_t vec_1[], vector_t vec_2[]);
+
+/**
+ * @brief Determine the equality of two vectors of type uint32_t.
+ *
+ * @param[in] length size of the vector.
+ * @param[in] vec_1[] pointer to the first vector.
+ * @param[in] vec_2[] pointer to the second vector.
+ *
+ * @return true, if the two vectors are equal.
+ * @return false, if not.
+ *
+ */
+bool vector_uint32_is_equal(uint32_t length, uint32_t vec_1[], uint32_t vec_2[]);
+
+/**
+ * @brief Determine the index of the vector elements before sorting.
+ * Determine the index of the elements of a sorted vector. These
+ * indices correspond to the positions of the elements in the unsorted
+ * vector.
+ *
+ * @param[in] k size of the unsorted vector.
+ * @param[in] n size of the sorted vector.
+ * @param[in] unsorted_vector[] pointer to the unsorted vector.
+ * @param[in] sorted_vector[] pointer to the sorted vector.
+ * @param[out] index_vector[] pointer to the index vector.
+ *
+ */
+void vector_get_index_vector(uint8_t k, uint8_t n, vector_t unsorted_vector[n],
+ vector_t sorted_vector[n], uint8_t index_vector[n]);
+
+/**
+ * @brief Compute the maximal value and its index of a vector.
+ *
+ *
+ * @param[in] length vector size.
+ * @param[in] vec[] pointer to the vector.
+ * @param[in] index pointer to the index.
+ *
+ * @return the maximal value of the vector.
+ *
+ */
+vector_t vector_get_max_and_index(uint8_t length, vector_t vec[],
+ uint8_t *index);
+
+/**
+ * @brief Compute the residual of two vectors.
+ *
+ *
+ * @param[in] length vector size.
+ * @param[in] a_vec[] pointer to the first vector.
+ * @param[in] b_vec[] pointer to the second vector.
+ *
+ * @return the residual of the two vectors.
+ *
+ */
+vector_t vector_get_residual(uint8_t length, vector_t a_vec[], vector_t b_vec[]);
+
+/**
+ * @brief Get the elements of the vector by an index vector.
+ *
+ * @param[in] src_vec[] pointer to source vector.
+ * @param[in] k size of the index vector.
+ * @param[in] index_vec[] pointer to the index vector.
+ * @param[out] dst_vec[] pointer to the destination vector
+ *
+ */
+void vector_get_elements(vector_t src_vec[], uint8_t k, uint8_t index_vec[],
+ vector_t dst_vec[]);
+
+/**
+ * @brief Display the values of the vector's elements.
+ *
+ * @param[in] length size of the vector to display.
+ * @param[in] arr pointer to the vector.
+ *
+ */
+void vector_print(uint32_t length, vector_t arr[]);
+
+/**
+ * @brief Display the values of the vector's elements of type uint8_t.
+ *
+ * @param[in] length size of the vector to display.
+ * @param[in] arr pointer to the vector.
+ *
+ */
+void vector_print_u8_array(uint32_t length, uint8_t arr[]);
+
+/**
+ * @brief Display the values of the vector's elements.
+ * This function allows the user to determine the precision as well as
+ * the with of the numbers to display.
+ *
+ * @param[in] length size of the vector to display.
+ * @param[in] arr pointer to the vector.
+ * @param[in] before_dot the number of digits to be printed before the
+ * decimal point.
+ * @param[in] after_dot the number of digits to be printed after the
+ * decimal point.
+ *
+ */
+void vector_flex_print(uint32_t length, vector_t arr[], uint8_t before_dot,
+ uint8_t after_dot);
+
+#endif /* VECTOR_H_ */
diff --git a/linear_algebra/basic_operations/matrix.c b/linear_algebra/basic_operations/matrix.c
new file mode 100644
index 0000000000000000000000000000000000000000..c9475df58bc540b5cb012312306545f4fa6eb9fb
--- /dev/null
+++ b/linear_algebra/basic_operations/matrix.c
@@ -0,0 +1,905 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Matrix computations.
+ *
+ * @details Matrix computations include operations such as addition,
+ * subtraction, and transposition.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdint.h>
+#include <string.h>
+#include <stdio.h>
+#include <math.h>
+#include <stdint.h>
+#include <stdbool.h>
+#include <float.h>
+
+#include "matrix.h"
+#include "utils.h"
+#include "svd.h"
+
+static int32_t max(int32_t a, int32_t b)
+{
+ if (a > b) {
+ return a;
+ }
+ else {
+ return b;
+ }
+}
+
+void matrix_clear(uint8_t m, uint8_t n, matrix_t matrix[m][n])
+{
+ memset(matrix, 0, sizeof(matrix[0][0]) * m * n);
+}
+
+void matrix_init(uint8_t m, uint8_t n, matrix_t matrix[m][n], matrix_t value)
+{
+ memset(matrix, value, sizeof(matrix[0][0]) * m * n);
+}
+
+void matrix_transpose(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
+ matrix_t dest_matrix[n][m])
+{
+ int i, j;
+
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < n; j++) {
+ dest_matrix[j][i] = src_matrix[i][j];
+ }
+ }
+}
+
+
+void matrix_in_place_transpose(uint8_t m, matrix_t matrix[][m])
+{
+ uint8_t i, j;
+ matrix_t tmp;
+
+ for (i = 0; i < m; i++) {
+ for (j = i + 1; j < m; j++) {
+ tmp = matrix[i][j];
+ matrix[i][j] = matrix[j][i];
+ matrix[j][i] = tmp;
+ }
+ }
+}
+
+void matrix_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
+ matrix_t dest_matrix[m][n])
+{
+ memcpy(dest_matrix, src_matrix, sizeof(matrix_t) * m * n);
+}
+
+void matrix_part_copy(uint8_t m, uint8_t n, matrix_t src_matrix[m][n],
+ uint8_t start_row_ind, uint8_t end_row_ind,
+ uint8_t start_col_ind, uint8_t end_col_ind,
+ uint8_t dest_row_num, uint8_t dest_col_num,
+ matrix_t dest_matrix[][dest_col_num])
+{
+ int16_t i, j;
+
+ if ((end_row_ind > start_row_ind) && (end_col_ind > start_col_ind)) {
+ for (i = start_row_ind;
+ (i <= end_row_ind)
+ && ((i - start_row_ind)
+ < dest_row_num);
+ i++) {
+ for (j = start_col_ind;
+ (j <= end_col_ind)
+ && ((j - start_col_ind)
+ < dest_col_num);
+ j++) {
+ dest_matrix[i - start_row_ind][j - start_col_ind] =
+ src_matrix[i][j];
+ }
+ }
+ }
+ else {
+ //copy one element
+ if ((start_row_ind == end_row_ind)
+ && (start_col_ind == end_col_ind)) {
+
+ dest_matrix[start_row_ind][start_row_ind] =
+ src_matrix[start_row_ind][start_row_ind];
+ }
+ //copy a row vector of the matrix
+ else if (start_row_ind == end_row_ind) {
+ i = 0;
+ for (j = start_col_ind; j <= end_col_ind; j++) {
+ dest_matrix[start_row_ind][i++] =
+ src_matrix[start_row_ind][j];
+ }
+ }
+ //copy a column vector of the matrix
+ else if (start_col_ind == end_col_ind) {
+ j = 0;
+ for (i = start_row_ind; i <= end_row_ind; i++) {
+ dest_matrix[j++][start_col_ind] =
+ src_matrix[i][start_col_ind];
+ }
+ }
+ }
+}
+
+
+void matrix_print(uint8_t m, uint8_t n, matrix_t matrix[m][n])
+{
+ int i, j;
+
+ puts("{");
+ for (i = 0; i < m; i++) {
+ printf("{");
+ for (j = 0; j < n; j++) {
+ printf("%7.4f", matrix[i][j]);
+ if (j < n - 1) {
+ printf(", ");
+ }
+ }
+ printf("}");
+ if (i < m - 1) {
+ printf(", ");
+ }
+ puts("");
+ }
+ printf("}");
+}
+
+
+void matrix_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t start_row_ind, uint8_t end_row_ind,
+ uint8_t start_col_ind, uint8_t end_col_ind)
+{
+ int16_t i, j;
+
+ if ((end_row_ind > start_row_ind) && (end_col_ind > start_col_ind)) {
+ puts("{");
+ for (i = start_row_ind; i <= end_row_ind; i++) {
+ printf("{");
+ for (j = start_col_ind; j <= end_col_ind; j++) {
+ printf("%7.4f", matrix[i][j]);
+ if (j < end_col_ind) {
+ printf(", ");
+ }
+ }
+ printf("}");
+ if (i < end_row_ind) {
+ printf(", ");
+ }
+ puts("");
+ }
+ printf("}");
+ }
+ else {
+ //print only one element of the matrix
+ if ((start_row_ind == end_row_ind)
+ & (start_col_ind == end_col_ind)) {
+
+ printf("{%3.4f}", matrix[start_row_ind][start_col_ind]);
+ }
+ // print row vector of the matrix
+ else if (start_row_ind == end_row_ind) {
+ printf("{");
+ for (j = start_col_ind; j <= end_col_ind; j++) {
+ printf("%3.4f", matrix[start_row_ind][j]);
+ if (j < end_col_ind) {
+ printf(", ");
+ }
+ }
+ printf("}");
+ }
+ // print column vector of the matrix
+ else if (start_col_ind == end_col_ind) {
+ puts("{");
+ for (i = start_row_ind; i <= end_row_ind; i++) {
+ printf("%3.4f", matrix[i][start_col_ind]);
+ if (i < end_row_ind) {
+ printf(",\n");
+ }
+ }
+ printf("}");
+ }
+ }
+}
+
+void matrix_flex_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t before_dot, uint8_t after_dot)
+{
+ uint16_t i, j;
+ char format_str_buff[13];
+
+ sprintf(format_str_buff, "%%%u.%uf", before_dot, after_dot);
+
+ puts("{");
+ for (i = 0; i < m; i++) {
+ printf("%11s", "{");
+ for (j = 0; j < n; j++) {
+ utils_printf(format_str_buff, matrix[i][j]);
+ if (j < n - 1) {
+ printf(", ");
+ }
+ }
+ printf("}");
+ if (i < m - 1) {
+ printf(", ");
+ }
+ puts("");
+ }
+ printf("%7s\n", "};");
+}
+
+
+void matrix_flex_part_print(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t start_row_ind, uint8_t end_row_ind,
+ uint8_t start_col_ind, uint8_t end_col_ind,
+ uint8_t before_dot, uint8_t after_dot)
+{
+ int16_t i, j;
+ char format_str_buff[13];
+
+ sprintf(format_str_buff, "%%%u.%uf", before_dot, after_dot);
+
+ if ((end_row_ind > start_row_ind) && (end_col_ind > start_col_ind)) {
+ puts("{");
+ for (i = start_row_ind; i <= end_row_ind; i++) {
+ printf("%11s", "{");
+ for (j = start_col_ind; j <= end_col_ind; j++) {
+ utils_printf(format_str_buff, matrix[i][j]);
+ if (j < end_col_ind) {
+ printf(", ");
+ }
+ }
+ printf("}");
+ if (i < end_row_ind) {
+ printf(", ");
+ }
+ puts("");
+ }
+ printf("%11s\n", "};");
+ }
+ else {
+ //print only one element of the matrix
+ if ((start_row_ind == end_row_ind)
+ & (start_col_ind == end_col_ind)) {
+ printf("{");
+ utils_printf(format_str_buff,
+ matrix[start_row_ind][start_col_ind]);
+ printf("}");
+ }
+ // print row vector of the matrix
+ else if (start_row_ind == end_row_ind) {
+ printf("{");
+ for (j = start_col_ind; j <= end_col_ind; j++) {
+ utils_printf(format_str_buff,
+ matrix[start_row_ind][j]);
+ if (j < end_col_ind) {
+ printf(", ");
+ }
+ }
+ printf("}");
+ }
+ // print column vector of the matrix
+ else if (start_col_ind == end_col_ind) {
+ puts("{");
+ for (i = start_row_ind; i <= end_row_ind; i++) {
+ utils_printf(format_str_buff,
+ matrix[i][start_col_ind]);
+ if (i < end_row_ind) {
+ printf(",\n");
+ }
+ }
+ printf("}");
+ }
+ }
+}
+
+
+uint8_t matrix_get_rank(uint8_t m, uint8_t n, matrix_t singl_values_arr[],
+ uint8_t length)
+{
+ uint8_t rank = 0;
+ int i;
+ double eps = pow(2.0, -52.0);
+ double tol = max(m, n) * singl_values_arr[0] * eps;
+
+ for (i = 0; i < length; i++) {
+ if (singl_values_arr[i] > tol) {
+ rank++;
+ }
+ }
+
+ return rank;
+}
+
+void matrix_sub(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n],
+ matrix_t A_minus_B[m][n])
+{
+ uint8_t i, j;
+
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < n; j++) {
+ A_minus_B[i][j] = A[i][j] - B[i][j];
+ }
+ }
+}
+
+void matrix_add(uint8_t m, uint8_t n, matrix_t A[m][n], matrix_t B[m][n],
+ matrix_t A_plus_B[m][n])
+{
+ uint8_t i, j;
+
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < n; j++) {
+ A_plus_B[i][j] = A[i][j] + B[i][j];
+ }
+ }
+}
+
+void matrix_add_to_diag(uint8_t n, matrix_t A[][n], uint8_t diag_el_num,
+ matrix_t value)
+{
+ uint8_t i;
+
+ for (i = 0; i < diag_el_num; i++) {
+ A[i][i] += value;
+ }
+}
+
+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])
+{
+
+ int i, j, k;
+
+ if (a_col_num != b_line_num) {
+ puts(
+ "the first matrix column number must be equal with the line number"
+ " of the second matrix !!!");
+ }
+ else {
+
+ for (i = 0; i < a_line_num; i++) {
+ for (j = 0; j < b_col_num; j++) {
+ dest_matrix[i][j] = 0;
+ for (k = 0; k < b_line_num; k++) {
+ dest_matrix[i][j] += a_matrix[i][k]
+ * b_matrix[k][j];
+ }
+ }
+ }
+ }
+}
+
+void matrix_part_mul(uint8_t a_col_num_max, matrix_t a_matrix[][a_col_num_max],
+ uint8_t b_col_num_max, matrix_t b_matrix[][b_col_num_max],
+ uint8_t a_start_row_ind, uint8_t a_end_row_ind,
+ uint8_t a_start_col_ind, uint8_t a_end_col_ind,
+ uint8_t b_start_row_ind, uint8_t b_end_row_ind,
+ uint8_t b_start_col_ind, uint8_t b_end_col_ind,
+ uint8_t dest_col_size,
+ matrix_t dest_matrix[][dest_col_size])
+{
+ uint8_t a_line_num;
+ uint8_t b_col_num;
+ uint8_t b_line_num;
+
+ if ((a_end_col_ind - a_start_col_ind)
+ != (b_end_row_ind - b_start_row_ind)) {
+ puts(
+ "the first matrix column number must be equal with the line number"
+ " of the second matrix !!!");
+ }
+ else {
+
+ a_line_num = a_end_row_ind - a_start_row_ind + 1;
+ b_line_num = b_end_row_ind - b_start_row_ind + 1;
+ b_col_num = b_end_col_ind - b_start_col_ind + 1;
+
+ for (uint8_t i = 0; i < a_line_num; i++) {
+ for (uint8_t j = 0; j < b_col_num; j++) {
+ dest_matrix[i][j] = 0;
+ for (uint8_t k = 0; k < b_line_num; k++) {
+ dest_matrix[i][j] +=
+ a_matrix[i
+ + a_start_row_ind][k
+ + a_start_col_ind]
+ * b_matrix[k
+ + b_start_row_ind][j
+ + b_start_col_ind];
+ }
+ }
+ }
+ }
+}
+
+
+void matrix_mul_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ matrix_t vec[n], matrix_t dst_arr[m])
+{
+
+ uint8_t i, j;
+
+ if (vec != dst_arr) {
+ memset(dst_arr, 0, m * sizeof(matrix_t));
+ }
+
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < n; j++) {
+ dst_arr[i] += matrix[i][j] * vec[j];
+ }
+ }
+}
+
+
+void matrix_vec_mul_matr(uint8_t m, uint8_t n, matrix_t vec[m],
+ matrix_t matrix[m][n], matrix_t dst_arr[n])
+{
+ uint8_t i, j;
+
+ //Make sense ?
+ if (vec != dst_arr) {
+ memset(dst_arr, 0, n * sizeof(matrix_t));
+ }
+
+ for (i = 0; i < n; i++) { //column
+ for (j = 0; j < m; j++) { //rows
+ dst_arr[i] += vec[j] * matrix[j][i];
+ }
+ }
+}
+
+
+void matrix_mul_scalar_vec_matr(uint8_t m, uint8_t n, matrix_t scalar,
+ matrix_t vec[m], matrix_t matrix[m][n],
+ matrix_t dst_arr[n])
+{
+ uint8_t i, j;
+
+ //Make sense ?
+ if (vec != dst_arr) {
+ memset(dst_arr, 0, n * sizeof(matrix_t));
+ }
+
+ for (i = 0; i < n; i++) { //column
+ for (j = 0; j < m; j++) { //rows
+ dst_arr[i] += scalar * vec[j] * matrix[j][i];
+ }
+ }
+}
+
+
+void matrix_part_mul_scalar_vec_matr(uint8_t max_m, uint8_t max_n,
+ matrix_t scalar, matrix_t vec[max_m],
+ matrix_t matrix[max_m][max_n],
+ uint8_t begin_row, uint8_t begin_column,
+ matrix_t dst_arr[max_n - begin_row])
+{
+ int16_t i, j;
+
+ if (vec != dst_arr) {
+ memset(dst_arr, 0, (max_n - begin_row) * sizeof(matrix_t));
+ }
+
+ for (i = begin_column; i < max_n; i++) { //column
+ for (j = begin_row; j < max_m; j++) { //rows
+ dst_arr[i - begin_column] += scalar * vec[j - begin_row]
+ * matrix[j][i];
+
+ }
+ }
+}
+
+
+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])
+{
+ uint8_t i, j;
+
+ if (m != b_size) {
+ puts(
+ "The vector size should be equal the raw number of the matrix !!!");
+ }
+ else {
+ for (j = 0; j < n; j++) {
+ c_vec[j] = 0;
+ for (i = 0; i < m; i++) {
+ c_vec[j] += A[i][j] * b_vec[i];
+ }
+ }
+ }
+}
+
+
+void matrix_mul_col_vec_row_vec(uint8_t m, matrix_t col_vec[m], uint8_t n,
+ matrix_t row_vec[n], uint8_t max_n,
+ matrix_t res_mat[][max_n])
+{
+ uint8_t i, j;
+
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < n; j++) {
+ res_mat[i][j] = col_vec[i] * row_vec[j];
+ }
+ }
+}
+
+
+void matrix_trans_mul_itself(uint8_t m, uint8_t n, matrix_t A[m][n],
+ matrix_t AT_mul_A[n][n])
+{
+ uint8_t i, j, k;
+ matrix_t column_vec[m];
+ matrix_t tmp_col_vec[m];
+
+ for (i = 0; i < n; i++) {
+ matrix_get_column_vec(m, n, A, i, column_vec);
+ for (j = 0; j < n; j++) {
+ matrix_get_column_vec(m, n, A, j, tmp_col_vec);
+ AT_mul_A[i][j] = 0;
+ for (k = 0; k < m; k++) {
+ AT_mul_A[i][j] += column_vec[k]
+ * tmp_col_vec[k];
+ }
+ }
+ }
+}
+
+void matrix_set_diag_elements(uint8_t m, uint8_t n, matrix_t value,
+ matrix_t diag_matrix[m][n])
+{
+ uint8_t max = 0;
+
+ if (m < n) {
+ max = m;
+ }
+ else {
+ max = n;
+ }
+ for (uint8_t i = 0; i < max; i++) {
+ diag_matrix[i][i] = value;
+ }
+}
+
+void matrix_get_diag_mat_new(uint8_t m, uint8_t n, matrix_t diag_matrix[m][n],
+ uint8_t length, matrix_t vec[])
+{
+ uint8_t i;
+
+ matrix_clear(m, n, diag_matrix);
+
+ for (i = 0; i < length; i++) {
+ diag_matrix[i][i] = vec[i];
+ }
+}
+
+
+void matrix_get_diag_mat(uint8_t m, uint8_t n, matrix_t value,
+ matrix_t diag_matrix[m][n])
+{
+ matrix_clear(m, n, diag_matrix);
+ matrix_set_diag_elements(m, n, value, diag_matrix);
+}
+
+
+void matrix_mul_scalar(uint8_t m, uint8_t n, matrix_t mat_src[m][n],
+ matrix_t value, matrix_t mat_dest[m][n])
+{
+ for (int i = 0; i < m; i++) {
+ for (int j = 0; j < n; j++) {
+ mat_dest[i][j] = value * mat_src[i][j];
+ }
+ }
+}
+
+void matrix_get_column_vec(uint8_t m, uint8_t n, matrix_t matrix[m][n],
+ uint8_t col_num, matrix_t col_vec[m])
+{
+ uint8_t i;
+
+ for (i = 0; i < m; i++) {
+ col_vec[i] = matrix[i][col_num];
+ }
+}
+
+
+void matrix_get_part_column_vec(uint8_t max_m, uint8_t max_n,
+ matrix_t matrix[max_m][max_n], uint8_t col_num,
+ uint8_t offset,
+ matrix_t col_vec[max_m - offset])
+{
+ uint8_t i;
+
+ for (i = offset; i < max_m; i++) {
+ col_vec[i - offset] = matrix[i][col_num];
+ }
+}
+
+matrix_t matrix_get_max_elem_in_column(uint8_t m, uint8_t n,
+ matrix_t matrix[m][n], uint8_t col_num)
+{
+ matrix_t elem_max;
+ uint8_t i;
+
+ elem_max = matrix[0][col_num];
+ for (i = 1; i < m; i++) {
+ if (matrix[i][col_num] > elem_max) {
+ elem_max = matrix[i][col_num];
+ }
+ }
+
+ return elem_max;
+}
+
+matrix_t matrix_get_abs_max_elem_in_column(uint8_t m, uint8_t n,
+ matrix_t matrix[m][n],
+ uint8_t col_num)
+{
+ matrix_t abs_max_elem;
+ uint8_t i;
+
+ abs_max_elem = fabs(matrix[0][col_num]);
+ for (i = 1; i < m; i++) {
+ if (fabs(matrix[i][col_num]) > abs_max_elem) {
+ abs_max_elem = fabs(matrix[i][col_num]);
+ }
+ }
+
+ return abs_max_elem;
+}
+
+matrix_t matrix_get_max_elem_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)
+{
+ matrix_t max_elem;
+ uint8_t i;
+
+ max_elem = matrix[row_num][col_num];
+ for (i = row_num + 1; i < max_m; i++) {
+ if (matrix[i][col_num] > max_elem) {
+ max_elem = matrix[i][col_num];
+ }
+ }
+
+ return max_elem;
+}
+
+matrix_t matrix_get_abs_max_elem_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)
+{
+ matrix_t abs_max_elem;
+ uint8_t i;
+
+ abs_max_elem = fabs(matrix[row_num][col_num]);
+ for (i = row_num + 1; i < max_m; i++) {
+ if (fabs(matrix[i][col_num]) > abs_max_elem) {
+ abs_max_elem = fabs(matrix[i][col_num]);
+ }
+ }
+
+ return abs_max_elem;
+}
+
+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)
+{
+ matrix_t abs_max_elem;
+ uint8_t i;
+
+ abs_max_elem = fabs(matrix[row_num][col_num]);
+ *index = row_num;
+ for (i = row_num + 1; i < max_m; i++) {
+ if (fabs(matrix[i][col_num]) > abs_max_elem) {
+ abs_max_elem = fabs(matrix[i][col_num]);
+ *index = i;
+ }
+
+ }
+
+ return abs_max_elem;
+}
+
+void matrix_swap_rows(uint8_t n, matrix_t matrix[][n], uint8_t i, uint8_t j)
+{
+ matrix_t tmp;
+
+ for (uint8_t k = 0; k < n; k++) {
+ tmp = matrix[i][k];
+ matrix[i][k] = matrix[j][k];
+ matrix[j][k] = tmp;
+ }
+}
+
+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)
+{
+ matrix_t tmp;
+
+ if (col_begin <= col_end) {
+ for (uint8_t k = col_begin; k < (col_end + 1); k++) {
+ tmp = matrix[i][k];
+ matrix[i][k] = matrix[j][k];
+ matrix[j][k] = tmp;
+ }
+ }
+
+}
+
+double matrix_get_two_norm(uint8_t m, uint8_t n, matrix_t A[][n])
+{
+ double two_norm = 0.0;
+ matrix_dim_t dim;
+ matrix_t copy_A[m][n];
+
+ svd_get_U_dim(m, n, &dim);
+ matrix_t U[dim.row_num][dim.col_num];
+ matrix_t S[n][n];
+ matrix_t V[n][n];
+ uint8_t s_length = svd_get_single_values_num(m, n);
+ matrix_t s[s_length];
+
+ matrix_copy(m, n, A, copy_A);
+ svd(m, n, copy_A, dim.row_num, dim.col_num, U, S, V, s_length, s);
+ printf("U1 = ");
+ matrix_print(dim.row_num, dim.col_num, U);
+ puts("");
+ printf("S1 = ");
+ matrix_print(n, n, S);
+ puts("");
+ printf("V1 = ");
+ matrix_print(n, n, V);
+ puts("");
+
+ two_norm = S[0][0];
+
+ return two_norm;
+}
+
+double matrix_get_frob_norm(uint8_t m, uint8_t n, matrix_t A[][n])
+{
+ double frob_norm = 0.0;
+
+ for (uint8_t i = 0; i < m; i++) {
+ for (uint8_t j = 0; j < n; j++) {
+ frob_norm += pow(A[i][j], 2);
+ }
+ }
+
+ return sqrt(frob_norm);
+}
+
+void matrix_get_inv_upp_triang(uint8_t m, uint8_t n, matrix_t U[][n],
+ matrix_t inv_U[][m])
+{
+ //Zero out the inv_U matrix
+ matrix_clear(n, m, inv_U);
+ for (uint8_t i = 0; i < n; i++) {
+ if (U[i][i] != 0) {
+ inv_U[i][i] = 1 / U[i][i];
+ }
+ else {
+ inv_U[i][i] = FLT_MAX;
+ }
+ for (uint8_t j = 0; j <= (i - 1); j++) {
+ matrix_t u = 0.0;
+ for (uint8_t k = j; k <= (i - 1); k++) {
+ u = u + U[k][i] * inv_U[j][k];
+ }
+ inv_U[j][i] = -u * inv_U[i][i];
+ }
+ }
+}
+
+void matrix_get_inv_low_triang(uint8_t m, uint8_t n, matrix_t L[][n],
+ matrix_t inv_L[][m])
+{
+ //Zero out the inv_L matrix
+ matrix_clear(n, m, inv_L);
+ for (uint8_t i = 0; i < n; i++) {
+ if ((i <= m) & (i <= n)) {
+ if (L[i][i] != 0) {
+ inv_L[i][i] = 1 / L[i][i];
+ }
+ else {
+ inv_L[i][i] = FLT_MAX;
+ }
+ for (uint8_t j = 0; j <= (i - 1); j++) {
+ matrix_t l = 0.0;
+ for (uint8_t k = j; k <= (i - 1); k++) {
+ l = l + L[i][k] * inv_L[k][j];
+ }
+ inv_L[i][j] = -l * inv_L[i][i];
+ }
+ } //if
+ } //for
+}
+
+void matrix_get_upp_triang(uint8_t m, uint8_t n, matrix_t A[][n],
+ matrix_t tr_up_A[][n])
+{
+ // Same matrix
+ if ((&A[0][0]) == (&tr_up_A[0][0])) {
+ for (uint8_t i = 1; i < m; i++) {
+ for (uint8_t j = 0; (j < i) && (j < n); j++) {
+ A[i][j] = 0;
+ }
+ }
+
+ }
+ else {
+ for (uint8_t i = 0; i < m; i++) {
+ for (uint8_t j = 0; j < n; j++) {
+ if (i <= j) {
+ tr_up_A[i][j] = A[i][j];
+ }
+ else {
+ tr_up_A[i][j] = 0;
+ }
+ }
+ }
+ }
+
+}
+
+void matrix_get_low_triang(uint8_t m, uint8_t n, matrix_t A[][n],
+ matrix_t tr_low_A[][n])
+{
+
+ if ((&A[0][0]) == (&tr_low_A[0][0])) {
+ // Same matrix
+ for (uint8_t i = 0; i < m; i++) {
+ for (uint8_t j = i + 1; (j > i) && (j < n); j++) {
+ A[i][j] = 0;
+ }
+ }
+
+ }
+ else {
+ for (uint8_t i = 0; i < m; i++) {
+ for (uint8_t j = 0; j < n; j++) {
+ if (i >= j) {
+ tr_low_A[i][j] = A[i][j];
+ }
+ else {
+ tr_low_A[i][j] = 0;
+ }
+ }
+ }
+ }
+}
+
+matrix_t matrix_read(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i,
+ uint8_t j)
+{
+ return matrix[i][j];
+}
+
+void matrix_write(uint8_t m, uint8_t n, matrix_t matrix[m][n], uint8_t i,
+ uint8_t j, matrix_t val)
+{
+ matrix[i][j] = val;
+}
diff --git a/linear_algebra/basic_operations/vector.c b/linear_algebra/basic_operations/vector.c
new file mode 100644
index 0000000000000000000000000000000000000000..b586c3a9f3c44f167495d86003185bc1a440ba4c
--- /dev/null
+++ b/linear_algebra/basic_operations/vector.c
@@ -0,0 +1,301 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Vector computations.
+ * Vector computations include operations such as addition,
+ * subtraction, and inner product (dot product).
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <inttypes.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdbool.h>
+
+#include "vector.h"
+#include "utils.h"
+
+void vector_clear(uint8_t n, vector_t arr[])
+{
+ memset(arr, 0, sizeof(vector_t) * n);
+}
+
+void vector_copy(uint8_t size, vector_t src_arr[], vector_t dest_arr[])
+{
+ memcpy(dest_arr, src_arr, size * sizeof(vector_t));
+}
+
+vector_t vector_get_norm2(uint8_t length, vector_t arr[])
+{
+ vector_t square_norm = 0.0;
+ uint8_t i = 0;
+
+ for (; i < length; i++) {
+ square_norm += arr[i] * arr[i];
+ }
+
+ return sqrt(square_norm);
+}
+
+vector_t vector_get_square_norm2(uint8_t length, vector_t arr[])
+{
+ vector_t square_norm = 0.0;
+ uint8_t i = 0;
+
+ for (; i < length; i++) {
+ square_norm += arr[i] * arr[i];
+ }
+
+ return square_norm;
+}
+
+vector_t vector_get_sum(uint8_t length, vector_t arr[])
+{
+ vector_t sum = 0.0;
+ uint8_t i;
+
+ for (i = 0; i < length; i++) {
+ sum += arr[i];
+ }
+
+ return sum;
+}
+
+vector_t vector_get_mean_value(uint8_t length, vector_t arr[])
+{
+ vector_t mean = 0.0;
+ uint8_t i;
+
+ for (i = 0; i < length; i++) {
+ mean += arr[i];
+ }
+
+ if (length != 0) {
+ mean /= length;
+ }
+
+ return mean;
+}
+
+void vector_sub(uint8_t size, vector_t a_vec[], vector_t b_vec[],
+ vector_t a_minus_b[])
+{
+ uint8_t i;
+
+ for (i = 0; i < size; i++) {
+ a_minus_b[i] = a_vec[i] - b_vec[i];
+ }
+}
+
+void vector_add(uint8_t size, vector_t a_vec[size], vector_t b_vec[size],
+ vector_t a_plus_b_vec[size])
+{
+ uint8_t i;
+
+ for (i = 0; i < size; i++) {
+ a_plus_b_vec[i] = a_vec[i] + b_vec[i];
+ }
+}
+
+void vector_mul(uint8_t size, vector_t a_vec[size], vector_t b_vec[size],
+ vector_t a_mul_b_vec[size])
+{
+ uint8_t i;
+
+ for (i = 0; i < size; i++) {
+ a_mul_b_vec[i] = a_vec[i] * b_vec[i];
+ }
+}
+
+void vector_square(uint8_t n, vector_t vec[n], vector_t square_vec[n])
+{
+ for (int i = 0; i < n; i++) {
+ square_vec[i] = vec[i] * vec[i];
+ }
+}
+
+void vector_in_place_scalar_mul(uint8_t size, vector_t a_vec[size],
+ vector_t scl)
+{
+ uint8_t i;
+
+ for (i = 0; i < size; i++) {
+ a_vec[i] = scl * a_vec[i];
+ }
+}
+
+void vector_scalar_mul(uint8_t size, vector_t src_vec[size], vector_t scl,
+ vector_t dest_vec[])
+{
+ uint8_t i;
+
+ for (i = 0; i < size; i++) {
+ dest_vec[i] = scl * src_vec[i];
+ }
+}
+
+void vector_scalar_div(uint8_t size, vector_t a_vec[size], vector_t scl)
+{
+ uint8_t i;
+
+ if (scl != 0) {
+ for (i = 0; i < size; i++) {
+ a_vec[i] = a_vec[i] / scl;
+ }
+ }
+}
+
+//Euclidean distance: d = sum((x-y).^2).^0.5
+vector_t vector_get_euclidean_distance(uint8_t length, vector_t vec1[],
+ vector_t vec2[])
+{
+ vector_t d = 0.0;
+ vector_t diff_vec[length];
+
+ vector_sub(length, vec1, vec2, diff_vec);
+ d = vector_get_norm2(length, diff_vec);
+
+ return d;
+}
+
+vector_t vector_get_max_and_index(uint8_t length, vector_t vec[],
+ uint8_t *index)
+{
+ vector_t max = vec[0];
+
+ *index = 0;
+ for (uint8_t i = 1; i < length; i++) {
+ if (vec[i] > max) {
+ max = vec[i];
+ *index = *index + 1;
+ }
+ }
+ return max;
+}
+
+vector_t vector_get_scalar_product(uint8_t n, vector_t vec1[n],
+ vector_t vec2[n])
+{
+ vector_t r = 0;
+
+ for (int i = 0; i < n; i++) {
+ r = r + (vec1[i] * vec2[i]);
+ }
+
+ return r;
+}
+
+bool vector_is_equal(uint16_t length, vector_t vec_1[], vector_t vec_2[])
+{
+
+ for (uint16_t i = 0; i < length; i++) {
+ if (vec_1[i] != vec_2[i]) {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+void vector_get_index_vector(uint8_t k, uint8_t n, vector_t unsorted_vector[n],
+ vector_t sorted_vector[n], uint8_t index_vector[n])
+{
+
+ for (int i = 0; i < k; i++) {
+ for (int j = 0; j < n; j++) {
+ if (unsorted_vector[j] == sorted_vector[i]) {
+ index_vector[i] = j;
+ }
+ }
+ }
+}
+
+vector_t vector_get_residual(uint8_t length, vector_t a_vec[], vector_t b_vec[])
+{
+ vector_t diff_vec[length];
+
+ vector_sub(length, a_vec, b_vec, diff_vec);
+ return vector_get_norm2(length, diff_vec);
+}
+
+void vector_get_elements(vector_t src_vec[], uint8_t k, uint8_t index_vec[],
+ vector_t dst_vec[])
+{
+ for (uint8_t i = 0; i < k; i++) {
+ dst_vec[i] = src_vec[index_vec[i]];
+ }
+}
+
+bool vector_uint32_is_equal(uint32_t length, uint32_t vec_1[], uint32_t vec_2[])
+{
+
+ for (uint32_t i = 0; i < length; i++) {
+ if (vec_1[i] != vec_2[i]) {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+void vector_print(uint32_t length, vector_t arr[])
+{
+ uint16_t i;
+
+ printf("{");
+ for (i = 0; i < length; i++) {
+ printf("%5.4f", arr[i]);
+ if (i < length - 1) {
+ printf(", ");
+ }
+ }
+ printf("}");
+}
+
+void vector_print_u8_array(uint32_t length, uint8_t arr[])
+{
+ uint32_t i;
+
+ printf("{");
+ for (i = 0; i < length; i++) {
+ printf("%u", arr[i]);
+ if (i < length - 1) {
+ printf(", ");
+ }
+ }
+ printf("}");
+}
+
+// This function is more memory-consuming than vector_print
+void vector_flex_print(uint32_t length, vector_t arr[], uint8_t before_dot,
+ uint8_t after_dot)
+{
+ uint32_t i;
+ char format_str_buff[13];
+
+ sprintf(format_str_buff, "%%%u.%uf", before_dot, after_dot);
+ printf("{");
+ for (i = 0; i < length; i++) {
+ //printf(format_str_buff, arr[i]);
+ utils_printf(format_str_buff, arr[i]);
+ if (i < length - 1) {
+ printf(", ");
+ }
+ }
+ printf("}");
+ //puts("");
+}
diff --git a/linear_algebra/doc.txt b/linear_algebra/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..91081ada6318963d5b3cd4222db5c63f294372aa
--- /dev/null
+++ b/linear_algebra/doc.txt
@@ -0,0 +1,21 @@
+/*
+ * Copyright (C) 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup linear_algebra LINEAR_ALGEBRA
+ *
+ * @brief Linear algebra operations
+ *
+ * @details The linear algebra module contains functions that are specific to vector and
+ * matrix operations, and other algebraic operations. This module is composed of
+ * five submodules: basic_operation, matrix_decompositions, pseudo_inverse,
+ * solve_linear_equations, and utilities submodules.
+ *
+ * @author Zakaria Kasmi
+ *
+ */
diff --git a/linear_algebra/matrix_decompositions/doc.txt b/linear_algebra/matrix_decompositions/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d8bee1408e32a491747feca02e5e99404dfe23f2
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/doc.txt
@@ -0,0 +1,20 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup matrix_decompositions MATRIX_DECOMPOSITIONS
+ * @ingroup linear_algebra
+ *
+ * @brief Matrix decomposition algorithms.
+ *
+ * @details This module provides algorithms of matrix decomposition such as
+ * Householder, Givens, or Golub-Kahan-Reinsch algorithms.
+ *
+ * @author Zakaria Kasmi
+ */
diff --git a/linear_algebra/matrix_decompositions/include/lu_decomp.h b/linear_algebra/matrix_decompositions/include/lu_decomp.h
new file mode 100644
index 0000000000000000000000000000000000000000..59f955c7f8b3b49c8762813f022a4bf4ff2d662d
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/include/lu_decomp.h
@@ -0,0 +1,46 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup matrix_decompositions
+ * @{
+ *
+ * @file
+ * @brief Computes the LU decomposition of the matrix.
+ * @details Computes the permutation matrix P such that: A = P'*L*U, where L is a lower
+ * triangular matrix and U is an upper triangular matrix. It implements the Gaussian
+ * Elimination (GE) with pivoting algorithm.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+#ifndef LU_DECOMP_H_
+#define LU_DECOMP_H_
+
+#include "matrix.h"
+
+/**
+ * @brief Computes the LU decomposition of the matrix.
+ * @details Computes the permutation matrix P such that: A = P'*L*U, where L is a lower triangular
+ * matrix and U is an upper triangular matrix. It implements the Gaussian Elimination with
+ * pivoting algorithm.
+ *
+ * @note Matrix U is stored in the matrix A.
+ *
+ * @param[in] n column number of the matrix.
+ * @param[in,out] A[][] pointer to the matrices A and U.
+ * @param[out] L[][] pointer to the L matrix.
+ * @param[out] P[][] pointer to the P matrix.
+ *
+ * @return the number of changes by computing the LU decomposition.
+ */
+uint8_t lu_decomp(uint8_t n, matrix_t A[][n], matrix_t L[][n], matrix_t P[][n]);
+
+#endif /*LU_DECOMP_H_ */
diff --git a/linear_algebra/matrix_decompositions/include/qr_common.h b/linear_algebra/matrix_decompositions/include/qr_common.h
new file mode 100644
index 0000000000000000000000000000000000000000..3edc335d75739849846ff5351d21baddafa239d4
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/include/qr_common.h
@@ -0,0 +1,62 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup matrix_decompositions
+ * @{
+ *
+ * @file
+ * @brief Common definitions and implementations for the QR-decomposition.
+ * Provide necessary methods to construct Q- and R- matrices using.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef QR_COMMON_H_
+#define QR_COMMON_H_
+
+#include <inttypes.h>
+#include "matrix.h"
+
+/**
+ * Possible algorithms to compute the QR-decomposition of a matrix.
+ */
+enum QR_ALGORITHM {
+ QR_Householder, QR_Givens
+};
+
+/**
+ * @brief Implements the backward substitution algorithm.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] U[][] pointer to the matrix U.
+ * @param[in] b[] pointer to the vector b.
+ * @param[out] x_sol[] pointer to the solution of the substitution algorithm.
+ *
+ */
+void qr_common_backward_subst(uint8_t m, uint8_t n, matrix_t U[][n],
+ matrix_t b[m], matrix_t x_sol[m]);
+/**
+ * @brief Compute the reduced form of the QR-decomposition algorithm.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] Q[][] pointer to the matrix Q.
+ * @param[in] R[][] pointer to the matrix R.
+ * @param[out] red_Q[][] pointer to the reduced matrix Q.
+ * @param[out] red_R[][] pointer to the reduced matrix R.
+ *
+ */
+void qr_common_get_reduced_QR(uint8_t m, uint8_t n, matrix_t Q[m][m],
+ matrix_t R[m][n], matrix_t red_Q[m][n], matrix_t red_R[n][n]);
+
+#endif /* QR_COMMON_H_ */
diff --git a/linear_algebra/matrix_decompositions/include/qr_givens.h b/linear_algebra/matrix_decompositions/include/qr_givens.h
new file mode 100644
index 0000000000000000000000000000000000000000..1708c332f3bfda8d1aaf6d3678008b1508ae4173
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/include/qr_givens.h
@@ -0,0 +1,71 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup matrix_decompositions
+ * @{
+ *
+ * @file
+ * @brief Givens algorithm for the QR-decomposition.
+ * Provide necessary methods to construct Q- and R- matrices using
+ * Givens rotations. A = QR, where Q is an (m \f$\times\f$ n)-matrix with
+ * orthonormal columns and R is an (n \f$\times\f$ n) upper triangular matrix.
+ *
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef QR_GIVENS_H_
+#define QR_GIVENS_H_
+
+#include <inttypes.h>
+
+#include "matrix.h"
+
+/**
+ * @brief Computes the QR decomposition of the matrix A by using the Givens
+ * algorithm.
+ * @details Gets a QR decomposition of an m-by-n matrix A such that A = Q*R.
+ * The compact as well as the full decomposition of the matrix can be computed.
+ *
+ * @note R is stored in the matrix A.
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in,out] A[][] pointer to the matrices A and R.
+ * @param[in] q_col_num column number of the matrix Q.
+ * @param[out] Q[][] pointer to the Q matrix.
+ * @param[in] reduced computes the compact form of the QR decomposition if true,
+ * otherwise the full version.
+ *
+ * @return 1, if computing the QR decomposition is successful.
+ * @return -1, if computing the QR decomposition is not successful.
+ *
+ */
+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);
+
+/**
+ * @brief Compute the Givens parameters.
+ * The computation of the parameters c, s, t, and r can have problems with
+ * overflow or underflow, therefore this algorithm employs a
+ * normalization procedure.
+ * The Givens parameters c, s, t and r are saved in a vector.
+ *
+ * @param[in] xjj value at the diagonal j of the matrix.
+ * @param[in] xij value at the index j of a column vector i.
+ * @param[out] c_s_t_r_vec[] pointer to the vector holding the c, s, t and r parameters.
+ *
+ */
+void qr_givens_get_params(matrix_t xjj, matrix_t xij, matrix_t c_s_t_r_vec[]);
+
+#endif /* QR_GIVENS_H_ */
diff --git a/linear_algebra/matrix_decompositions/include/qr_householder.h b/linear_algebra/matrix_decompositions/include/qr_householder.h
new file mode 100644
index 0000000000000000000000000000000000000000..97e9707586aa9850beae315f831bf9ec4cd7043e
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/include/qr_householder.h
@@ -0,0 +1,52 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup matrix_decompositions
+ * @{
+ *
+ * @file
+ * @brief Householder algorithm for the QR-decomposition.
+ * @details Provide necessary methods to construct Q- and R- matrices using
+ * Householder reflections. A = QR, where Q is an (m \f$\times\f$ n)-matrix with
+ * orthonormal columns and R is an (n \f$\times\f$ n) upper triangular matrix.
+ *
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+#ifndef QR_HOUSEHOLDER_H_
+#define QR_HOUSEHOLDER_H_
+
+#include <inttypes.h>
+
+#include "matrix.h"
+
+/**
+ * @brief Computes the QR decomposition of the matrix A by using the
+ * Householder algorithm.
+ * @note R is stored in the matrix A.
+ *
+ * @param[in] m row number of the matrix to decompose in QR.
+ * @param[in] n column number of the matrix to decompose in QR.
+ * @param[in,out] A[][] pointer to the matrices A and R.
+ * @param[in] q_col_num column number of the matrix Q.
+ * @param[in,out] Q[][] pointer to the matrix Q.
+ * @param[in] reduced computes the compact form of the QR decomposition if true,
+ * otherwise the full version.
+ *
+ * @return 1, if computing the QR decomposition is successful.
+ * @return -1, if computing the QR decomposition is not successful.
+ *
+ */
+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);
+#endif /* QR_HOUSEHOLDER_H_ */
diff --git a/linear_algebra/matrix_decompositions/include/svd.h b/linear_algebra/matrix_decompositions/include/svd.h
new file mode 100644
index 0000000000000000000000000000000000000000..109bdce87e6e691b1be9ff052856beb3597db695
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/include/svd.h
@@ -0,0 +1,150 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup matrix_decompositions
+ * @{
+ *
+ * @file
+ * @brief Algorithm for the Singular Value Decomposition (SVD).
+ * @details Provide necessary methods to compute the compact SVD of
+ * a matrix. A = U*S*V, where U is a (m x l) orthogonal matrix,
+ * S is a (l x l) diagonal matrix, V is a (l x n) orthogonal
+ * matrix, and l = min(m,n). The SVD is computed by using the
+ * Golub--Kahan--Reinsch algorithm that works in two phases:
+ * bidiagonalization and a reduction to the diagonal form phase.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef SVD_H_
+#define SVD_H_
+
+#include "matrix.h"
+
+/* Define the cases of the Golub-Reinsch Algorithm */
+
+/**
+ * The case of computing negligible values.
+ */
+#define SVD_COMPUTE_NEGLIGIBLE_VALUES 1
+
+/**
+ * The case of splitting at negligible values.
+ */
+#define SVD_SPLIT_AT_NEGLIGIBLE_VALUES 2
+
+/**
+ * The case of the QR-step.
+ */
+#define SVD_QR_STEP 3
+
+/**
+ * The case of the order of absolute singular values.
+ */
+#define SVD_ORDER_ABSOLUTE_SING_VALUES 4
+
+/**
+ * @brief Compute the Singular-Value Decomposition (SVD) of a matrix.
+ *
+ * @details Matrix A is transformed to A = U*S*V, where U and V are unitary
+ * matrices, and S is a diagonal matrix.
+ *
+ *
+ * @param[in] m row number of the matrix A.
+ * @param[in] n column number of the matrix A.
+ * @param[in, out] A[][] pointer to the matrix A.
+ * @param[in] u_m row number of the matrix U.
+ * @param[in] u_n column number of the matrix U.
+ * @param[in, out] U[][] pointer to the matrix U.
+ * @param[in, out] S[][] pointer to the matrix S.
+ * @param[in, out] V[][] pointer to the matrix V.
+ * @param[in] sing_vec_length length of the singular vector.
+ * @param[in, out] singl_values_vec[] pointer to the vector saving the singular values.
+ *
+ */
+void svd(uint8_t m, uint8_t n, matrix_t A[m][n],
+ uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n],
+ matrix_t S[u_n][n], matrix_t V[n][n],
+ uint8_t sing_vec_length, matrix_t singl_values_vec[]);
+
+/**
+ * @brief Calculate the dimension of the matrix U.
+ *
+ * @param[in] m row number of the matrix to decompose.
+ * @param[in] n column number of the matrix to decompose.
+ * @param[out] u_dim pointer to the dimension struct.
+ *
+ */
+void svd_get_U_dim(uint8_t m, uint8_t n, matrix_dim_t *u_dim);
+
+/**
+ * @brief Calculate the dimension of the matrix S.
+ *
+ * @param[in] m row number of the matrix to decompose.
+ * @param[in] n column number of the matrix to decompose.
+ * @param[out] s_dim pointer to the dimension struct.
+ *
+ */
+void svd_get_S_dim(uint8_t m, uint8_t n, matrix_dim_t *s_dim);
+
+/**
+ * @brief Calculate the dimension of the matrix V.
+ *
+ * @param[in] m row number of the matrix to decompose.
+ * @param[in] n column number of the matrix to decompose.
+ * @param[out] v_dim pointer to the dimension struct.
+ *
+ */
+void svd_get_V_dim(uint8_t m, uint8_t n, matrix_dim_t *v_dim);
+
+/**
+ * @brief Calculate the number of the singular values.
+ *
+ * @param[in] m row number of the matrix to decompose.
+ * @param[in] n column number of the matrix to decompose.
+ *
+ * @return the number of the singular values.
+ *
+ */
+uint8_t svd_get_single_values_num(uint8_t m, uint8_t n);
+
+/**
+ * @brief Compute the reciprocal singular values.
+ *
+ * @details This method is based on the singular values.
+ *
+ * @param[in] m row number of the matrix to transform in SVD.
+ * @param[in] n column number of the matrix to transform in SVD.
+ * @param[in] length length of the array of single values.
+ * @param[in] singl_values_arr pointer to the array of single values.
+ * @param[out] recip_singl_values_arr pointer to the array of the reciprocal
+ * singular values.
+ *
+ */
+void svd_get_reciproc_singular_values(uint8_t m, uint8_t n, uint8_t length,
+ matrix_t singl_values_arr[],
+ matrix_t recip_singl_values_arr[]
+ );
+/**
+ * @brief Compute and print the SVD of a matrix.
+ *
+ *
+ * @param[in] m row number of the matrix to transform in SVD.
+ * @param[in] n column number of the matrix to transform in SVD.
+ * @param[in] matrix_arr[][] pointer to the matrix.
+ * @param[in] i label.
+ *
+ */
+void svd_compute_print_U_S_V_s(uint8_t m, uint8_t n, matrix_t matrix_arr[m][n],
+ uint8_t i);
+
+#endif /* SVD_H_ */
diff --git a/linear_algebra/matrix_decompositions/lu_decomp.c b/linear_algebra/matrix_decompositions/lu_decomp.c
new file mode 100644
index 0000000000000000000000000000000000000000..9be7fc7580348c0cfcb94db6cbdd27ad7be32b39
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/lu_decomp.c
@@ -0,0 +1,85 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Computes the LU decomposition of the matrix.
+ * @details Computes the permutation matrix P such that: A = P'*L*U, where L is a lower
+ * triangular matrix and U is an upper triangular matrix. It implements the Gaussian
+ * Elimination (GE) with pivoting algorithm.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <float.h>
+#include <stdio.h>
+
+#include "matrix.h"
+#include "vector.h"
+
+// U will be saved in A
+uint8_t lu_decomp(uint8_t n, matrix_t A[][n], matrix_t L[][n], matrix_t P[][n])
+{
+ matrix_t abs_max_col_elem;
+ uint8_t k;
+ uint8_t pivot_index;
+ uint8_t changes;
+ matrix_t multipliers_vec[n - 1];
+
+ matrix_get_diag_mat(n, n, 1, L);
+ matrix_get_diag_mat(n, n, 1, P);
+ changes = 0;
+ for (uint8_t i = 0; i < n - 1; i++) {
+ abs_max_col_elem =
+ matrix_get_abs_max_elem_and_index_in_part_column(
+ n, n, A, i, i, &k);
+ if (abs_max_col_elem <= FLT_EPSILON) {
+ continue;
+ }
+
+ pivot_index = k;
+
+ if (pivot_index != i) {
+ //Swap the rows at the position i and pivot_index of the matrix A(i, i:n)
+ matrix_part_swap_rows(n, A, i, k, i, n - 1);
+
+ //Swap the matrix P at the rows i and pivot_index
+ matrix_swap_rows(n, P, i, k);
+
+ //Swap the rows of L in columns 0 through i-1. //FIXED
+ if (i >= 1) {
+
+ matrix_part_swap_rows(n, L, i, k, 0, i - 1);
+ }
+ changes++;
+ } //if
+
+ for (uint8_t j = i + 1; j < n; j++) {
+ //compute multipliers
+ multipliers_vec[j - 1] = A[j][i] / A[i][i];
+ //update the matrix A
+ for (uint8_t l = i + 1; l < n; l++) {
+ A[j][l] = A[j][l]
+ - multipliers_vec[j - 1]
+ * A[i][l];
+ }
+ //zero out the elements of column i, from row i+1 to n-1.
+ A[j][i] = 0;
+ // Set the matrix L with the multipliers
+ L[j][i] = multipliers_vec[j - 1];
+ }
+
+ } //for
+
+ return changes;
+}
diff --git a/linear_algebra/matrix_decompositions/qr_common.c b/linear_algebra/matrix_decompositions/qr_common.c
new file mode 100644
index 0000000000000000000000000000000000000000..27091fa79ea1480411389e9f2d58e81312d1e937
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/qr_common.c
@@ -0,0 +1,71 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Common definitions and implementations for the QR-decomposition.
+ * Provide necessary methods to construct Q- and R- matrices using.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <inttypes.h>
+#include <stdio.h>
+#include <string.h>
+
+#include "matrix.h"
+
+void qr_common_backward_subst(uint8_t m, uint8_t n, matrix_t U[][n],
+ matrix_t b[m], matrix_t x_sol[m])
+{
+ int8_t i;
+ uint8_t j;
+ matrix_t sum;
+
+ if (m != n) {
+ puts("The matrix should be square !!!");
+ return;
+ }
+
+ //clear x_sol
+ memset(x_sol, 0, sizeof(matrix_t) * m);
+
+ if (U[m - 1][m - 1] != 0) {
+ x_sol[m - 1] = b[m - 1] / U[m - 1][m - 1];
+ }
+
+ for (i = m - 2; i >= 0; i--) {
+ sum = 0.0;
+ for (j = i + 1; j < m; j++) {
+ sum += U[i][j] * x_sol[j];
+ }
+
+ if (U[i][i] != 0) {
+ x_sol[i] = (b[i] - sum) / U[i][i];
+ }
+ }
+}
+
+void qr_common_get_reduced_QR(uint8_t m, uint8_t n, matrix_t Q[m][m],
+ matrix_t R[m][n], matrix_t reduc_Q[m][n],
+ matrix_t reduc_R[n][n])
+{
+ if (m >= n) {
+ matrix_part_copy(m, m, Q, 0, m - 1, 0, n - 1, m, n, reduc_Q);
+ matrix_part_copy(m, n, R, 0, n - 1, 0, n - 1, n, n, reduc_R);
+ }
+ else {
+ puts(
+ "The row number of A should be greater than the column number of A");
+ }
+}
diff --git a/linear_algebra/matrix_decompositions/qr_givens.c b/linear_algebra/matrix_decompositions/qr_givens.c
new file mode 100644
index 0000000000000000000000000000000000000000..ff9d30a64a871a09c9734a97887da60a4cae15e8
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/qr_givens.c
@@ -0,0 +1,127 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Givens algorithm for the QR-decomposition.
+ * Provide necessary methods to construct Q- and R- matrices using
+ * Givens rotations. A = QR, where Q is an (m \f$\times \f$ n)-matrix with
+ * orthonormal columns and R is an (n \f$\times\f$ n) upper triangular matrix.
+ *
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+
+#include "matrix.h"
+#include "qr_givens.h"
+
+/*R will be stored in A */
+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)
+{
+ if (m < n) {
+ puts("A has more columns than rows");
+ return -1;
+ }
+
+ for (uint8_t j = 0; j < n; j++) {
+ for (uint8_t i = j + 1; i < m; i++) {
+ matrix_t c_s_t_r_vec[4];
+ qr_givens_get_params(A[j][j], A[i][j], c_s_t_r_vec);
+ A[i][j] = c_s_t_r_vec[2]; //A[i][j] = t
+ A[j][j] = c_s_t_r_vec[3]; //A[j][j] = r
+ if (j < n) {
+ for (uint8_t k = j + 1; k < n; k++) {
+ matrix_t c, s;
+ c = c_s_t_r_vec[0];
+ s = c_s_t_r_vec[1];
+ matrix_t tmp = A[j][k];
+ A[j][k] = c * A[j][k] + s * A[i][k];
+ A[i][k] = -s * tmp + c * A[i][k];
+ }
+ }
+
+ }
+ }
+
+ uint8_t max_m;
+ if (reduced) {
+ max_m = n;
+ }
+ else {
+ max_m = m;
+ }
+
+ if (!reduced) {
+ matrix_clear(m, max_m, Q);
+ }
+
+ matrix_part_copy(m, n, A, 0, m - 1, 0, n - 1, m, max_m, Q);
+
+ // Build R
+ matrix_get_upp_triang(m, n, A, A);
+
+ // Build Q
+ for (int16_t j = max_m - 1; j >= 0; j--) {
+ // Zero out column j from row 0 up to row j
+ for (uint8_t k = 0; k <= j; k++) {
+ Q[k][j] = 0;
+ }
+ Q[j][j] = 1;
+
+ for (int16_t i = m - 1; i >= j + 1; i--) {
+ matrix_t t = Q[i][j];
+ Q[i][j] = 0;
+ matrix_t c = 1 / sqrt(1 + t * t);
+ matrix_t s;
+ if (c != 0) {
+ s = c * t;
+ }
+ else {
+ s = 1;
+ }
+
+ for (uint8_t k = j; k < max_m; k++) {
+ matrix_t tmp = Q[j][k];
+ Q[j][k] = c * Q[j][k] - s * Q[i][k];
+ Q[i][k] = s * tmp + c * Q[i][k];
+ }
+ } //for
+ }
+
+ return 1;
+}
+
+void qr_givens_get_params(matrix_t xjj, matrix_t xij, matrix_t c_s_t_r_vec[])
+{
+ if (xjj != 0) {
+ matrix_t u;
+
+ c_s_t_r_vec[2] = xij / xjj; //t = xij/xjj
+ u = sqrt(1 + c_s_t_r_vec[2] * c_s_t_r_vec[2]); //u = sqrt(1 + t*t)
+ c_s_t_r_vec[0] = 1 / u; //c = 1/u
+ c_s_t_r_vec[1] = c_s_t_r_vec[0] * c_s_t_r_vec[2]; //s = c*t
+ c_s_t_r_vec[3] = u * xjj; //r = u*xjj
+ }
+ else {
+ c_s_t_r_vec[0] = 0; //c = 0
+ c_s_t_r_vec[1] = 1; //s = 1
+ c_s_t_r_vec[2] = FLT_MAX; //t = infinity
+ c_s_t_r_vec[3] = xij; //r = xij
+ }
+}
diff --git a/linear_algebra/matrix_decompositions/qr_householder.c b/linear_algebra/matrix_decompositions/qr_householder.c
new file mode 100644
index 0000000000000000000000000000000000000000..b0da7346dc42f3785762bb5cb25061bd60294e0d
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/qr_householder.c
@@ -0,0 +1,125 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Householder algorithm for the QR-decomposition.
+ * @details Provide necessary methods to construct Q- and R- matrices using
+ * Householder reflections. A = QR, where Q is an (m \f$\times\f$ n)-matrix with
+ * orthonormal columns and R is an (n \f$\times\f$ n) upper triangular matrix.
+ *
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <stdbool.h>
+#include <math.h>
+
+#include "matrix.h"
+#include "utils.h"
+
+/*R will be stored in A */
+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)
+{
+ if (m < n) {
+ puts("A has more columns than rows");
+ return -1;
+ }
+ matrix_t diag_R_vec[n];
+
+ for (uint8_t k = 0; k < n; k++) {
+ //compute the two-norm of the k-th column of the matrix A
+ double norm_2 = 0.0;
+ for (uint8_t i = k; i < m; i++) {
+ norm_2 = utils_get_save_square_root(norm_2, A[i][k]);
+
+ }
+
+ if (norm_2 != 0.0) {
+ // Form k-th Householder vector.
+ if (A[k][k] < 0) {
+ norm_2 = -norm_2;
+ }
+ for (int i = k; i < m; i++) {
+ A[i][k] /= norm_2;
+ }
+
+ A[k][k] += 1.0;
+
+ // Apply the Householder transformation to the rest columns of the matrix.
+ for (uint8_t j = k + 1; j < n; j++) {
+ double sigma = 0.0;
+ for (uint8_t i = k; i < m; i++) {
+ sigma += A[i][k] * A[i][j];
+ }
+ sigma = -sigma / A[k][k];
+ for (uint8_t i = k; i < m; i++) {
+ A[i][j] += sigma * A[i][k];
+ }
+ }
+
+ } // if
+
+ diag_R_vec[k] = -norm_2;
+
+ } //for
+
+ // Build Q
+ uint8_t max_n, max_m;
+ if (!reduced) {
+ max_n = m;
+ max_m = m;
+ }
+ else {
+ max_n = n;
+ max_m = n;
+ }
+
+ for (int k = max_n - 1; k >= 0; k--) {
+ for (int i = 0; i < m; i++) {
+ Q[i][k] = 0.0;
+ } // clear matrix Q
+ Q[k][k] = 1.0;
+ for (int j = k; j < max_n; j++) {
+ if (k < n) {
+ if (A[k][k] != 0) {
+ double s = 0.0;
+ for (int i = k; i < m; i++) {
+ s += A[i][k] * Q[i][j];
+ }
+ s = -s / A[k][k];
+ for (int i = k; i < m; i++) {
+ Q[i][j] += s * A[i][k];
+ }
+ } //if
+ } //if
+ }
+ }
+
+ // Build R
+ for (int i = 0; i < max_m; i++) {
+ for (int j = 0; j < n; j++) {
+ if (i == j) {
+ A[i][j] = diag_R_vec[i];
+ }
+ else if (i > j) {
+ A[i][j] = 0;
+ }
+ }
+ }
+
+ return 1;
+}
diff --git a/linear_algebra/matrix_decompositions/svd.c b/linear_algebra/matrix_decompositions/svd.c
new file mode 100644
index 0000000000000000000000000000000000000000..c9ed22f568e1d0f5408fbddf7ef5c89363fcd3ea
--- /dev/null
+++ b/linear_algebra/matrix_decompositions/svd.c
@@ -0,0 +1,809 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Algorithm for the Singular Value Decomposition (SVD).
+ * Provide necessary methods to compute the compact SVD of
+ * a matrix. A = U*S*V, where U is a (m x l) orthogonal matrix,
+ * S is a (l x l) diagonal matrix, V is a (l x n) orthogonal
+ * matrix, and l = min(m,n). The SVD is computed by using the
+ * Golub--Kahan--Reinsch algorithm that works in two phases:
+ * bidiagonalization and a reduction to the diagonal form phase.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdbool.h>
+#include <inttypes.h>
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+
+#include "svd.h"
+#include "matrix.h"
+#include "vector.h"
+
+static void svd_bidiagonal_trans_and_store_diag_elem(uint8_t m, uint8_t n,
+ matrix_t A[m][n],
+ uint8_t u_m, uint8_t u_n,
+ matrix_t U0[u_m][u_n],
+ matrix_t V0[n][n],
+ matrix_t sup_diag_elem_vec[],
+ matrix_t diag_elem_vec[]
+ );
+
+static void svd_setup_bi_diagonal_matrix(uint8_t m, uint8_t n, matrix_t A[m][n],
+ uint8_t u_n, matrix_t U1[m][u_n],
+ matrix_t V1[n][n],
+ matrix_t sup_diag_elem_vec[],
+ matrix_t diag_elem_vec[]
+ );
+
+static void svd_compute_singular_values_U_V(uint8_t m, uint8_t n, uint8_t u_n,
+ matrix_t U[m][u_n],
+ matrix_t V[n][n],
+ matrix_t sup_diag_elem_vec[],
+ matrix_t diag_elem_vec[]);
+
+static void svd_generate_S(uint8_t m, uint8_t n, matrix_t S[n][n],
+ uint8_t sing_vec_length, matrix_t singl_values_vec[]);
+
+static int32_t min(int32_t a, int32_t b)
+{
+ if (a < b) {
+ return a;
+ }
+ else {
+ return b;
+ }
+}
+
+static int32_t max(int32_t a, int32_t b)
+{
+ if (a > b) {
+ return a;
+ }
+ else {
+ return b;
+ }
+}
+
+static matrix_t get_abs_max(matrix_t arr[], uint8_t length)
+{
+ uint8_t i;
+ matrix_t max = fabs(arr[0]);
+
+ for (i = 1; i < length; i++) {
+ if (fabs(arr[i]) > max) {
+ max = fabs(arr[i]);
+ }
+ }
+
+ return max;
+}
+
+void svd_get_U_dim(uint8_t m, uint8_t n, matrix_dim_t *u_dim)
+{
+ u_dim->row_num = m;
+ u_dim->col_num = min(m, n);
+}
+
+void svd_get_S_dim(uint8_t m, uint8_t n, matrix_dim_t *s_dim)
+{
+ s_dim->row_num = min(m, n);
+ s_dim->col_num = n;
+}
+
+void svd_get_V_dim(uint8_t m, uint8_t n, matrix_dim_t *v_dim)
+{
+ v_dim->row_num = m;
+ v_dim->col_num = n;
+}
+
+uint8_t svd_get_single_values_num(uint8_t m, uint8_t n)
+{
+ return min(m + 1, n);
+}
+
+void svd(uint8_t m, uint8_t n, matrix_t A[m][n],
+ uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n],
+ matrix_t S[u_n][n], matrix_t V[n][n],
+ uint8_t sing_vec_length, matrix_t singl_values_vec[])
+{
+
+ int8_t u_n_min = min(m, n);
+ matrix_t sup_diag_elem_vec[n];
+
+ matrix_clear(m, u_n_min, U);
+ matrix_clear(n, n, S);
+ matrix_clear(n, n, V);
+ memset(singl_values_vec, 0, sizeof(matrix_t) * sing_vec_length);
+
+ svd_bidiagonal_trans_and_store_diag_elem(m, n, A, u_m, u_n, U, V,
+ sup_diag_elem_vec, singl_values_vec);
+
+ svd_setup_bi_diagonal_matrix(m, n, A, u_n, U, V, sup_diag_elem_vec,
+ singl_values_vec);
+
+ svd_compute_singular_values_U_V(m, n, u_n, U, V, sup_diag_elem_vec,
+ singl_values_vec);
+
+ svd_generate_S(u_n, n, S, sing_vec_length, singl_values_vec);
+}
+
+/**
+ * @brief Reduce the matrix A to a bi-diagonal form.
+ *
+ * @details The bi-diagonal transformation is placed in the U0 and V0 matrices
+ * for a subsequent back multiplication. The super-diagonal and the
+ * diagonal elements are saved in the sup_diag_elem_vec and
+ * diag_elem_vec vectors.
+ *
+ * @param[in] m row number of the matrix A.
+ * @param[in] n column number of the matrix A.
+ * @param[in, out] A[][n] pointer to the matrix A.
+ * @param[in] u_m row number of the matrix U0.
+ * @param[in] u_n column number of the matrix U0.
+ * @param[out] U0[][u_n] pointer to the matrix U0.
+ * @param[out] V0[][u_n] pointer to the matrix V0.
+ * @param[out] sup_diag_elem_vec[] pointer to the vector saving the
+ * super-diagonal elements.
+ * @param[out] diag_elem_vec[] pointer to the vector saving the
+ * diagonal elements.
+ *
+ */
+static void svd_bidiagonal_trans_and_store_diag_elem(uint8_t m, uint8_t n,
+ matrix_t A[m][n],
+ uint8_t u_m, uint8_t u_n,
+ matrix_t U0[u_m][u_n],
+ matrix_t V0[n][n],
+ matrix_t sup_diag_elem_vec[],
+ matrix_t diag_elem_vec[])
+{
+ int col_trans_num = min(m - 1, n);
+ int row_trans_num = max(0, min(n - 2, m));
+ int upper_lim = max(col_trans_num, row_trans_num);
+
+ for (int index = 0; index < upper_lim; index++) {
+ if (index < col_trans_num) {
+ diag_elem_vec[index] = 0;
+ for (int row_num = index; row_num < m; row_num++) {
+ diag_elem_vec[index] = hypot(
+ diag_elem_vec[index],
+ A[row_num][index]);
+ }
+ if (diag_elem_vec[index] != 0.0) {
+ if (A[index][index] < 0.0) {
+ diag_elem_vec[index] =
+ -diag_elem_vec[index];
+ }
+ for (int row_num = index; row_num < m;
+ row_num++) {
+ A[row_num][index] /=
+ diag_elem_vec[index];
+ }
+ A[index][index] += 1.0;
+ }
+ diag_elem_vec[index] = -diag_elem_vec[index];
+ }
+ for (int col_num = index + 1; col_num < n; col_num++) {
+ if ((index < col_trans_num)
+ & (diag_elem_vec[index] != 0.0)) {
+
+ /* Transformation */
+ matrix_t transf_val = 0;
+ for (int i = index; i < m; i++) {
+ transf_val += A[i][index]
+ * A[i][col_num];
+ }
+ transf_val = -transf_val / A[index][index];
+ for (int i = index; i < m; i++) {
+ A[i][col_num] += transf_val
+ * A[i][index];
+ }
+ }
+ sup_diag_elem_vec[col_num] = A[index][col_num];
+ }
+ if (index < col_trans_num) {
+
+ /* Place a part of the bi-diagonal transformation in U0 */
+ for (int row_num = index; row_num < m; row_num++) {
+ U0[row_num][index] = A[row_num][index];
+ }
+ }
+ if (index < row_trans_num) {
+ sup_diag_elem_vec[index] = 0;
+ for (int i = index + 1; i < n; i++) {
+ sup_diag_elem_vec[index] = hypot(
+ sup_diag_elem_vec[index],
+ sup_diag_elem_vec[i]);
+ }
+ if (sup_diag_elem_vec[index] != 0.0) {
+ if (sup_diag_elem_vec[index + 1] < 0.0) {
+ sup_diag_elem_vec[index] =
+ -sup_diag_elem_vec[index];
+ }
+ for (int i = index + 1; i < n; i++) {
+ sup_diag_elem_vec[i] /=
+ sup_diag_elem_vec[index];
+ }
+ sup_diag_elem_vec[index + 1] += 1.0;
+ }
+ sup_diag_elem_vec[index] = -sup_diag_elem_vec[index];
+ if ((index + 1 < m)
+ & (sup_diag_elem_vec[index] != 0.0)) {
+ matrix_t transf_vec[m];
+
+ /* Compute the transformation vector */
+ for (int i = index + 1; i < m; i++) {
+ transf_vec[i] = 0.0;
+ }
+ for (int j = index + 1; j < n; j++) {
+ for (int i = index + 1; i < m; i++) {
+ transf_vec[i] +=
+ sup_diag_elem_vec[j]
+ * A[i][j];
+ }
+ }
+
+ /* Reduce the matrix */
+ for (int col_num = index + 1; col_num < n;
+ col_num++) {
+ matrix_t t =
+ -sup_diag_elem_vec[col_num]
+ / sup_diag_elem_vec[index
+ + 1];
+ for (int row_num = index + 1;
+ row_num < m;
+ row_num++) {
+ A[row_num][col_num] +=
+ t
+ * transf_vec[row_num];
+ }
+ }
+ }
+
+ /* Place a part of the bi-diagonal transformation in V0 */
+ for (int row_num = index + 1; row_num < n; row_num++) {
+ V0[row_num][index] = sup_diag_elem_vec[row_num];
+ }
+ }
+ }
+}
+
+/**
+ * @brief Setup the final bi-diagonal transformation of the matrix.
+ *
+ * @details The final transformation is placed in the U1 and V1 matrices.
+ * The final transformation is based on the reduced bi-diagonal
+ * matrix computed in a previous step.
+ *
+ * @param[in] m row number of the matrix A.
+ * @param[in] n column number of the matrix A.
+ * @param[in] A[][n] pointer to the matrix A.
+ * @param[in] u_n column number of the matrix U1.
+ * @param[out] U1[][u_n] pointer to the matrix U1.
+ * @param[out] V1[][n] pointer to the matrix V1.
+ * @param[in,out] sup_diag_elem_vec[] pointer to the vector saving the
+ * super-diagonal elements.
+ * @param[in,out] diag_elem_vec[] pointer to the vector saving the
+ * diagonal elements.
+ *
+ */
+static void svd_setup_bi_diagonal_matrix(uint8_t m, uint8_t n, matrix_t A[m][n],
+ uint8_t u_n, matrix_t U1[m][u_n],
+ matrix_t V1[n][n],
+ matrix_t sup_diag_elem_vec[],
+ matrix_t diag_elem_vec[])
+{
+
+ int col_trans_num = min(m - 1, n);
+ int row_trans_num = max(0, min(n - 2, m));
+ int8_t upp_col_num = min(m, n);
+ int p_index = min(n, m + 1);
+
+ if (col_trans_num < n) {
+ diag_elem_vec[col_trans_num] = A[col_trans_num][col_trans_num];
+ }
+ if (m < p_index) {
+ diag_elem_vec[p_index - 1] = 0.0;
+ }
+ if (row_trans_num + 1 < p_index) {
+ sup_diag_elem_vec[row_trans_num] =
+ A[row_trans_num][p_index - 1];
+ }
+ sup_diag_elem_vec[p_index - 1] = 0.0;
+
+ /* Generate the matrix U1 */
+ for (int j = col_trans_num; j < upp_col_num; j++) {
+ for (int i = 0; i < m; i++) {
+ U1[i][j] = 0.0;
+ }
+ U1[j][j] = 1.0;
+ }
+ for (int k = col_trans_num - 1; k >= 0; k--) {
+ if (diag_elem_vec[k] != 0.0) {
+ for (int j = k + 1; j < upp_col_num; j++) {
+ matrix_t factor_t = 0;
+ for (int i = k; i < m; i++) {
+ factor_t += U1[i][k] * U1[i][j];
+ }
+ factor_t = -factor_t / U1[k][k];
+ for (int i = k; i < m; i++) {
+ U1[i][j] += factor_t * U1[i][k];
+ }
+ }
+ for (int i = k; i < m; i++) {
+ U1[i][k] = -U1[i][k];
+ }
+ U1[k][k] += 1.0;
+ for (int i = 0; i < k - 1; i++) {
+ U1[i][k] = 0.0;
+ }
+ }
+ else {
+ for (int i = 0; i < m; i++) {
+ U1[i][k] = 0.0;
+ }
+ U1[k][k] = 1.0;
+ }
+ }
+
+ /* Generate the matrix V1 */
+ for (int k = n - 1; k >= 0; k--) {
+ if ((k < row_trans_num) & (sup_diag_elem_vec[k] != 0.0)) {
+ for (int j = k + 1; j < n; j++) {
+ matrix_t t = 0;
+ for (int i = k + 1; i < n; i++) {
+ t += V1[i][k] * V1[i][j];
+ }
+ t = -t / V1[k + 1][k];
+ for (int i = k + 1; i < n; i++) {
+ V1[i][j] += t * V1[i][k];
+ }
+ }
+ }
+ for (int i = 0; i < n; i++) {
+ V1[i][k] = 0.0;
+ }
+ V1[k][k] = 1.0;
+ }
+}
+
+/**
+ * @brief Compute the singular values and the matrices U and V.
+ *
+ * @details The computed bi-diagonal matrix in the first phase is reduced to
+ * a diagonal form. In the bi-diagonal procedure, we distinguish
+ * between the splitting, the cancellation and the negligibility
+ * steps. The singular values will be stored in the vector of
+ * diagonal elements.
+ *
+ *
+ * @param[in] m row number of the matrix U.
+ * @param[in] n row and column number of the matrix
+ * V.
+ * @param[in] u_n column number of the matrix U.
+ * @param[in,out] U[][u_n] pointer to the matrix U.
+ * @param[in,out] V[][n] pointer to the matrix V.
+ * @param[in,out] sup_diag_elem_vec[] pointer to the vector saving the
+ * super-diagonal elements.
+ * @param[in,out] diag_elem_vec[] pointer to the vector saving the
+ * diagonal and singular elements.
+ *
+ */
+static void svd_compute_singular_values_U_V(uint8_t m, uint8_t n, uint8_t u_n,
+ matrix_t U[m][u_n],
+ matrix_t V[n][n],
+ matrix_t sup_diag_elem_vec[],
+ matrix_t diag_elem_vec[])
+{
+ int p_index = min(n, m + 1);
+ int max_k = p_index - 1;
+ int iter_num = 0;
+ matrix_t eps = pow(2.0, -52.0);
+ matrix_t tol = pow(2.0, -966.0);
+
+ while (p_index > 0) {
+ int k_index, case_num;
+
+ for (k_index = p_index - 2; k_index >= -1; k_index--) {
+ if (k_index == -1) {
+ break;
+ }
+ if (fabs(sup_diag_elem_vec[k_index])
+ <=
+ tol
+ + eps
+ * (fabs(
+ diag_elem_vec[k_index])
+ + fabs(
+ diag_elem_vec[k_index
+ + 1]))) {
+ sup_diag_elem_vec[k_index] = 0.0;
+ break;
+ }
+ }
+ if (k_index == p_index - 2) {
+ case_num = SVD_ORDER_ABSOLUTE_SING_VALUES; //4
+ }
+ else {
+ int ks_index;
+ for (ks_index = p_index - 1; ks_index >= k_index;
+ ks_index--) {
+ if (ks_index == k_index) {
+ break;
+ }
+ matrix_t t =
+ (ks_index != p_index ?
+ fabs(
+ sup_diag_elem_vec[ks_index]) :
+ 0.)
+ +
+ (ks_index
+ != k_index
+ + 1 ?
+ fabs(
+ sup_diag_elem_vec[ks_index
+ - 1]) :
+ 0.);
+ if (fabs(diag_elem_vec[ks_index])
+ <= tol + eps * t) {
+ diag_elem_vec[ks_index] = 0.0;
+ break;
+ }
+ }
+ if (ks_index == k_index) {
+ case_num = SVD_QR_STEP; //3
+ }
+ else if (ks_index == p_index - 1) {
+ case_num = SVD_COMPUTE_NEGLIGIBLE_VALUES; //1
+ }
+ else {
+ case_num = SVD_SPLIT_AT_NEGLIGIBLE_VALUES; //2
+ k_index = ks_index;
+ }
+ }
+ k_index++;
+
+ switch (case_num) {
+
+ /* Compute negligible values */
+ case SVD_COMPUTE_NEGLIGIBLE_VALUES: {
+ matrix_t f_val = sup_diag_elem_vec[p_index - 2];
+ sup_diag_elem_vec[p_index - 2] = 0.0;
+ for (int j = p_index - 2; j >= k_index; j--) {
+ matrix_t t_val = hypot(diag_elem_vec[j], f_val);
+ matrix_t split_cs_fac = diag_elem_vec[j]
+ / t_val;
+ matrix_t split_sn_fac = f_val / t_val;
+ diag_elem_vec[j] = t_val;
+ if (j != k_index) {
+ f_val = -split_sn_fac
+ * sup_diag_elem_vec[j
+ - 1];
+ sup_diag_elem_vec[j - 1] = split_cs_fac
+ * sup_diag_elem_vec[j
+ - 1];
+ }
+ /* Compute the V matrix */
+ for (int i = 0; i < n; i++) {
+ t_val =
+ split_cs_fac * V[i][j]
+ + split_sn_fac
+ * V[i][p_index
+ - 1];
+ V[i][p_index - 1] =
+ -split_sn_fac * V[i][j]
+ + split_cs_fac
+ * V[i][p_index
+ - 1];
+ V[i][j] = t_val;
+ }
+ }
+ }
+ break;
+
+ /* Split at negligible values sup_diag_elem_vec(k) */
+ case SVD_SPLIT_AT_NEGLIGIBLE_VALUES: {
+ matrix_t f_val = sup_diag_elem_vec[k_index - 1];
+ sup_diag_elem_vec[k_index - 1] = 0.0;
+ for (int j = k_index; j < p_index; j++) {
+ matrix_t t_val = hypot(diag_elem_vec[j], f_val);
+ matrix_t split_cs_fac = diag_elem_vec[j]
+ / t_val;
+ matrix_t split_sn_fac = f_val / t_val;
+ diag_elem_vec[j] = t_val;
+ f_val = -split_sn_fac * sup_diag_elem_vec[j];
+ sup_diag_elem_vec[j] = split_cs_fac
+ * sup_diag_elem_vec[j];
+
+ /* Compute the U matrix */
+ for (int i = 0; i < m; i++) {
+ t_val =
+ split_cs_fac * U[i][j]
+ + split_sn_fac
+ * U[i][k_index
+ - 1];
+ U[i][k_index - 1] =
+ -split_sn_fac * U[i][j]
+ + split_cs_fac
+ * U[i][k_index
+ - 1];
+ U[i][j] = t_val;
+ }
+ }
+ }
+ break;
+
+ /* The QR Step */
+ case SVD_QR_STEP: {
+
+ matrix_t buff[5] = { diag_elem_vec[p_index - 1],
+ diag_elem_vec[p_index - 2],
+ sup_diag_elem_vec[p_index - 2],
+ diag_elem_vec[k_index],
+ sup_diag_elem_vec[k_index] };
+ matrix_t scale = get_abs_max(buff, 5);
+
+ matrix_t p_diag_elem = diag_elem_vec[p_index - 1]
+ / scale;
+ matrix_t p_minus_diag_elem = diag_elem_vec[p_index - 2]
+ / scale;
+ matrix_t p_minus_sup_diag_elem =
+ sup_diag_elem_vec[p_index - 2] / scale;
+ matrix_t k_diag_elem = diag_elem_vec[k_index] / scale;
+ matrix_t k_sup_diag_elem = sup_diag_elem_vec[k_index]
+ / scale;
+ matrix_t d_val = ((p_minus_diag_elem + p_diag_elem)
+ * (p_minus_diag_elem - p_diag_elem)
+ + p_minus_sup_diag_elem
+ * p_minus_sup_diag_elem)
+ / 2.0;
+ matrix_t e_val = (p_diag_elem * p_minus_sup_diag_elem)
+ * (p_diag_elem * p_minus_sup_diag_elem);
+ matrix_t wilkinson_shift = 0.0;
+ if ((d_val != 0.0) | (e_val != 0.0)) {
+ wilkinson_shift = sqrt(d_val * d_val + e_val);
+ if (d_val < 0.0) {
+ wilkinson_shift = -wilkinson_shift;
+ }
+ wilkinson_shift = e_val
+ / (d_val + wilkinson_shift);
+ }
+ matrix_t f_val = (k_diag_elem + p_diag_elem)
+ * (k_diag_elem - p_diag_elem)
+ + wilkinson_shift;
+ matrix_t g_val = k_diag_elem * k_sup_diag_elem;
+
+ for (int j = k_index; j < p_index - 1; j++) {
+ matrix_t t = hypot(f_val, g_val);
+ matrix_t split_cs_fac = f_val / t;
+ matrix_t split_sn_fac = g_val / t;
+ if (j != k_index) {
+ sup_diag_elem_vec[j - 1] = t;
+ }
+ f_val =
+ split_cs_fac * diag_elem_vec[j]
+ + split_sn_fac
+ * sup_diag_elem_vec[j];
+ sup_diag_elem_vec[j] =
+ split_cs_fac
+ * sup_diag_elem_vec[j]
+ - split_sn_fac
+ * diag_elem_vec[j];
+ g_val = split_sn_fac * diag_elem_vec[j + 1];
+ diag_elem_vec[j + 1] = split_cs_fac
+ * diag_elem_vec[j + 1];
+
+ /* Compute the V matrix */
+ for (int i = 0; i < n; i++) {
+ t =
+ split_cs_fac * V[i][j]
+ + split_sn_fac
+ * V[i][j
+ + 1];
+ V[i][j + 1] =
+ -split_sn_fac * V[i][j]
+ + split_cs_fac
+ * V[i][j
+ + 1];
+ V[i][j] = t;
+ }
+
+ t = hypot(f_val, g_val);
+ split_cs_fac = f_val / t;
+ split_sn_fac = g_val / t;
+ diag_elem_vec[j] = t;
+ f_val =
+ split_cs_fac
+ * sup_diag_elem_vec[j]
+ + split_sn_fac
+ * diag_elem_vec[j
+ + 1];
+ diag_elem_vec[j + 1] =
+ -split_sn_fac
+ * sup_diag_elem_vec[j]
+ + split_cs_fac
+ * diag_elem_vec[j
+ + 1];
+ g_val = split_sn_fac * sup_diag_elem_vec[j + 1];
+ sup_diag_elem_vec[j + 1] = split_cs_fac
+ * sup_diag_elem_vec[j + 1];
+ if (j < m - 1) {
+ for (int i = 0; i < m; i++) {
+ t =
+ split_cs_fac
+ * U[i][j]
+ + split_sn_fac
+ * U[i][j
+ + 1];
+ U[i][j + 1] =
+ -split_sn_fac
+ * U[i][j]
+ + split_cs_fac
+ * U[i][j
+ + 1];
+ U[i][j] = t;
+ }
+ }
+ }
+ sup_diag_elem_vec[p_index - 2] = f_val;
+ iter_num++;
+ }
+ break;
+
+ /* Order the absolute singular values */
+ case SVD_ORDER_ABSOLUTE_SING_VALUES: {
+
+ if (diag_elem_vec[k_index] <= 0.0) {
+ diag_elem_vec[k_index] =
+ (diag_elem_vec[k_index] < 0.0 ?
+ -diag_elem_vec[k_index] :
+ 0.0);
+
+ for (int i = 0; i < n; i++) {
+ V[i][k_index] = -V[i][k_index];
+ }
+ }
+
+ while (k_index < max_k) {
+ if (diag_elem_vec[k_index]
+ >= diag_elem_vec[k_index + 1]) {
+ break;
+ }
+ matrix_t k_diag_elem = diag_elem_vec[k_index];
+ diag_elem_vec[k_index] = diag_elem_vec[k_index
+ + 1];
+ diag_elem_vec[k_index + 1] = k_diag_elem;
+ if (k_index < n - 1) {
+ for (int i = 0; i < n; i++) {
+ k_diag_elem = V[i][k_index + 1];
+ V[i][k_index + 1] =
+ V[i][k_index];
+ V[i][k_index] = k_diag_elem;
+ }
+ }
+ if (k_index < m - 1) {
+ for (int i = 0; i < m; i++) {
+ k_diag_elem = U[i][k_index + 1];
+ U[i][k_index + 1] =
+ U[i][k_index];
+ U[i][k_index] = k_diag_elem;
+ }
+ }
+ k_index++;
+ }
+ iter_num = 0;
+ p_index--;
+ }
+ break;
+ }
+ }
+}
+
+/**
+ * @brief Compute the rectangular diagonal matrix S.
+ *
+ * @details The calculation of the matrix S is based on the three previous
+ * steps: the bi-diagonal reduction, the final bi-diagonal
+ * transformation, and the computation of the singular values.
+ *
+ * @param[in] m row number of the matrix to transform in SVD.
+ * @param[in] n row and column number of the matrix S.
+ * @param[in] S[][n] pointer to the matrix S.
+ * @param[in] sing_vec_length length of the vector saving the singular
+ * values.
+ * @param[in] singl_values_vec[] pointer to the vector saving the singular
+ * elements.
+ *
+ */
+static void svd_generate_S(uint8_t m, uint8_t n, matrix_t S[n][n],
+ uint8_t sing_vec_length, matrix_t singl_values_vec[])
+{
+ uint8_t i;
+
+ for (i = 0; i < m; i++) {
+ for (uint8_t j = 0; j < n; j++) {
+ S[i][j] = 0.0;
+ }
+ if (i < sing_vec_length) {
+ S[i][i] = singl_values_vec[i];
+ }
+ }
+}
+
+void svd_get_reciproc_singular_values(uint8_t m, uint8_t n, uint8_t length,
+ matrix_t singl_values_arr[],
+ matrix_t recip_singl_values_arr[]
+ )
+{
+ uint8_t i;
+ matrix_t eps = pow(2.0, -52.0);
+ matrix_t tol = max(m, n) * singl_values_arr[0] * eps;
+
+ for (i = 0; i < length; i++) {
+ if (fabs(singl_values_arr[i]) >= tol) {
+ recip_singl_values_arr[i] = 1.0 / singl_values_arr[i];
+ }
+ else {
+ recip_singl_values_arr[i] = 0;
+ }
+ }
+}
+
+void svd_compute_print_U_S_V_s(uint8_t m, uint8_t n, matrix_t matrix_arr[m][n],
+ uint8_t i)
+{
+ printf(
+ "########################## Test %d ##########################\n",
+ i);
+ matrix_dim_t u_dim;
+ svd_get_U_dim(m, n, &u_dim);
+ matrix_t U[u_dim.row_num][u_dim.col_num];
+ matrix_t S[u_dim.col_num][n];
+ matrix_t V[n][n];
+ uint8_t s_length = svd_get_single_values_num(m, n);
+
+ printf("matrix%d =\n", i);
+ matrix_print(m, n, matrix_arr);
+ puts("");
+ matrix_t s[s_length];
+ svd(m, n, matrix_arr, u_dim.row_num, u_dim.col_num, U, S, V, s_length,
+ s);
+
+ printf("U%d =\n", i);
+ matrix_print(u_dim.row_num, u_dim.col_num, U);
+ puts("");
+
+ printf("S%d =\n", i);
+ matrix_print(u_dim.col_num, n, S);
+ puts("");
+
+ printf("V%d =\n", i);
+ matrix_print(n, n, V);
+ puts("");
+
+ printf("s%d = ", i);
+ printf("{");
+ for (uint8_t i = 0; i < s_length; i++) {
+ printf("%5.4f ", s[i]);
+ if (i < s_length - 1) {
+ printf(", ");
+ }
+ }
+ printf("}");
+
+ puts("");
+}
diff --git a/linear_algebra/pseudo_inverse/doc.txt b/linear_algebra/pseudo_inverse/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c95aa897927459dec98fd27f9905f9d2e7351c7b
--- /dev/null
+++ b/linear_algebra/pseudo_inverse/doc.txt
@@ -0,0 +1,20 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup pseudo_inverse PSEUDO_INVERSE
+ * @ingroup linear_algebra
+ *
+ * @brief Algorithms to calculate the pseudo-inverse of a matrix.
+ *
+ * @details The pseudo-inverse matrix can be computed using QR-decomposition or SVD algorithms.
+ *
+ * @author Zakaria Kasmi
+ *
+ */
diff --git a/linear_algebra/pseudo_inverse/include/moore_penrose_pseudo_inverse.h b/linear_algebra/pseudo_inverse/include/moore_penrose_pseudo_inverse.h
new file mode 100644
index 0000000000000000000000000000000000000000..a53da965169b9ec30173224f26ef89007a620aa4
--- /dev/null
+++ b/linear_algebra/pseudo_inverse/include/moore_penrose_pseudo_inverse.h
@@ -0,0 +1,93 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup pseudo_inverse
+ * @{
+ *
+ * @file
+ * @brief Moore--Penrose algorithm to compute the pseudo-inverse of a matrix.
+ *
+ * @details The computation of the pseudo-inverse is based on the
+ * Singular Value Decomposition (SVD).
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef MOORE_PENROSE_PSEUDO_INVERSE_H_
+#define MOORE_PENROSE_PSEUDO_INVERSE_H_
+
+/**
+ * The maximal row number allowed.
+ */
+#define MAX_ROW_NUM 23
+
+/**
+ * The maximal column number allowed.
+ */
+#define MAX_COL_NUM 23
+
+/* define error numbers */
+/**
+ * The Moore--Penrose inverse is successfully completed.
+ */
+#define MOORE_PENROSE_PSEUDO_COMP_SUCCESS 1
+
+/**
+ * The maximal row number allowed is exceeded.
+ */
+#define MOORE_PENROSE_PSEUDO_MAX_ALLOW_ROW_COL_EXCEEED -1
+
+/**
+ * The transposed matrix should be delivered.
+ */
+#define MOORE_PENROSE_PSEUDO_GIVE_MATRIX_TRANSPOSE -2
+
+/**
+ * Invalid rank value of a matrix.
+ */
+#define MOORE_PENROSE_INVALID_RANK_VALUE -3
+
+#include "matrix.h"
+
+/**
+ * @brief Calculate the Moore--Penrose inverse of a rectangular matrix.
+ * The computation of the Moore--Penrose inverse is based on the
+ * Golub--Kahan--Reinsch SVD algorithm.
+ *
+ * @param[in] m row number of the matrix to inverse.
+ * @param[in] n column number of the matrix to inverse.
+ * @param[in] A[][] pointer to the matrix A.
+ * @param[out] pinv_A[][] pointer to the pseudo-inverse matrix.
+ *
+ * @return 1, if the computation of the Moore-Penrose inverse is successful.
+ * @return -1, if the maximal, allowed column or row number is exceeded.
+ * @return -2, if the matrix is underdetermined (m<n).
+ * @return -3, if the rank of the matrix is equal to 0.
+ *
+ */
+int8_t moore_penrose_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
+ matrix_t pinv_A[n][m]);
+
+/**
+ * @brief Compute and print the Moore--Penrose inverse of a matrix.
+ *
+ *
+ * @param[in] m row number of the matrix.
+ * @param[in] n column number of the matrix.
+ * @param[in] matrix[][] pointer to the matrix.
+ * @param[in] i label.
+ *
+ */
+void moore_penrose_pinv_compute_print(uint8_t m, uint8_t n,
+ matrix_t matrix[m][n], uint8_t i);
+
+#endif /* MOORE_PENROSE_PSEUDO_INVERSE_H_ */
diff --git a/linear_algebra/pseudo_inverse/include/pseudo_inverse.h b/linear_algebra/pseudo_inverse/include/pseudo_inverse.h
new file mode 100644
index 0000000000000000000000000000000000000000..fa9f3d7fc5b2ce73b2ecbcba794a64c3f78520d8
--- /dev/null
+++ b/linear_algebra/pseudo_inverse/include/pseudo_inverse.h
@@ -0,0 +1,37 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup pseudo_inverse
+ * @{
+ *
+ * @file
+ * @brief Compute the pseudo-inverse of a matrix.
+ *
+ * @details The pseudo-inverse matrix can be computed using the Singular Value Decomposition
+ * (SVD), Householder, or Givens algorithms.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+#ifndef PSEUDO_INVERSE_H_
+#define PSEUDO_INVERSE_H_
+
+/**
+ * Possible algorithms to compute the pseudo-inverse matrix.
+ */
+enum ALGORITHM {
+ Moore_Penrose, /**< Moore--Penrose algorithm */
+ Householder, /**< Householder algorithm */
+ Givens, /**< Givens algorithm */
+ Gauss /**< Gaussian elimination with pivoting algorithm */
+};
+
+#endif /* PSEUDO_INVERSE_H_ */
diff --git a/linear_algebra/pseudo_inverse/include/qr_pseudo_inverse.h b/linear_algebra/pseudo_inverse/include/qr_pseudo_inverse.h
new file mode 100644
index 0000000000000000000000000000000000000000..b04af0cc04c826752621c29b17f072e162a4ff39
--- /dev/null
+++ b/linear_algebra/pseudo_inverse/include/qr_pseudo_inverse.h
@@ -0,0 +1,51 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup pseudo_inverse
+ * @{
+ *
+ * @file
+ * @brief QR decomposition algorithms to compute the pseudo-inverse of a matrix.
+ *
+ * @details The computation of the pseudo-inverse is implemented using the Householder or
+ * the Givens algorithms.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef QR_GET_PINV_H_
+#define QR_GET_PINV_H_
+
+#include <inttypes.h>
+
+#include "matrix.h"
+#include "qr_common.h"
+
+/**
+ * @brief Calculate the pseudo inverse of a rectangular matrix using the QR decomposition.
+ * @details The computation of the pseudo inverse is based on the Householder or Givens
+ * algorithm.
+ *
+ * @param[in] m row number of the matrix to inverse.
+ * @param[in] n column number of the matrix to inverse.
+ * @param[in] A[][] pointer to the matrix A.
+ * @param[out] pinv_A[][] pointer to the pseudo-inverse matrix.
+ * @param[in] algo choice between the Householder or Givens algorithms.
+ *
+ * @return 1, if computing the pseudo-inverse matrix is successful.
+ * @return -1, if computing the pseudo-inverse matrix is not successful.
+ *
+ */
+int8_t qr_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
+ matrix_t pinv_A[n][m], enum QR_ALGORITHM algo);
+
+#endif /* QR_GET_PINV_H_ */
diff --git a/linear_algebra/pseudo_inverse/moore_penrose_pseudo_inverse.c b/linear_algebra/pseudo_inverse/moore_penrose_pseudo_inverse.c
new file mode 100644
index 0000000000000000000000000000000000000000..13ea72050932c3902267f86819bb4f7e830e8897
--- /dev/null
+++ b/linear_algebra/pseudo_inverse/moore_penrose_pseudo_inverse.c
@@ -0,0 +1,176 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Moore--Penrose algorithm to compute the pseudo-inverse of a
+ * rectangular matrix.
+ * @details The computation of the pseudo-inverse is based on the
+ * Singular Value Decomposition (SVD).
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdbool.h>
+#include <stdio.h>
+
+#include "moore_penrose_pseudo_inverse.h"
+#include "matrix.h"
+#include "svd.h"
+
+static int32_t min(int32_t a, int32_t b);
+static int8_t moore_penrose_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
+ uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n],
+ matrix_t S[u_n][n], matrix_t V[n][n],
+ uint8_t s_length, matrix_t s[],
+ matrix_t pinv_A[n][m]);
+
+int8_t moore_penrose_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
+ matrix_t pinv_A[n][m])
+{
+ matrix_dim_t udim;
+ int8_t answer = 0;
+
+ if (m < n) { /*
+ * first compute the pseudo inverse of the transposed of the
+ * matrix A, than transpose the computed pseudo inverse.
+ */
+ uint8_t trans_m = n;
+ uint8_t trans_n = m;
+ matrix_t trans_A[n][m];
+ matrix_t trans_A_pinv[trans_n][trans_m];
+
+ matrix_transpose(m, n, A, trans_A);
+ svd_get_U_dim(trans_m, trans_n, &udim);
+ matrix_t U[udim.row_num][udim.col_num];
+
+ matrix_t S[udim.col_num][trans_n];
+ matrix_t V[trans_n][trans_n];
+ uint8_t s_length = svd_get_single_values_num(trans_m, trans_n);
+ matrix_t s[s_length];
+
+ answer = moore_penrose_pinv(trans_m, trans_n, trans_A,
+ udim.row_num, udim.col_num, U, S, V,
+ s_length, s, trans_A_pinv);
+ matrix_transpose(trans_n, trans_m, trans_A_pinv, pinv_A);
+
+ return answer;
+ }
+ else { /* compute the pseudo inverse */
+ svd_get_U_dim(m, n, &udim);
+ matrix_t U[udim.row_num][udim.col_num];
+ matrix_t S[udim.col_num][n];
+ matrix_t V[n][n];
+ uint8_t s_length = svd_get_single_values_num(m, n);
+ matrix_t s[s_length];
+ answer = moore_penrose_pinv(m, n, A,
+ udim.row_num, udim.col_num, U, S, V,
+ s_length, s, pinv_A);
+
+ return answer;
+ }
+
+}
+
+/**
+ * @brief Calculate the Moore--Penrose inverse of a (m x n) matrix with (m>n).
+ * @details To compute the pseudo-inverse of a matrix with (m<n), the
+ * transposed of the matrix must be delivered.
+ *
+ * @param[in] m row number of the matrix to inverse.
+ * @param[in] n column number of the matrix to inverse.
+ * @param[in] A[][n] pointer to the matrix A.
+ * @param[out] pinv_A[][m] pointer to the pseudo-inverse matrix.
+ *
+ * @return 1, if the computation of the Moore-Penrose inverse is successful.
+ * @return -1, if the maximal, allowed column or row number is exceeded.
+ * @return -2, if the matrix is underdetermined (m<n).
+ * @return -3, if the rank of the matrix is equal to 0.
+ *
+ */
+static int8_t moore_penrose_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
+ uint8_t u_m, uint8_t u_n, matrix_t U[u_m][u_n],
+ matrix_t S[u_n][n], matrix_t V[n][n],
+ uint8_t s_length, matrix_t s[],
+ matrix_t pinv_A[n][m])
+{
+
+ matrix_dim_t u_dim;
+ matrix_t rec_s[s_length];
+
+ uint8_t i, j, k, col_min;
+
+ if ((m > MAX_ROW_NUM) | (n > MAX_COL_NUM)) {
+
+ printf(
+ "The maximal, allowed number of the rows and columns is %d\n",
+ MAX_ROW_NUM);
+ return MOORE_PENROSE_PSEUDO_MAX_ALLOW_ROW_COL_EXCEEED;
+ }
+
+ if (m < n) {
+
+ puts("Please, give the transposed matrix");
+ return MOORE_PENROSE_PSEUDO_GIVE_MATRIX_TRANSPOSE;
+ }
+ else {
+ svd_get_U_dim(m, n, &u_dim);
+ s_length = svd_get_single_values_num(m, n);
+
+ svd(m, n, A, u_dim.row_num, u_dim.col_num, U, S, V, s_length,
+ s);
+ }
+
+ if (matrix_get_rank(m, n, s, s_length) < 1) {
+ puts("The minimum rank-value of a matrix is one");
+ return MOORE_PENROSE_INVALID_RANK_VALUE;
+ }
+
+ svd_get_reciproc_singular_values(m, n, s_length, s, rec_s);
+
+ col_min = min(n, u_dim.col_num);
+
+ matrix_clear(n, m, pinv_A);
+
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < u_dim.row_num; j++) {
+ for (k = 0; k < col_min; k++) {
+ pinv_A[i][j] += V[i][k] * rec_s[k] * U[j][k];
+
+ }
+ }
+ }
+
+ return MOORE_PENROSE_PSEUDO_COMP_SUCCESS;
+}
+
+static int32_t min(int32_t a, int32_t b)
+{
+ if (a < b) {
+ return a;
+ }
+ else {
+ return b;
+ }
+}
+
+void moore_penrose_pinv_compute_print(uint8_t m, uint8_t n,
+ matrix_t matrix[m][n], uint8_t i)
+{
+ matrix_t pinv_mat[n][m];
+
+ moore_penrose_get_pinv(m, n, matrix, pinv_mat);
+ printf("pinv%d = ", i);
+ matrix_flex_print(n, m, pinv_mat, 11, 7);
+ puts("");
+}
diff --git a/linear_algebra/pseudo_inverse/qr_pseudo_inverse.c b/linear_algebra/pseudo_inverse/qr_pseudo_inverse.c
new file mode 100644
index 0000000000000000000000000000000000000000..cbb06419cb0377dde2b2f17f4c2f2cc92c360869
--- /dev/null
+++ b/linear_algebra/pseudo_inverse/qr_pseudo_inverse.c
@@ -0,0 +1,90 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief QR decomposition algorithms to compute the pseudo-inverse of a matrix.
+ *
+ * @details The computation of the pseudo-inverse is implemented using the Householder or
+ * the Givens algorithms.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+
+#include "matrix.h"
+#include "qr_common.h"
+#include "qr_givens.h"
+#include "qr_householder.h"
+#include "pseudo_inverse.h"
+
+int8_t qr_get_pinv(uint8_t m, uint8_t n, matrix_t A[m][n],
+ matrix_t pinv_A[n][m], enum QR_ALGORITHM algo)
+{
+ int8_t status;
+
+ switch (algo) {
+ case QR_Householder:
+ {
+ puts("Householder !!!!");
+ matrix_t Q[m][m];
+ matrix_t R[m][n];
+ matrix_copy(m, n, A, R);
+ status = qr_householder_decomp(m, n, R, m, Q, false);
+ matrix_t inv_R[n][m];
+ matrix_get_inv_upp_triang(m, n, R, inv_R);
+
+ // A+ = R^-1*Q'
+ matrix_in_place_transpose(m, Q);
+ matrix_mul(n, m, inv_R, m, m, Q, pinv_A);
+
+ break;
+ }
+
+ case QR_Givens:
+ {
+ puts("Givens !!!");
+
+ matrix_t Q[m][m];
+ matrix_t R[m][n];
+ matrix_copy(m, n, A, R);
+ status = qr_givens_decomp(m, n, R, m, Q, false);
+ matrix_t inv_R[n][m];
+ matrix_get_inv_upp_triang(m, n, R, inv_R);
+
+ // A+ = R^-1*Q'
+ matrix_in_place_transpose(m, Q);
+ matrix_mul(n, m, inv_R, m, m, Q, pinv_A);
+
+ break;
+ }
+
+ default:
+ {
+ matrix_t Q[m][m];
+ matrix_t R[m][n];
+ matrix_copy(m, n, A, R);
+ status = qr_householder_decomp(m, n, R, m, Q, false);
+ matrix_t inv_R[n][m];
+ matrix_get_inv_upp_triang(m, n, R, inv_R);
+
+ // A+ = R^-1*Q'
+ matrix_in_place_transpose(m, Q);
+ matrix_mul(n, m, inv_R, m, m, Q, pinv_A);
+
+ }
+ }
+
+ return status;
+}
diff --git a/linear_algebra/solve_linear_equations/doc.txt b/linear_algebra/solve_linear_equations/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..83fdda83ced14816bfd8f9cc75a69f3b89ee35ef
--- /dev/null
+++ b/linear_algebra/solve_linear_equations/doc.txt
@@ -0,0 +1,21 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup solve_linear_equations SOLVE_LINEAR_EQUATIONS
+ * @ingroup linear_algebra
+ *
+ * @brief Enables to solve systems of linear equations Ax = b for x.
+ *
+ * @details The user can select various algorithm such as the Moore--Penrose
+ * inverse, the Givens or the Householder algorithm for the
+ * QR-decomposition to solve the systems of linear equations.
+ *
+ * @author Zakaria Kasmi
+ */
diff --git a/linear_algebra/solve_linear_equations/include/solve.h b/linear_algebra/solve_linear_equations/include/solve.h
new file mode 100644
index 0000000000000000000000000000000000000000..ece6ed3bba7815a45aa7880f7ee96beadee7267c
--- /dev/null
+++ b/linear_algebra/solve_linear_equations/include/solve.h
@@ -0,0 +1,102 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup solve_linear_equations
+ * @{
+ *
+ * @file
+ * @brief Enables to solve systems of linear equations Ax = b for x.
+ * @details The user can select various algorithm such as the Moore--Penrose
+ * inverse, the Givens or the Householder algorithm for the QR-decomposition.
+ * The user can also choose the Gaussian Elimination with pivoting algorithm to solve
+ * the systems of linear equations.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef SOLVE_H_
+#define SOLVE_H_
+
+#include "matrix.h"
+#include "pseudo_inverse.h"
+
+/**
+ * @brief Solve an (m \f$\times\f$ n) linear system Ax = b by using the Moore--Penrose,
+ * Householder, or the Givens algorithm.
+ *
+ * @param[in] m row number of the matrix A.
+ * @param[in] n column number of the matrix A.
+ * @param[in] A[][] pointer to the matrix A.
+ * @param[in] b[] pointer to the vector b.
+ * @param[out] x_sol[] pointer to the solution vector.
+ * @param[in] algo specifies the algorithm to use (e.g. the Householder method).
+ *
+ * @return 1, if solving the linear equation system is successful.
+ * @return -1, if solving the linear equation system is not successful.
+ *
+ */
+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);
+
+/**
+ * @brief Solve an (m \f$\times\f$ n) linear system Ax = b, using the Householder
+ * algorithm.
+ *
+ * @param[in] m row number of the matrix A.
+ * @param[in] n column number of the matrix A.
+ * @param[in] A[][] pointer to the matrix A.
+ * @param[in] b[] pointer to the vector b.
+ * @param[out] x_sol[] pointer to the solution vector.
+ *
+ * @return 1, if solving the linear equation system is successful.
+ * @return -1, if solving the linear equation system is not successful.
+ *
+ */
+int8_t solve_householder(uint8_t m, uint8_t n, matrix_t A[][n],
+ matrix_t b[m],
+ matrix_t x_sol[n]);
+
+/**
+ * @brief Solve an (m \f$\times\f$ n) linear system Ax = b, using the Givens algorithm.
+ *
+ * @param[in] m row number of the matrix A.
+ * @param[in] n column number of the matrix A.
+ * @param[in] A[][] pointer to the matrix A.
+ * @param[in] b[] pointer to the vector b.
+ * @param[out] x_sol[] pointer to the solution vector.
+ *
+ * @return 1, if solving the linear equation system is successful.
+ * @return -1, if solving the linear equation system is not successful.
+ * @return -1, if the linear system is not solvable.
+ *
+ */
+int8_t solve_givens(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m],
+ matrix_t x_sol[n]);
+
+/**
+ * @brief Solve an (m \f$\times\f$ n) linear system Ax = b, using the Gaussian Elimination with
+ * pivoting algorithm.
+ *
+ * @param[in] m row number of the matrix A.
+ * @param[in] n column number of the matrix A.
+ * @param[in] A[][] pointer to the matrix A.
+ * @param[in] b[] pointer to the vector b.
+ * @param[out] x_sol[] pointer to the solution vector.
+ *
+ * @return 1, if solving the linear equation system is successful.
+ * @return -1, if solving the linear equation system is not successful.
+ * @return -2, if the linear system is not solvable.
+ *
+ */
+int8_t solve_lu_decomp(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m],
+ matrix_t x_sol[n]);
+#endif /* SOLVE_H_ */
diff --git a/linear_algebra/solve_linear_equations/solve.c b/linear_algebra/solve_linear_equations/solve.c
new file mode 100644
index 0000000000000000000000000000000000000000..44976490c243498187cdae8353a80096c2a6a318
--- /dev/null
+++ b/linear_algebra/solve_linear_equations/solve.c
@@ -0,0 +1,152 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ *
+ * @{
+ *
+ * @file
+ * @brief Enables to solve systems of linear equations Ax = b for x.
+ * @details The user can select various algorithm such as the Moore--Penrose
+ * inverse, the Givens or the Householder algorithm for the
+ * QR-decomposition to solve the systems of linear equations.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+
+#include "solve.h"
+#include "matrix.h"
+#include "qr_givens.h"
+#include "qr_common.h"
+#include "lu_decomp.h"
+#include "qr_householder.h"
+#include "moore_penrose_pseudo_inverse.h"
+
+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)
+{
+ int8_t status;
+
+ switch (algo) {
+ case Moore_Penrose:
+ {
+ //puts("Moore Penrose !!!!");
+ matrix_t pinv_A[n][m];
+ status = moore_penrose_get_pinv(m, n, A, pinv_A);
+ matrix_mul_vec(n, m, pinv_A, b, x_sol);
+ break;
+ }
+ case Householder:
+ //puts("Householder !!!!");
+ status = solve_householder(m, n, A, b, x_sol);
+ break;
+
+ case Givens:
+ //puts("Givens !!!");
+ status = solve_givens(m, n, A, b, x_sol);
+ break;
+
+ case Gauss:
+ //puts("Gauss !!!");
+ status = solve_lu_decomp(m, n, A, b, x_sol);
+ break;
+
+ default:
+ {
+ //puts("Default: Moore Penrose");
+ matrix_t pinv_A[n][m];
+ status = moore_penrose_get_pinv(m, n, A, pinv_A);
+ matrix_mul_vec(n, m, pinv_A, b, x_sol);
+ }
+
+ }
+
+ return status;
+}
+
+int8_t solve_householder(uint8_t m, uint8_t n, matrix_t A[][n],
+ matrix_t b[m],
+ matrix_t x_sol[n])
+{
+
+ matrix_t Q[m][n];
+ matrix_t qt_b[n];
+
+ if (m >= n) {
+ int8_t status;
+ status = qr_householder_decomp(m, n, A, n, Q, true);
+
+ /* qt_b = Q'*b --> Rx = Q'b (R is stored in A) */
+ matrix_trans_mul_vec(m, n, Q, m, b, qt_b);
+ qr_common_backward_subst(n, n, A, qt_b, x_sol);
+
+ return status;
+ }
+ else {
+ puts("[solve_householder]: The equation is not solvable !!!");
+
+ return -2;
+ }
+}
+
+int8_t solve_givens(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m],
+ matrix_t x_sol[n])
+{
+ matrix_t Q[m][n];
+ matrix_t qt_b[n];
+
+ if (m >= n) {
+ int8_t status;
+ status = qr_givens_decomp(m, n, A, n, Q, true);
+
+ /* qt_b = Q'*b --> Rx = Q'b (R is stored in A) */
+ matrix_trans_mul_vec(m, n, Q, m, b, qt_b);
+ qr_common_backward_subst(n, n, A, qt_b, x_sol);
+
+ return status;
+ }
+ else {
+ puts("[solve_givens]: The equation is not solvable !!!");
+ return -2;
+ }
+}
+
+int8_t solve_lu_decomp(uint8_t m, uint8_t n, matrix_t A[][n], matrix_t b[m],
+ matrix_t x_sol[n])
+{
+
+ if (m == n) {
+ //puts("solve_lu_decomp !!!");
+ matrix_t L[m][m];
+ matrix_t P[m][m];
+ lu_decomp(m, A, L, P);
+
+ /* lu_b = inv(L)*P*b */
+ matrix_t inv_L[m][m];
+ matrix_get_inv_low_triang(m, m, L, inv_L);
+ matrix_t tmp[m][m];
+ matrix_mul(m, m, inv_L, m, m, P, tmp);
+ matrix_t lu_b[m];
+ matrix_mul_vec(m, m, tmp, b, lu_b);
+
+ /* U*x = inv(L)*P*b = lu_b*/
+ qr_common_backward_subst(n, n, A, lu_b, x_sol);
+ return 1;
+ }
+
+ else {
+ puts(
+ "[solve_lu_decomp]: Only quadratic linear equation systems are supported !!!");
+ return -1;
+ }
+}
diff --git a/linear_algebra/utilities/combinatorics.c b/linear_algebra/utilities/combinatorics.c
new file mode 100644
index 0000000000000000000000000000000000000000..214322b97ee6df205b0dc6f9a442b7273dee8ae4
--- /dev/null
+++ b/linear_algebra/utilities/combinatorics.c
@@ -0,0 +1,77 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ *
+ * @{
+ *
+ * @file
+ * @brief Calculate possible \f$ \binom{n}{k} combinations \f$ without repetition in ascending
+ * order.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "combinatorics.h"
+#include "vector.h"
+
+uint8_t combinatorics_init(uint8_t n, uint8_t k, uint8_t comb_arr[])
+{
+ if (k > n) {
+ return (COMBI_ERROR);
+ }
+
+ if (k == 0) {
+ return (COMBI_EMPTY);
+ }
+
+
+ // Initialization of the combinatoric array with k-values.
+ for (uint32_t i = 0; i < k; i++) {
+ comb_arr[i] = i;
+ }
+
+ return (COMBI_SUCCESS);
+}
+
+uint8_t combinatorics_get_next_without_rep(uint8_t n, uint8_t k,
+ uint8_t comb_arr[])
+{
+
+ if (comb_arr[k - 1] < n - 1) {
+ comb_arr[k - 1] = comb_arr[k - 1] +1;
+ return (COMBI_SUCCESS);
+ }
+
+ int32_t i;
+ for (i = k - 2; i >= 0; i--) {
+ if (comb_arr[i] < n - k + i) {
+ break;
+ }
+ }
+
+ //Break, if comb_arr[0] == n - k
+ if (i < 0) {
+ return (COMBI_END);
+ }
+
+ comb_arr[i] = comb_arr[i] + 1;
+
+ while (i < k - 1) {
+ comb_arr[i + 1] = comb_arr[i] + 1;
+ i++;
+ }
+
+ return (COMBI_SUCCESS);
+}
diff --git a/linear_algebra/utilities/doc.txt b/linear_algebra/utilities/doc.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7754019330097308ed27ab6e8fa2bfa32bba0b1d
--- /dev/null
+++ b/linear_algebra/utilities/doc.txt
@@ -0,0 +1,20 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser
+ * General Public License v2.1. See the file LICENSE in the top level
+ * directory for more details.
+ */
+
+/**
+ * @defgroup utilities UTILITIES
+ * @ingroup linear_algebra
+ *
+ * @brief Utilities for linear algebra.
+ * @details Utility-functions are needed by linear algebra-module as well as
+ * other modules such as the position algorithm-module.
+ *
+ * @author Zakaria Kasmi
+ *
+ */
diff --git a/linear_algebra/utilities/include/combinatorics.h b/linear_algebra/utilities/include/combinatorics.h
new file mode 100644
index 0000000000000000000000000000000000000000..0a3ec4a9c896e44136ebb6f4218e3c17c59c43bd
--- /dev/null
+++ b/linear_algebra/utilities/include/combinatorics.h
@@ -0,0 +1,77 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup utilities
+ * @{
+ *
+ * @file
+ * @brief Calculate possible \f$ \binom{n}{k} combinations \f$ without repetition in ascending
+ * order.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef COMBINATORICS_H_
+#define COMBINATORICS_H_
+
+#include<stdint.h>
+
+/**
+ * Case of an error.
+ */
+#define COMBI_ERROR -1
+
+/**
+ * Case of an empty combination set.
+ */
+#define COMBI_EMPTY 0
+
+/**
+ * Case of successfully calculated combination set.
+ */
+#define COMBI_SUCCESS 1
+
+/**
+ * Case of completion of calculating combination sets.
+ */
+#define COMBI_END 2
+
+
+/**
+ * @brief Initialize the combinations generator.
+ *
+ * @param[in] n size of the set.
+ * @param[in] k size of the sub-set.
+ * @param[out] comb_arr[] pointer to the combination set.
+ *
+ * return @ref COMBI_ERROR, if k > n.
+ * return @ref COMBI_EMPTY, if k =0.
+ * return @ref COMBI_SUCCESS, if successful.
+ *
+ */
+uint8_t combinatorics_init(uint8_t n, uint8_t k, uint8_t comb_arr[]);
+
+/**
+ * @brief Generate the next combination.
+ *
+ * @param[in] n size of the set.
+ * @param[in] k size of the sub-set.
+ * @param[in, out] comb_arr[] pointer to the combination set.
+ *
+ * return @ref COMBI_END, if the last combination is generated.
+ * return @ref COMBI_SUCCESS, if successful.
+ *
+ */
+uint8_t combinatorics_get_next_without_rep(uint8_t n, uint8_t k, uint8_t comb_arr[]);
+
+
+#endif /*COMBINATORICS_H_ */
diff --git a/linear_algebra/utilities/include/norm_dist_rnd_generator.h b/linear_algebra/utilities/include/norm_dist_rnd_generator.h
new file mode 100644
index 0000000000000000000000000000000000000000..369c1a6e578f8f059942c90e1f9c6c330009103a
--- /dev/null
+++ b/linear_algebra/utilities/include/norm_dist_rnd_generator.h
@@ -0,0 +1,39 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup utilities
+ * @{
+ *
+ * @file
+ * @brief Generating normally distributed random numbers.
+ * @details The generation of normally distributed random numbers is implemented by using
+ * the Box--Muller method.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+// COMBINATORICS_H_
+
+#ifndef NORM_DIST_RND_GENERATOR_H_
+#define NORM_DIST_RND_GENERATOR_H_
+
+/*!
+ \def PI
+ Pi, the ratio of a circle's circumference to its diameter.
+ */
+#define PI 3.14159265
+
+double get_norm_distr_rand_num(double mean, double std_dev);
+double get_rand_num(int seed);
+
+
+#endif /* NORM_DIST_RND_GENERATOR_H_ */
diff --git a/linear_algebra/utilities/include/shell_sort.h b/linear_algebra/utilities/include/shell_sort.h
new file mode 100644
index 0000000000000000000000000000000000000000..6e0df2c8b3f20ba70ea25e63684a7af9c60db7f1
--- /dev/null
+++ b/linear_algebra/utilities/include/shell_sort.h
@@ -0,0 +1,49 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup utilities
+ * @{
+ *
+ * @file
+ * @brief Implement the Shell sort algorithm.
+ * @details The Shell sort algorithm is more convenient for devices with limited storage
+ * capacity.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef SHELL_SORT_H_
+#define SHELL_SORT_H_
+
+#include <stdint.h>
+#include "vector.h"
+
+/**
+ * @brief Sort a data set of integers by using the Shell sort algorithm.
+ *
+ * @param[in] array[] pointer to the data set.
+ * @param[in] length size of the data set.
+ *
+ */
+void int_shell_sort(int *array, int length);
+
+/**
+ * @brief Sort a data set of type utils_t by using the Shell sort algorithm.
+ *
+ * @param[in] arr[] pointer to the data set.
+ * @param[in] length size of the data set.
+ *
+ */
+void shell_sort(vector_t *arr, uint8_t length);
+
+
+#endif /* SHELL_SORT_H_ */
diff --git a/linear_algebra/utilities/include/shell_sort.h.bak b/linear_algebra/utilities/include/shell_sort.h.bak
new file mode 100644
index 0000000000000000000000000000000000000000..cf241333ed0cb14f3d7b3704fb6dd51ac2c26bc9
--- /dev/null
+++ b/linear_algebra/utilities/include/shell_sort.h.bak
@@ -0,0 +1,60 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup utils
+ * @{
+ *
+ * @file
+ * @brief Implement the Shell sort algorithm.
+ * @details The Shell sort algorithm is more convenient for devices with limited storage
+ * capacity.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef SHELL_SORT_H_
+#define SHELL_SORT_H_
+
+#include <stdint.h>
+#include "vector.h"
+
+/**
+ * @brief Sort a data set of integers by using the Shell sort algorithm.
+ *
+ * @param[in] array[] pointer to the data set.
+ * @param[in] length size of the data set.
+ *
+ */
+void int_shell_sort(int* array, int length);
+
+/**
+ * @brief Sort a data set of type utils_t by using the Shell sort algorithm.
+ *
+ * @param[in] arr[] pointer to the data set.
+ * @param[in] length size of the data set.
+ *
+ */
+void shell_sort(vector_t* arr, uint8_t length);
+
+/**
+ * @brief Sort a data set of type utils_t by using the Shell sort algorithm.
+ *
+ * @param[in] arr[] pointer to the data set.
+ * @param[in] length size of the data set.
+ *
+ */
+void double_shell_sort(vector_t* arr, uint8_t length);
+
+
+
+
+
+#endif /* SHELL_SORT_H_ */
diff --git a/linear_algebra/utilities/include/utils.h b/linear_algebra/utilities/include/utils.h
new file mode 100644
index 0000000000000000000000000000000000000000..9d6b25535096e14ebb603cd4d1b8dd01387c2a00
--- /dev/null
+++ b/linear_algebra/utilities/include/utils.h
@@ -0,0 +1,172 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup utilities
+ * @{
+ *
+ * @file
+ * @brief Utilities for linear algebra.
+ * @details Utility-functions are needed by linear algebra-module as well as
+ * other modules such as the position algorithm-module.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef UTILS_H_
+#define UTILS_H_
+
+#include <inttypes.h>
+#include "vector.h"
+
+/** @brief Define the pi-constant */
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+/**
+ * @brief Convert the angle from degrees to radians.
+ *
+ * @param[in] deg_angle angle in degrees.
+ *
+ * @return angle in radians.
+ *
+ */
+double utils_to_radian(double deg_angle);
+
+/**
+ * @brief Compute the sine of a variable in degrees.
+ * Calculate the sine of the variable deg_angle, which is expressed
+ * in degrees.
+ *
+ * @param[in] deg_angle angle in degrees.
+ *
+ * @return sine value.
+ *
+ */
+double utils_sind(double deg_angle);
+
+/**
+ * @brief Interchange the values of two variables of type uint8_t.
+ *
+ * @param[in] *a pointer to first variable.
+ * @param[in] *b pointer to second variable.
+ *
+ */
+void utils_swap(uint8_t *a, uint8_t *b);
+
+/**
+ * @brief Returns the greater of two real numbers.
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the greater of a and b.
+ *
+ */
+double utils_max(double a, double b);
+
+/**
+ * @brief Returns the smaller of two real numbers.
+ *
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the smaller of a and b.
+ *
+ */
+double utils_min(double a, double b);
+
+/**
+ * @brief Returns the greater of two numbers from type uint8_t.
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the greater of a and b.
+ *
+ */
+uint8_t utils_u8_max(uint8_t a, uint8_t b);
+
+/**
+ * @brief Returns the smaller of two numbers from type uint8_t.
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the smaller of a and b.
+ *
+ */
+uint8_t utils_u8_min(uint8_t a, uint8_t b);
+
+// Enables to use variable format string as well as argument lists
+// Fixing error: format not a string literal.
+
+/**
+ * @brief Print by using variable format string as well as argument lists.
+ * This function enables to print data by using a variable format
+ * string as well as argument list. Furthermore, it avoids the error:
+ * "format not a string literal", if printf is used.
+ *
+ * @param[in] *format_str format string.
+ * @param[in] ... argument list.
+ *
+ */
+void utils_printf(char *format_str, ...);
+
+/**
+ * @brief Compute the mean value of a data set.
+ *
+ * @param[in] arr_size size of the data set.
+ * @param[in] in_arr[] pointer to the data set.
+ *
+ * @return the mean value of the data set.
+ *
+ */
+double utils_mean(uint8_t arr_size, vector_t in_arr[]);
+
+/**
+ * @brief Compute the moving average of a data set.
+ *
+ * @param[in] arr_size size of the data set.
+ * @param[in] in_arr[] pointer to the data set.
+ * @param[in] window_size window size.
+ * @param[out] out_arr pointer to the values of the moving average.
+ *
+ */
+void utils_moving_average(uint8_t arr_size, vector_t in_arr[],
+ uint8_t window_size,
+ vector_t out_arr[]);
+
+/**
+ * @brief Compute the median of a finite array of numbers.
+ *
+ * @param[in] arr[] pointer to the data set.
+ * @param[in] length size of the data set.
+ *
+ * @return the median value of the data set.
+ *
+ */
+double utils_get_median(vector_t arr[], uint8_t length);
+
+/**
+ * @brief Compute the square root without under/overflow.
+ *
+ * @param[in] x first value.
+ * @param[in] y second value.
+ *
+ * @return save square root.
+ *
+ */
+double utils_get_save_square_root(double x, double y);
+
+#endif /* UTILS_H_ */
diff --git a/linear_algebra/utilities/include/utils.h.bak b/linear_algebra/utilities/include/utils.h.bak
new file mode 100644
index 0000000000000000000000000000000000000000..4316862c16abaaee8717f2fdb3800d83c2dcfac3
--- /dev/null
+++ b/linear_algebra/utilities/include/utils.h.bak
@@ -0,0 +1,176 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @ingroup utils
+ * @{
+ *
+ * @file
+ * @brief Utilities for linear algebra.
+ * @details Utility-functions are needed by linear algebra-module as well as
+ * other modules such as the position algorithm-module.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#ifndef UTILS_H_
+#define UTILS_H_
+
+#include <inttypes.h>
+
+/** @brief Define the data type */
+#define utils_t double
+
+#ifndef utils_t
+#define utils_t double
+#endif
+
+/** @brief Define the pi-constant */
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
+
+/**
+ * @brief Convert the angle from degrees to radians.
+ *
+ * @param[in] deg_angle angle in degrees.
+ *
+ * @return angle in radians.
+ *
+ */
+double utils_to_radian(double deg_angle);
+
+/**
+ * @brief Compute the sine of a variable in degrees.
+ * Calculate the sine of the variable deg_angle, which is expressed
+ * in degrees.
+ *
+ * @param[in] deg_angle angle in degrees.
+ *
+ * @return sine value.
+ *
+ */
+double utils_sind(double deg_angle);
+
+/**
+ * @brief Interchange the values of two variables of type uint8_t.
+ *
+ * @param[in] *a pointer to first variable.
+ * @param[in] *b pointer to second variable.
+ *
+ */
+void utils_swap(uint8_t *a, uint8_t *b);
+
+/**
+ * @brief Returns the greater of two real numbers.
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the greater of a and b.
+ *
+ */
+double utils_max(double a, double b);
+
+/**
+ * @brief Returns the smaller of two real numbers.
+ *
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the smaller of a and b.
+ *
+ */
+double utils_min(double a, double b);
+
+/**
+ * @brief Returns the greater of two numbers from type uint8_t.
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the greater of a and b.
+ *
+ */
+uint8_t utils_u8_max(uint8_t a, uint8_t b);
+
+/**
+ * @brief Returns the smaller of two numbers from type uint8_t.
+ *
+ * @param[in] a the first value to compare.
+ * @param[in] b the second value to compare.
+ *
+ * @return the smaller of a and b.
+ *
+ */
+uint8_t utils_u8_min(uint8_t a, uint8_t b);
+
+// Enables to use variable format string as well as argument lists
+// Fixing error: format not a string literal.
+
+/**
+ * @brief Print by using variable format string as well as argument lists.
+ * This function enables to print data by using a variable format
+ * string as well as argument list. Furthermore, it avoids the error:
+ * "format not a string literal", if printf is used.
+ *
+ * @param[in] *format_str format string.
+ * @param[in] ... argument list.
+ *
+ */
+void utils_printf(char *format_str, ...);
+
+/**
+ * @brief Compute the mean value of a data set.
+ *
+ * @param[in] arr_size size of the data set.
+ * @param[in] in_arr[] pointer to the data set.
+ *
+ * @return the mean value of the data set.
+ *
+ */
+double utils_mean(uint8_t arr_size, double in_arr[]);
+
+/**
+ * @brief Compute the moving average of a data set.
+ *
+ * @param[in] arr_size size of the data set.
+ * @param[in] in_arr[] pointer to the data set.
+ * @param[in] window_size window size.
+ * @param[out] out_arr pointer to the values of the moving average.
+ *
+ */
+void utils_moving_average(uint8_t arr_size, double in_arr[], uint8_t window_size, double out_arr[]);
+
+
+/**
+ * @brief Compute the median of a finite array of numbers.
+ *
+ * @param[in] arr[] pointer to the data set.
+ * @param[in] length size of the data set.
+ *
+ * @return the median value of the data set.
+ *
+ */
+double utils_get_median(double arr[], uint8_t length);
+
+/**
+ * @brief Compute the square root without under/overflow.
+ *
+ * @param[in] x first value.
+ * @param[in] y second value.
+ *
+ * @return save square root.
+ *
+ */
+double utils_get_save_square_root(double x, double y);
+
+#endif /* UTILS_H_ */
diff --git a/linear_algebra/utilities/norm_dist_rnd_generator.c b/linear_algebra/utilities/norm_dist_rnd_generator.c
new file mode 100644
index 0000000000000000000000000000000000000000..d7766a767b4cc7ff59b673c81f9b811af5b36592
--- /dev/null
+++ b/linear_algebra/utilities/norm_dist_rnd_generator.c
@@ -0,0 +1,108 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ *
+ * @{
+ *
+ * @file
+ * @brief Generating normally distributed random numbers.
+ * @details The generation of normally distributed random numbers is implemented by using
+ * the Box--Muller method.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdlib.h>
+#include <math.h>
+
+#include "norm_dist_rnd_generator.h"
+
+
+
+/**
+ * @brief Get a normally distributed random number by applying the Box--Muller method.
+ *
+ * @param[in] mean_val mean value.
+ * @param[in] std_dev_val standard deviation.
+ *
+ * @return normally distributed random number.
+ *
+ */
+double get_norm_distr_rand_num(double mean_val, double std_dev_val)
+{
+ double u, r, theta;
+ double x;
+ double adj_norm_rv;
+
+ // Generate a random value of u
+ u = 0.0;
+ while (u == 0.0) {
+ u = get_rand_num(0);
+ }
+
+ r = sqrt(-2.0 * log(u));
+
+ // Generate a random value of theta
+ theta = 0.0;
+ while (theta == 0.0) {
+ theta = 2.0 * PI * get_rand_num(0);
+ }
+
+ x = r * cos(theta);
+
+ adj_norm_rv = (x * std_dev_val) + mean_val;
+
+ return (adj_norm_rv);
+}
+
+/**
+ * @brief Generate uniform (0.0, 1.0) random numbers by using the Linear Congruential
+ * Generator (LGC) algorithm.
+ *
+ * @note The LGC is implemented based on the following book:
+ * R. Jain, "The Art of Computer Systems Performance Analysis,"
+ * John Wiley & Sons, 1991.
+ *
+ * @param[in] initial_seed_val initial seed value.
+ *
+ * @return random number between 0.0 and 1.0.
+ *
+ */
+double get_rand_num(int initial_seed_val)
+{
+ const long a = 16807;
+ const long m = 2147483647;
+ const long q = 127773;
+ const long r = 2836;
+ static long x;
+ long x_div_q;
+ long x_mod_q;
+ long new_x;
+
+ // Setting of the seed value
+ if (initial_seed_val > 0) {
+ x = initial_seed_val;
+ return (0.0);
+ }
+
+ x_div_q = x / q;
+ x_mod_q = x % q;
+ new_x = (a * x_mod_q) - (r * x_div_q);
+ if (new_x > 0) {
+ x = new_x;
+ }
+ else {
+ x = new_x + m;
+ }
+
+ return ((double)x / m);
+}
diff --git a/linear_algebra/utilities/shell_sort.c b/linear_algebra/utilities/shell_sort.c
new file mode 100644
index 0000000000000000000000000000000000000000..c64d662f49680cc28134c3b9f405e15311e88b29
--- /dev/null
+++ b/linear_algebra/utilities/shell_sort.c
@@ -0,0 +1,63 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Implement the Shell sort algorithm.
+ * @details The Shell sort algorithm is more convenient for devices with limited storage
+ * capacity.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <stdio.h>
+#include "utils.h"
+#include "vector.h"
+
+/* Shell Sort Program */
+void int_shell_sort(int *array, int length)
+{
+ int n = length;
+ int i, j, gap, temp;
+
+ for (gap = n / 2; gap > 0; gap /= 2) {
+ for (i = gap; i < n; i++) {
+ temp = array[i];
+ for (j = i; j >= gap; j -= gap) {
+ if (temp < array[j - gap]) {
+ array[j] = array[j - gap];
+ }
+ else {
+ break;
+ }
+ }
+ array[j] = temp;
+ }
+ }
+}
+
+void shell_sort(vector_t *arr, uint8_t length)
+{
+ int32_t j, k;
+ vector_t temp;
+
+ for (int32_t i = length / 2; i > 0; i = i / 2) {
+ for (j = i; j < length; j++) {
+ temp = arr[j];
+ for (k = j - 1; k >= 0 && arr[k] > temp; k--) {
+ arr[k + 1] = arr[k];
+ }
+ arr[k + 1] = temp;
+ }
+ }
+}
diff --git a/linear_algebra/utilities/utils.c b/linear_algebra/utilities/utils.c
new file mode 100644
index 0000000000000000000000000000000000000000..b0f6673eda562d584af06ee21adf30c9e2b98516
--- /dev/null
+++ b/linear_algebra/utilities/utils.c
@@ -0,0 +1,185 @@
+/*
+ * Copyright (C) 2020 Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ * 2020 Freie Universität Berlin
+ *
+ * This file is subject to the terms and conditions of the GNU Lesser General
+ * Public License v2.1. See the file LICENSE in the top level directory for more
+ * details.
+ */
+
+/**
+ * @{
+ *
+ * @file
+ * @brief Utilities for linear algebra.
+ * Utility-functions are needed by the linear algebra-module as
+ * well as other modules such as the position algorithm-module.
+ *
+ * @author Zakaria Kasmi <zkasmi@inf.fu-berlin.de>
+ *
+ * @}
+ */
+
+#include <math.h>
+#include <stdarg.h>
+#include <stdio.h>
+
+#include "utils.h"
+#include "shell_sort.h"
+
+// from degree to radian
+double utils_to_radian(double deg_angle)
+{
+ double radians_angle = deg_angle * M_PI / 180.0;
+
+ return radians_angle;
+}
+
+//Sine of argument in degrees
+double utils_sind(double deg_angle)
+{
+ double radians_angle = deg_angle * M_PI / 180.0;
+
+ return sin(radians_angle);
+}
+
+void utils_swap(uint8_t *a, uint8_t *b)
+{
+ uint8_t tmp;
+
+ tmp = *a;
+ *a = *b;
+ *b = tmp;
+}
+
+double utils_max(double a, double b)
+{
+ if (a < b) {
+ return b;
+ }
+ else {
+ return a;
+ }
+}
+
+double utils_min(double a, double b)
+{
+ if (a < b) {
+ return a;
+ }
+ else {
+ return b;
+ }
+}
+
+uint8_t utils_u8_max(uint8_t a, uint8_t b)
+{
+ if (a < b) {
+ return b;
+ }
+ else {
+ return a;
+ }
+}
+
+uint8_t utils_u8_min(uint8_t a, uint8_t b)
+{
+ if (a < b) {
+ return a;
+ }
+ else {
+ return b;
+ }
+}
+
+// Enables to use variable format strings as well as argument lists
+void utils_printf(char *format_str, ...)
+{
+ va_list args;
+
+ va_start(args, format_str);
+ vprintf(format_str, args);
+ va_end(args);
+}
+
+double utils_mean(uint8_t arr_size, vector_t in_arr[])
+{
+ double mean = 0.0;
+
+ if (arr_size <= 0) {
+ puts("Not valid input array size !!!");
+ return mean;
+ }
+
+ for (uint8_t i = 0; i < arr_size; i++) {
+ mean += in_arr[i];
+ }
+
+ return mean / arr_size;
+}
+
+void utils_moving_average(uint8_t arr_size, vector_t in_arr[],
+ uint8_t window_size,
+ vector_t out_arr[])
+{
+ double sum = 0.0;
+ double trail_value = 0;
+ uint8_t pos = 0;
+
+ if (arr_size <= 0) {
+ puts("Not valid input array size !!!");
+ return;
+ }
+
+ if (window_size <= 0) {
+ puts("Not valid window size !!!");
+ return;
+ }
+
+ for (uint8_t i = 0; i < arr_size; i++) {
+ sum = sum - trail_value + in_arr[i];
+ out_arr[i] = sum / window_size;
+
+ pos++;
+ if (pos >= window_size) {
+
+ trail_value = in_arr[pos - window_size];
+
+ }
+ }
+}
+
+double utils_get_median(vector_t arr[], uint8_t length)
+{
+ int middle = length / 2;
+
+ shell_sort(arr, length);
+
+ if (length % 2 == 1) {
+ return arr[middle];
+ }
+ else {
+ return (arr[middle - 1] + arr[middle]) / 2.0;
+ }
+
+}
+
+/** sqrt(a^2 + b^2) without under/overflow. **/
+double utils_get_save_square_root(double x, double y)
+{
+ double sqr_root;
+
+ if (fabs(x) > fabs(y)) {
+ sqr_root = y / x;
+ sqr_root = fabs(x) * sqrt(1 + sqr_root * sqr_root);
+ }
+ else if (y != 0) {
+ sqr_root = x / y;
+ sqr_root = fabs(y) * sqrt(1 + sqr_root * sqr_root);
+ }
+ else {
+ sqr_root = 0.0;
+ }
+
+ return sqr_root;
+}