1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_mat_cholesky_f32.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 @defgroup MatrixChol Cholesky and LDLT decompositions
37
38 Computes the Cholesky or LDL^t decomposition of a matrix.
39
40
41 If the input matrix does not have a decomposition, then the
42 algorithm terminates and returns error status ARM_MATH_DECOMPOSITION_FAILURE.
43 */
44
45 /**
46 @addtogroup MatrixChol
47 @{
48 */
49
50 /**
51 * @brief Floating-point Cholesky decomposition of positive-definite matrix.
52 * @param[in] pSrc points to the instance of the input floating-point matrix structure.
53 * @param[out] pDst points to the instance of the output floating-point matrix structure.
54 * @return The function returns ARM_MATH_SIZE_MISMATCH, if the dimensions do not match.
55 * @return execution status
56 - \ref ARM_MATH_SUCCESS : Operation successful
57 - \ref ARM_MATH_SIZE_MISMATCH : Matrix size check failed
58 - \ref ARM_MATH_DECOMPOSITION_FAILURE : Input matrix cannot be decomposed
59 * @par
60 * If the matrix is ill conditioned or only semi-definite, then it is better using the LDL^t decomposition.
61 * The decomposition of A is returning a lower triangular matrix U such that A = U U^t
62 */
63
64 #if defined(ARM_MATH_MVEF) && !defined(ARM_MATH_AUTOVECTORIZE)
65
66 #include "arm_helium_utils.h"
67
arm_mat_cholesky_f32(const arm_matrix_instance_f32 * pSrc,arm_matrix_instance_f32 * pDst)68 arm_status arm_mat_cholesky_f32(
69 const arm_matrix_instance_f32 * pSrc,
70 arm_matrix_instance_f32 * pDst)
71 {
72
73 arm_status status; /* status of matrix inverse */
74
75
76 #ifdef ARM_MATH_MATRIX_CHECK
77
78 /* Check for matrix mismatch condition */
79 if ((pSrc->numRows != pSrc->numCols) ||
80 (pDst->numRows != pDst->numCols) ||
81 (pSrc->numRows != pDst->numRows) )
82 {
83 /* Set status as ARM_MATH_SIZE_MISMATCH */
84 status = ARM_MATH_SIZE_MISMATCH;
85 }
86 else
87
88 #endif /* #ifdef ARM_MATH_MATRIX_CHECK */
89
90 {
91 int i,j,k;
92 int n = pSrc->numRows;
93 float32_t invSqrtVj;
94 float32_t *pA,*pG;
95 int kCnt;
96
97 mve_pred16_t p0;
98
99 f32x4_t acc, acc0, acc1, acc2, acc3;
100 f32x4_t vecGi;
101 f32x4_t vecGj,vecGj0,vecGj1,vecGj2,vecGj3;
102
103
104 pA = pSrc->pData;
105 pG = pDst->pData;
106
107 for(i=0 ;i < n ; i++)
108 {
109 for(j=i ; j+3 < n ; j+=4)
110 {
111 pG[(j + 0) * n + i] = pA[(j + 0) * n + i];
112 pG[(j + 1) * n + i] = pA[(j + 1) * n + i];
113 pG[(j + 2) * n + i] = pA[(j + 2) * n + i];
114 pG[(j + 3) * n + i] = pA[(j + 3) * n + i];
115
116 kCnt = i;
117 acc0 = vdupq_n_f32(0.0f);
118 acc1 = vdupq_n_f32(0.0f);
119 acc2 = vdupq_n_f32(0.0f);
120 acc3 = vdupq_n_f32(0.0f);
121
122 for(k=0; k < i ; k+=4)
123 {
124 p0 = vctp32q(kCnt);
125
126 vecGi=vldrwq_z_f32(&pG[i * n + k],p0);
127
128 vecGj0=vldrwq_z_f32(&pG[(j + 0) * n + k],p0);
129 vecGj1=vldrwq_z_f32(&pG[(j + 1) * n + k],p0);
130 vecGj2=vldrwq_z_f32(&pG[(j + 2) * n + k],p0);
131 vecGj3=vldrwq_z_f32(&pG[(j + 3) * n + k],p0);
132
133 acc0 = vfmaq_m(acc0, vecGi, vecGj0, p0);
134 acc1 = vfmaq_m(acc1, vecGi, vecGj1, p0);
135 acc2 = vfmaq_m(acc2, vecGi, vecGj2, p0);
136 acc3 = vfmaq_m(acc3, vecGi, vecGj3, p0);
137
138 kCnt -= 4;
139 }
140 pG[(j + 0) * n + i] -= vecAddAcrossF32Mve(acc0);
141 pG[(j + 1) * n + i] -= vecAddAcrossF32Mve(acc1);
142 pG[(j + 2) * n + i] -= vecAddAcrossF32Mve(acc2);
143 pG[(j + 3) * n + i] -= vecAddAcrossF32Mve(acc3);
144 }
145
146 for(; j < n ; j++)
147 {
148 pG[j * n + i] = pA[j * n + i];
149
150 kCnt = i;
151 acc = vdupq_n_f32(0.0f);
152
153 for(k=0; k < i ; k+=4)
154 {
155 p0 = vctp32q(kCnt);
156
157 vecGi=vldrwq_z_f32(&pG[i * n + k],p0);
158 vecGj=vldrwq_z_f32(&pG[j * n + k],p0);
159
160 acc = vfmaq_m(acc, vecGi, vecGj,p0);
161
162 kCnt -= 4;
163 }
164 pG[j * n + i] -= vecAddAcrossF32Mve(acc);
165 }
166
167 if (pG[i * n + i] <= 0.0f)
168 {
169 return(ARM_MATH_DECOMPOSITION_FAILURE);
170 }
171
172 invSqrtVj = 1.0f/sqrtf(pG[i * n + i]);
173 for(j=i; j < n ; j++)
174 {
175 pG[j * n + i] = pG[j * n + i] * invSqrtVj ;
176 }
177 }
178
179 status = ARM_MATH_SUCCESS;
180
181 }
182
183
184 /* Return to application */
185 return (status);
186 }
187
188 #else
189 #if defined(ARM_MATH_NEON) && !defined(ARM_MATH_AUTOVECTORIZE)
190
arm_mat_cholesky_f32(const arm_matrix_instance_f32 * pSrc,arm_matrix_instance_f32 * pDst)191 arm_status arm_mat_cholesky_f32(
192 const arm_matrix_instance_f32 * pSrc,
193 arm_matrix_instance_f32 * pDst)
194 {
195
196 arm_status status; /* status of matrix inverse */
197
198
199 #ifdef ARM_MATH_MATRIX_CHECK
200
201 /* Check for matrix mismatch condition */
202 if ((pSrc->numRows != pSrc->numCols) ||
203 (pDst->numRows != pDst->numCols) ||
204 (pSrc->numRows != pDst->numRows) )
205 {
206 /* Set status as ARM_MATH_SIZE_MISMATCH */
207 status = ARM_MATH_SIZE_MISMATCH;
208 }
209 else
210
211 #endif /* #ifdef ARM_MATH_MATRIX_CHECK */
212
213 {
214 int i,j,k;
215 int n = pSrc->numRows;
216 float32_t invSqrtVj;
217 float32_t *pA,*pG;
218 int kCnt;
219
220
221 f32x4_t acc, acc0, acc1, acc2, acc3;
222 f32x4_t vecGi;
223 f32x4_t vecGj,vecGj0,vecGj1,vecGj2,vecGj3;
224 f32x2_t tmp = vdup_n_f32(0);
225 float32_t sum=0.0f;
226 float32_t sum0=0.0f,sum1=0.0f,sum2=0.0f,sum3=0.0f;
227
228
229 pA = pSrc->pData;
230 pG = pDst->pData;
231
232 for(i=0 ;i < n ; i++)
233 {
234 for(j=i ; j+3 < n ; j+=4)
235 {
236 pG[(j + 0) * n + i] = pA[(j + 0) * n + i];
237 pG[(j + 1) * n + i] = pA[(j + 1) * n + i];
238 pG[(j + 2) * n + i] = pA[(j + 2) * n + i];
239 pG[(j + 3) * n + i] = pA[(j + 3) * n + i];
240
241 acc0 = vdupq_n_f32(0.0f);
242 acc1 = vdupq_n_f32(0.0f);
243 acc2 = vdupq_n_f32(0.0f);
244 acc3 = vdupq_n_f32(0.0f);
245
246 kCnt = i >> 2;
247 k=0;
248 while(kCnt > 0)
249 {
250
251 vecGi=vld1q_f32(&pG[i * n + k]);
252
253 vecGj0=vld1q_f32(&pG[(j + 0) * n + k]);
254 vecGj1=vld1q_f32(&pG[(j + 1) * n + k]);
255 vecGj2=vld1q_f32(&pG[(j + 2) * n + k]);
256 vecGj3=vld1q_f32(&pG[(j + 3) * n + k]);
257
258 acc0 = vfmaq_f32(acc0, vecGi, vecGj0);
259 acc1 = vfmaq_f32(acc1, vecGi, vecGj1);
260 acc2 = vfmaq_f32(acc2, vecGi, vecGj2);
261 acc3 = vfmaq_f32(acc3, vecGi, vecGj3);
262
263 kCnt--;
264 k+=4;
265 }
266
267 #if __aarch64__
268 sum0 = vpadds_f32(vpadd_f32(vget_low_f32(acc0), vget_high_f32(acc0)));
269 sum1 = vpadds_f32(vpadd_f32(vget_low_f32(acc1), vget_high_f32(acc1)));
270 sum2 = vpadds_f32(vpadd_f32(vget_low_f32(acc2), vget_high_f32(acc2)));
271 sum3 = vpadds_f32(vpadd_f32(vget_low_f32(acc3), vget_high_f32(acc3)));
272
273 #else
274 tmp = vpadd_f32(vget_low_f32(acc0), vget_high_f32(acc0));
275 sum0 = vget_lane_f32(tmp, 0) + vget_lane_f32(tmp, 1);
276
277 tmp = vpadd_f32(vget_low_f32(acc1), vget_high_f32(acc1));
278 sum1 = vget_lane_f32(tmp, 0) + vget_lane_f32(tmp, 1);
279
280 tmp = vpadd_f32(vget_low_f32(acc2), vget_high_f32(acc2));
281 sum2 = vget_lane_f32(tmp, 0) + vget_lane_f32(tmp, 1);
282
283 tmp = vpadd_f32(vget_low_f32(acc3), vget_high_f32(acc3));
284 sum3 = vget_lane_f32(tmp, 0) + vget_lane_f32(tmp, 1);
285 #endif
286
287 kCnt = i & 3;
288 while(kCnt > 0)
289 {
290
291 sum0 = sum0 + pG[i * n + k] * pG[(j + 0) * n + k];
292 sum1 = sum1 + pG[i * n + k] * pG[(j + 1) * n + k];
293 sum2 = sum2 + pG[i * n + k] * pG[(j + 2) * n + k];
294 sum3 = sum3 + pG[i * n + k] * pG[(j + 3) * n + k];
295 kCnt--;
296 k++;
297 }
298
299 pG[(j + 0) * n + i] -= sum0;
300 pG[(j + 1) * n + i] -= sum1;
301 pG[(j + 2) * n + i] -= sum2;
302 pG[(j + 3) * n + i] -= sum3;
303 }
304
305 for(; j < n ; j++)
306 {
307 pG[j * n + i] = pA[j * n + i];
308
309 acc = vdupq_n_f32(0.0f);
310
311 kCnt = i >> 2;
312 k=0;
313 while(kCnt > 0)
314 {
315
316 vecGi=vld1q_f32(&pG[i * n + k]);
317 vecGj=vld1q_f32(&pG[j * n + k]);
318
319 acc = vfmaq_f32(acc, vecGi, vecGj);
320
321 kCnt--;
322 k+=4;
323 }
324
325 #if __aarch64__
326 sum = vpadds_f32(vpadd_f32(vget_low_f32(acc), vget_high_f32(acc)));
327 #else
328 tmp = vpadd_f32(vget_low_f32(acc), vget_high_f32(acc));
329 sum = vget_lane_f32(tmp, 0) + vget_lane_f32(tmp, 1);
330 #endif
331
332 kCnt = i & 3;
333 while(kCnt > 0)
334 {
335 sum = sum + pG[i * n + k] * pG[(j + 0) * n + k];
336
337
338 kCnt--;
339 k++;
340 }
341
342 pG[j * n + i] -= sum;
343 }
344
345 if (pG[i * n + i] <= 0.0f)
346 {
347 return(ARM_MATH_DECOMPOSITION_FAILURE);
348 }
349
350 invSqrtVj = 1.0f/sqrtf(pG[i * n + i]);
351 for(j=i; j < n ; j++)
352 {
353 pG[j * n + i] = pG[j * n + i] * invSqrtVj ;
354 }
355 }
356
357 status = ARM_MATH_SUCCESS;
358
359 }
360
361
362 /* Return to application */
363 return (status);
364 }
365
366 #else
arm_mat_cholesky_f32(const arm_matrix_instance_f32 * pSrc,arm_matrix_instance_f32 * pDst)367 arm_status arm_mat_cholesky_f32(
368 const arm_matrix_instance_f32 * pSrc,
369 arm_matrix_instance_f32 * pDst)
370 {
371
372 arm_status status; /* status of matrix inverse */
373
374
375 #ifdef ARM_MATH_MATRIX_CHECK
376
377 /* Check for matrix mismatch condition */
378 if ((pSrc->numRows != pSrc->numCols) ||
379 (pDst->numRows != pDst->numCols) ||
380 (pSrc->numRows != pDst->numRows) )
381 {
382 /* Set status as ARM_MATH_SIZE_MISMATCH */
383 status = ARM_MATH_SIZE_MISMATCH;
384 }
385 else
386
387 #endif /* #ifdef ARM_MATH_MATRIX_CHECK */
388
389 {
390 int i,j,k;
391 int n = pSrc->numRows;
392 float32_t invSqrtVj;
393 float32_t *pA,*pG;
394
395 pA = pSrc->pData;
396 pG = pDst->pData;
397
398
399 for(i=0 ; i < n ; i++)
400 {
401 for(j=i ; j < n ; j++)
402 {
403 pG[j * n + i] = pA[j * n + i];
404
405 for(k=0; k < i ; k++)
406 {
407 pG[j * n + i] = pG[j * n + i] - pG[i * n + k] * pG[j * n + k];
408 }
409 }
410
411 if (pG[i * n + i] <= 0.0f)
412 {
413 return(ARM_MATH_DECOMPOSITION_FAILURE);
414 }
415
416 invSqrtVj = 1.0f/sqrtf(pG[i * n + i]);
417 for(j=i ; j < n ; j++)
418 {
419 pG[j * n + i] = pG[j * n + i] * invSqrtVj ;
420 }
421 }
422
423 status = ARM_MATH_SUCCESS;
424
425 }
426
427
428 /* Return to application */
429 return (status);
430 }
431 #endif /* #if defined(ARM_MATH_NEON) */
432 #endif /* defined(ARM_MATH_MVEF) && !defined(ARM_MATH_AUTOVECTORIZE) */
433
434 /**
435 @} end of MatrixChol group
436 */
437
