1 /* ----------------------------------------------------------------------
2 * Project: CMSIS DSP Library
3 * Title: arm_mat_inverse_f32.c
4 * Description: Floating-point matrix inverse
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 #include "dsp/matrix_utils.h"
31
32
33 /**
34 @ingroup groupMatrix
35 */
36
37 /**
38 @defgroup MatrixInv Matrix Inverse
39
40 Computes the inverse of a matrix.
41
42 The inverse is defined only if the input matrix is square and non-singular (the determinant is non-zero).
43 The function checks that the input and output matrices are square and of the same size.
44
45 Matrix inversion is numerically sensitive and the CMSIS DSP library only supports matrix
46 inversion of floating-point matrices.
47
48 @par Algorithm
49 The Gauss-Jordan method is used to find the inverse.
50 The algorithm performs a sequence of elementary row-operations until it
51 reduces the input matrix to an identity matrix. Applying the same sequence
52 of elementary row-operations to an identity matrix yields the inverse matrix.
53 If the input matrix is singular, then the algorithm terminates and returns error status
54 <code>ARM_MATH_SINGULAR</code>.
55
56 @par Matrix Inverse of a 3 x 3 matrix using Gauss-Jordan Method
57
58 \f[
59 \begin{pmatrix}
60 a_{1,1} & a_{1,2} & a_{1,3} & | & 1 & 0 & 0\\
61 a_{2,1} & a_{2,2} & a_{2,3} & | & 0 & 1 & 0\\
62 a_{3,1} & a_{3,2} & a_{3,3} & | & 0 & 0 & 1\\
63 \end{pmatrix}
64 \rightarrow
65 \begin{pmatrix}
66 1 & 0 & 0 & | & x_{1,1} & x_{2,1} & x_{3,1} \\
67 0 & 1 & 0 & | & x_{1,2} & x_{2,2} & x_{3,2} \\
68 0 & 0 & 1 & | & x_{1,3} & x_{2,3} & x_{3,3} \\
69 \end{pmatrix}
70 \f]
71 */
72
73 /**
74 @addtogroup MatrixInv
75 @{
76 */
77
78 /**
79 @brief Floating-point matrix inverse.
80 @param[in] pSrc points to input matrix structure. The source matrix is modified by the function.
81 @param[out] pDst points to output matrix structure
82 @return execution status
83 - \ref ARM_MATH_SUCCESS : Operation successful
84 - \ref ARM_MATH_SIZE_MISMATCH : Matrix size check failed
85 - \ref ARM_MATH_SINGULAR : Input matrix is found to be singular (non-invertible)
86 */
arm_mat_inverse_f32(const arm_matrix_instance_f32 * pSrc,arm_matrix_instance_f32 * pDst)87 ARM_DSP_ATTRIBUTE arm_status arm_mat_inverse_f32(
88 const arm_matrix_instance_f32 * pSrc,
89 arm_matrix_instance_f32 * pDst)
90 {
91 float32_t *pIn = pSrc->pData; /* input data matrix pointer */
92 float32_t *pOut = pDst->pData; /* output data matrix pointer */
93
94 float32_t *pTmp;
95 uint32_t numRows = pSrc->numRows; /* Number of rows in the matrix */
96 uint32_t numCols = pSrc->numCols; /* Number of Cols in the matrix */
97
98
99 float32_t pivot = 0.0f, newPivot=0.0f; /* Temporary input values */
100 uint32_t selectedRow,pivotRow,i, rowNb, rowCnt, flag = 0U, j,column; /* loop counters */
101 arm_status status; /* status of matrix inverse */
102
103 #ifdef ARM_MATH_MATRIX_CHECK
104
105 /* Check for matrix mismatch condition */
106 if ((pSrc->numRows != pSrc->numCols) ||
107 (pDst->numRows != pDst->numCols) ||
108 (pSrc->numRows != pDst->numRows) )
109 {
110 /* Set status as ARM_MATH_SIZE_MISMATCH */
111 status = ARM_MATH_SIZE_MISMATCH;
112 }
113 else
114
115 #endif /* #ifdef ARM_MATH_MATRIX_CHECK */
116
117 {
118 /*--------------------------------------------------------------------------------------------------------------
119 * Matrix Inverse can be solved using elementary row operations.
120 *
121 * Gauss-Jordan Method:
122 *
123 * 1. First combine the identity matrix and the input matrix separated by a bar to form an
124 * augmented matrix as follows:
125 * _ _ _ _
126 * | a11 a12 | 1 0 | | X11 X12 |
127 * | | | = | |
128 * |_ a21 a22 | 0 1 _| |_ X21 X21 _|
129 *
130 * 2. In our implementation, pDst Matrix is used as identity matrix.
131 *
132 * 3. Begin with the first row. Let i = 1.
133 *
134 * 4. Check to see if the pivot for row i is zero.
135 * The pivot is the element of the main diagonal that is on the current row.
136 * For instance, if working with row i, then the pivot element is aii.
137 * If the pivot is zero, exchange that row with a row below it that does not
138 * contain a zero in column i. If this is not possible, then an inverse
139 * to that matrix does not exist.
140 *
141 * 5. Divide every element of row i by the pivot.
142 *
143 * 6. For every row below and row i, replace that row with the sum of that row and
144 * a multiple of row i so that each new element in column i below row i is zero.
145 *
146 * 7. Move to the next row and column and repeat steps 2 through 5 until you have zeros
147 * for every element below and above the main diagonal.
148 *
149 * 8. Now an identical matrix is formed to the left of the bar(input matrix, pSrc).
150 * Therefore, the matrix to the right of the bar is our solution(pDst matrix, pDst).
151 *----------------------------------------------------------------------------------------------------------------*/
152
153 /* Working pointer for destination matrix */
154 pTmp = pOut;
155
156 /* Loop over the number of rows */
157 rowCnt = numRows;
158
159 /* Making the destination matrix as identity matrix */
160 while (rowCnt > 0U)
161 {
162 /* Writing all zeroes in lower triangle of the destination matrix */
163 j = numRows - rowCnt;
164 while (j > 0U)
165 {
166 *pTmp++ = 0.0f;
167 j--;
168 }
169
170 /* Writing all ones in the diagonal of the destination matrix */
171 *pTmp++ = 1.0f;
172
173 /* Writing all zeroes in upper triangle of the destination matrix */
174 j = rowCnt - 1U;
175 while (j > 0U)
176 {
177 *pTmp++ = 0.0f;
178 j--;
179 }
180
181 /* Decrement loop counter */
182 rowCnt--;
183 }
184
185 /* Loop over the number of columns of the input matrix.
186 All the elements in each column are processed by the row operations */
187
188 /* Index modifier to navigate through the columns */
189 for(column = 0U; column < numCols; column++)
190 {
191 /* reset flag */
192 flag = 0;
193
194 /* Check if the pivot element is zero..
195 * If it is zero then interchange the row with non zero row below.
196 * If there is no non zero element to replace in the rows below,
197 * then the matrix is Singular. */
198
199 pivotRow = column;
200
201 /* Temporary variable to hold the pivot value */
202 pTmp = ELEM(pSrc,column,column) ;
203 pivot = *pTmp;
204 selectedRow = column;
205
206 /* Find maximum pivot in column */
207
208 /* Loop over the number rows present below */
209
210 for (rowNb = column+1; rowNb < numRows; rowNb++)
211 {
212 /* Update the input and destination pointers */
213 pTmp = ELEM(pSrc,rowNb,column);
214 newPivot = *pTmp;
215 if (fabsf(newPivot) > fabsf(pivot))
216 {
217 selectedRow = rowNb;
218 pivot = newPivot;
219 }
220 }
221
222 /* Check if there is a non zero pivot element to
223 * replace in the rows below */
224 if ((pivot != 0.0f) && (selectedRow != column))
225 {
226
227 SWAP_ROWS_F32(pSrc,column, pivotRow,selectedRow);
228 SWAP_ROWS_F32(pDst,0, pivotRow,selectedRow);
229
230
231 /* Flag to indicate whether exchange is done or not */
232 flag = 1U;
233 }
234
235
236
237
238
239 /* Update the status if the matrix is singular */
240 if ((flag != 1U) && (pivot == 0.0f))
241 {
242 return ARM_MATH_SINGULAR;
243 }
244
245
246 /* Pivot element of the row */
247 pivot = 1.0f / pivot;
248
249 SCALE_ROW_F32(pSrc,column,pivot,pivotRow);
250 SCALE_ROW_F32(pDst,0,pivot,pivotRow);
251
252
253 /* Replace the rows with the sum of that row and a multiple of row i
254 * so that each new element in column i above row i is zero.*/
255
256 rowNb = 0;
257 for (;rowNb < pivotRow; rowNb++)
258 {
259 pTmp = ELEM(pSrc,rowNb,column) ;
260 pivot = *pTmp;
261
262 MAS_ROW_F32(column,pSrc,rowNb,pivot,pSrc,pivotRow);
263 MAS_ROW_F32(0 ,pDst,rowNb,pivot,pDst,pivotRow);
264
265
266 }
267
268 for (rowNb = pivotRow + 1; rowNb < numRows; rowNb++)
269 {
270 pTmp = ELEM(pSrc,rowNb,column) ;
271 pivot = *pTmp;
272
273 MAS_ROW_F32(column,pSrc,rowNb,pivot,pSrc,pivotRow);
274 MAS_ROW_F32(0 ,pDst,rowNb,pivot,pDst,pivotRow);
275
276 }
277
278 }
279
280 /* Set status as ARM_MATH_SUCCESS */
281 status = ARM_MATH_SUCCESS;
282
283 if ((flag != 1U) && (pivot == 0.0f))
284 {
285 pIn = pSrc->pData;
286 for (i = 0; i < numRows * numCols; i++)
287 {
288 if (pIn[i] != 0.0f)
289 break;
290 }
291
292 if (i == numRows * numCols)
293 status = ARM_MATH_SINGULAR;
294 }
295 }
296
297 /* Return to application */
298 return (status);
299 }
300 /**
301 @} end of MatrixInv group
302 */
303
