1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_householder_f32.c
4 * Description: 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 @defgroup MatrixHouseholder Householder transform of a vector
45
46 Computes the Householder transform of a vector x.
47
48 The Householder transform of x is a vector v with
49
50 \f[
51 v_0 = 1
52 \f]
53
54 and a scalar \f$\beta\f$ such that:
55
56 \f[
57 P = I - \beta v v^T
58 \f]
59
60 is an orthogonal matrix and
61
62 \f[
63 P x = ||x||_2 e_1
64 \f]
65
66 So P is an hyperplane reflection such that the image of x
67 is proportional to \f$e_1\f$.
68
69 \f$e_1\f$ is the vector of coordinates:
70
71 \f[
72 \begin{pmatrix}
73 1 \\
74 0 \\
75 \vdots \\
76 \end{pmatrix}
77 \f]
78
79 If x is already proportional to \f$e_1\f$ then
80 the matrix P should be the identity.
81
82 Thus, \f$\beta\f$ should be 0 and in this case the vector v
83 can also be null.
84
85 But how do we detect that x is already proportional to
86 \f$e_1\f$.
87
88 If x
89 \f[
90 x =
91 \begin{pmatrix}
92 x_0 \\
93 xr \\
94 \end{pmatrix}
95 \f]
96
97 where \f$xr\f$ is a vector.
98
99 The algorithm is computing the norm squared of this vector:
100
101 \f[
102 ||xr||^2
103 \f]
104
105 and this value is compared to a `threshold`. If the value
106 is smaller than the `threshold`, the algorithm is
107 returning 0 for \f$\beta\f$ and the householder vector.
108
109 This `threshold` is an argument of the function.
110
111 Default values are provided in the header
112 `dsp/matrix_functions.h` like for instance
113 `DEFAULT_HOUSEHOLDER_THRESHOLD_F32`
114
115
116
117 */
118
119 /**
120 @addtogroup MatrixHouseholder
121 @{
122 */
123
124 /**
125 @brief Householder transform of a floating point vector.
126 @param[in] pSrc points to the input vector.
127 @param[in] threshold norm2 threshold.
128 @param[in] blockSize dimension of the vector space.
129 @param[out] pOut points to the output vector.
130 @return beta return the scaling factor beta
131 */
132
133
134
135
arm_householder_f32(const float32_t * pSrc,const float32_t threshold,uint32_t blockSize,float32_t * pOut)136 float32_t arm_householder_f32(
137 const float32_t * pSrc,
138 const float32_t threshold,
139 uint32_t blockSize,
140 float32_t * pOut
141 )
142
143 {
144 uint32_t i;
145 float32_t epsilon;
146 float32_t x1norm2,alpha;
147 float32_t beta,tau,r;
148
149 epsilon = threshold;
150
151 alpha = pSrc[0];
152
153 for(i=1; i < blockSize; i++)
154 {
155 pOut[i] = pSrc[i];
156 }
157 pOut[0] = 1.0f;
158
159 arm_dot_prod_f32(pSrc+1,pSrc+1,blockSize-1,&x1norm2);
160
161 if (x1norm2<=epsilon)
162 {
163 tau = 0.0f;
164 memset(pOut,0,blockSize * sizeof(float32_t));
165 }
166 else
167 {
168 beta = alpha * alpha + x1norm2;
169 (void)arm_sqrt_f32(beta,&beta);
170
171 if (alpha > 0.0f)
172 {
173 beta = -beta;
174 }
175
176 r = 1.0f / (alpha -beta);
177 arm_scale_f32(pOut,r,pOut,blockSize);
178 pOut[0] = 1.0f;
179
180
181 tau = (beta - alpha) / beta;
182
183 }
184
185 return(tau);
186
187 }
188
189
190 /**
191 @} end of MatrixHouseholder group
192 */
193
