1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_mat_cholesky_f64.c
4 * Description: Floating-point Cholesky decomposition
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 @addtogroup MatrixChol
37 @{
38 */
39
40 /**
41 * @brief Floating-point Cholesky decomposition of positive-definite matrix.
42 * @param[in] pSrc points to the instance of the input floating-point matrix structure.
43 * @param[out] pDst points to the instance of the output floating-point matrix structure.
44 * @return The function returns ARM_MATH_SIZE_MISMATCH, if the dimensions do not match.
45 * @return execution status
46 - \ref ARM_MATH_SUCCESS : Operation successful
47 - \ref ARM_MATH_SIZE_MISMATCH : Matrix size check failed
48 - \ref ARM_MATH_DECOMPOSITION_FAILURE : Input matrix cannot be decomposed
49 * @par
50 * If the matrix is ill conditioned or only semi-definite, then it is better using the LDL^t decomposition.
51 * The decomposition of A is returning a lower triangular matrix U such that A = U U^t
52 */
53
54
arm_mat_cholesky_f64(const arm_matrix_instance_f64 * pSrc,arm_matrix_instance_f64 * pDst)55 arm_status arm_mat_cholesky_f64(
56 const arm_matrix_instance_f64 * pSrc,
57 arm_matrix_instance_f64 * pDst)
58 {
59
60 arm_status status; /* status of matrix inverse */
61
62
63 #ifdef ARM_MATH_MATRIX_CHECK
64
65 /* Check for matrix mismatch condition */
66 if ((pSrc->numRows != pSrc->numCols) ||
67 (pDst->numRows != pDst->numCols) ||
68 (pSrc->numRows != pDst->numRows) )
69 {
70 /* Set status as ARM_MATH_SIZE_MISMATCH */
71 status = ARM_MATH_SIZE_MISMATCH;
72 }
73 else
74
75 #endif /* #ifdef ARM_MATH_MATRIX_CHECK */
76
77 {
78 int i,j,k;
79 int n = pSrc->numRows;
80 float64_t invSqrtVj;
81 float64_t *pA,*pG;
82
83 pA = pSrc->pData;
84 pG = pDst->pData;
85
86
87 for(i=0 ; i < n ; i++)
88 {
89 for(j=i ; j < n ; j++)
90 {
91 pG[j * n + i] = pA[j * n + i];
92
93 for(k=0; k < i ; k++)
94 {
95 pG[j * n + i] = pG[j * n + i] - pG[i * n + k] * pG[j * n + k];
96 }
97 }
98
99 if (pG[i * n + i] <= 0.0f)
100 {
101 return(ARM_MATH_DECOMPOSITION_FAILURE);
102 }
103
104 invSqrtVj = 1.0/sqrt(pG[i * n + i]);
105 for(j=i ; j < n ; j++)
106 {
107 pG[j * n + i] = pG[j * n + i] * invSqrtVj ;
108 }
109 }
110
111 status = ARM_MATH_SUCCESS;
112
113 }
114
115
116 /* Return to application */
117 return (status);
118 }
119
120 /**
121 @} end of MatrixChol group
122 */
123
