1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_householder_f64.c
4 * Description: Double floating-point Householder transform
5 *
6 * $Date: 15 June 2022
7 * $Revision: V1.11.0
8 *
9 * Target Processor: Cortex-M and Cortex-A cores
10 * -------------------------------------------------------------------- */
11 /*
12 * Copyright (C) 2010-2022 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 #include "dsp/basic_math_functions.h"
31 #include "dsp/fast_math_functions.h"
32
33 #include "dsp/matrix_utils.h"
34
35 #include <math.h>
36
37
38
39 /**
40 @ingroup groupMatrix
41 */
42
43
44 /**
45 @addtogroup MatrixHouseholder
46 @{
47 */
48
49 /**
50 @brief Householder transform of a double floating point vector.
51 @param[in] pSrc points to the input vector.
52 @param[in] threshold norm2 threshold.
53 @param[in] blockSize dimension of the vector space.
54 @param[out] pOut points to the output vector.
55 @return beta return the scaling factor beta
56 */
57
58
59
60
arm_householder_f64(const float64_t * pSrc,const float64_t threshold,uint32_t blockSize,float64_t * pOut)61 ARM_DSP_ATTRIBUTE float64_t arm_householder_f64(
62 const float64_t * pSrc,
63 const float64_t threshold,
64 uint32_t blockSize,
65 float64_t * pOut
66 )
67
68 {
69 uint32_t i;
70 float64_t epsilon;
71 float64_t x1norm2,alpha;
72 float64_t beta,tau,r;
73
74 epsilon = threshold;
75
76 alpha = pSrc[0];
77
78 for(i=1; i < blockSize; i++)
79 {
80 pOut[i] = pSrc[i];
81 }
82 pOut[0] = 1.0;
83
84 arm_dot_prod_f64(pSrc+1,pSrc+1,blockSize-1,&x1norm2);
85
86 if (x1norm2<=epsilon)
87 {
88 tau = 0.0;
89 memset(pOut,0,blockSize * sizeof(float64_t));
90 }
91 else
92 {
93 beta = alpha * alpha + x1norm2;
94 beta=sqrt(beta);
95
96 if (alpha > 0.0)
97 {
98 beta = -beta;
99 }
100
101 r = 1.0 / (alpha -beta);
102 arm_scale_f64(pOut,r,pOut,blockSize);
103 pOut[0] = 1.0;
104
105
106 tau = (beta - alpha) / beta;
107
108 }
109
110 return(tau);
111
112 }
113
114
115 /**
116 @} end of MatrixHouseholder group
117 */
118
