1 /* ---------------------------------------------------------------------- 2 * Project: CMSIS DSP Library 3 * Title: arm_mat_solve_lower_triangular_f64.c 4 * Description: Solve linear system LT X = A with LT lower triangular matrix 5 * 6 * $Date: 23 April 2021 7 * $Revision: V1.9.0 8 * 9 * Target Processor: Cortex-M and Cortex-A cores 10 * -------------------------------------------------------------------- */ 11 /* 12 * Copyright (C) 2010-2021 ARM Limited or its affiliates. All rights reserved. 13 * 14 * SPDX-License-Identifier: Apache-2.0 15 * 16 * Licensed under the Apache License, Version 2.0 (the License); you may 17 * not use this file except in compliance with the License. 18 * You may obtain a copy of the License at 19 * 20 * www.apache.org/licenses/LICENSE-2.0 21 * 22 * Unless required by applicable law or agreed to in writing, software 23 * distributed under the License is distributed on an AS IS BASIS, WITHOUT 24 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 25 * See the License for the specific language governing permissions and 26 * limitations under the License. 27 */ 28 29 #include "dsp/matrix_functions.h" 30 31 /** 32 @ingroup groupMatrix 33 */ 34 35 36 /** 37 @addtogroup MatrixInv 38 @{ 39 */ 40 41 42 /** 43 * @brief Solve LT . X = A where LT is a lower triangular matrix 44 * @param[in] lt The lower triangular matrix 45 * @param[in] a The matrix a 46 * @param[out] dst The solution X of LT . X = A 47 * @return The function returns ARM_MATH_SINGULAR, if the system can't be solved. 48 */ arm_mat_solve_lower_triangular_f64(const arm_matrix_instance_f64 * lt,const arm_matrix_instance_f64 * a,arm_matrix_instance_f64 * dst)49 arm_status arm_mat_solve_lower_triangular_f64( 50 const arm_matrix_instance_f64 * lt, 51 const arm_matrix_instance_f64 * a, 52 arm_matrix_instance_f64 * dst) 53 { 54 arm_status status; /* status of matrix inverse */ 55 56 57 #ifdef ARM_MATH_MATRIX_CHECK 58 59 /* Check for matrix mismatch condition */ 60 if ((lt->numRows != lt->numCols) || 61 (a->numRows != a->numCols) || 62 (lt->numRows != a->numRows) ) 63 { 64 /* Set status as ARM_MATH_SIZE_MISMATCH */ 65 status = ARM_MATH_SIZE_MISMATCH; 66 } 67 else 68 69 #endif /* #ifdef ARM_MATH_MATRIX_CHECK */ 70 71 { 72 /* a1 b1 c1 x1 = a1 73 b2 c2 x2 a2 74 c3 x3 a3 75 76 x3 = a3 / c3 77 x2 = (a2 - c2 x3) / b2 78 79 */ 80 int i,j,k,n; 81 82 n = dst->numRows; 83 84 float64_t *pX = dst->pData; 85 float64_t *pLT = lt->pData; 86 float64_t *pA = a->pData; 87 88 float64_t *lt_row; 89 float64_t *a_col; 90 91 for(j=0; j < n; j ++) 92 { 93 a_col = &pA[j]; 94 95 for(i=0; i < n ; i++) 96 { 97 lt_row = &pLT[n*i]; 98 99 float64_t tmp=a_col[i * n]; 100 101 for(k=0; k < i; k++) 102 { 103 tmp -= lt_row[k] * pX[n*k+j]; 104 } 105 106 if (lt_row[i]==0.0f) 107 { 108 return(ARM_MATH_SINGULAR); 109 } 110 tmp = tmp / lt_row[i]; 111 pX[i*n+j] = tmp; 112 } 113 114 } 115 status = ARM_MATH_SUCCESS; 116 117 } 118 119 /* Return to application */ 120 return (status); 121 } 122 /** 123 @} end of MatrixInv group 124 */ 125