QuakeGod
2023-02-01 6126f6a78b14297cefb02f06ba58806767d424b5
提交 | 用户 | age
bfc108 1 /* ----------------------------------------------------------------------    
Q 2 * Copyright (C) 2010-2014 ARM Limited. All rights reserved.    
3 *    
4 * $Date:        19. March 2015
5 * $Revision:     V.1.4.5
6 *    
7 * Project:         CMSIS DSP Library    
8 * Title:        arm_mat_inverse_f32.c    
9 *    
10 * Description:    Floating-point matrix inverse.    
11 *    
12 * Target Processor: Cortex-M4/Cortex-M3/Cortex-M0
13 *  
14 * Redistribution and use in source and binary forms, with or without 
15 * modification, are permitted provided that the following conditions
16 * are met:
17 *   - Redistributions of source code must retain the above copyright
18 *     notice, this list of conditions and the following disclaimer.
19 *   - Redistributions in binary form must reproduce the above copyright
20 *     notice, this list of conditions and the following disclaimer in
21 *     the documentation and/or other materials provided with the 
22 *     distribution.
23 *   - Neither the name of ARM LIMITED nor the names of its contributors
24 *     may be used to endorse or promote products derived from this
25 *     software without specific prior written permission.
26 *
27 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
28 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
29 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
30 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 
31 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
32 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
33 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
34 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
35 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
36 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
37 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
38 * POSSIBILITY OF SUCH DAMAGE.    
39 * -------------------------------------------------------------------- */
40
41 #include "arm_math.h"
42
43 /**    
44  * @ingroup groupMatrix    
45  */
46
47 /**    
48  * @defgroup MatrixInv Matrix Inverse    
49  *    
50  * Computes the inverse of a matrix.    
51  *    
52  * The inverse is defined only if the input matrix is square and non-singular (the determinant    
53  * is non-zero). The function checks that the input and output matrices are square and of the    
54  * same size.    
55  *    
56  * Matrix inversion is numerically sensitive and the CMSIS DSP library only supports matrix    
57  * inversion of floating-point matrices.    
58  *    
59  * \par Algorithm    
60  * The Gauss-Jordan method is used to find the inverse.    
61  * The algorithm performs a sequence of elementary row-operations until it    
62  * reduces the input matrix to an identity matrix. Applying the same sequence    
63  * of elementary row-operations to an identity matrix yields the inverse matrix.    
64  * If the input matrix is singular, then the algorithm terminates and returns error status    
65  * <code>ARM_MATH_SINGULAR</code>.    
66  * \image html MatrixInverse.gif "Matrix Inverse of a 3 x 3 matrix using Gauss-Jordan Method"    
67  */
68
69 /**    
70  * @addtogroup MatrixInv    
71  * @{    
72  */
73
74 /**    
75  * @brief Floating-point matrix inverse.    
76  * @param[in]       *pSrc points to input matrix structure    
77  * @param[out]      *pDst points to output matrix structure    
78  * @return             The function returns    
79  * <code>ARM_MATH_SIZE_MISMATCH</code> if the input matrix is not square or if the size    
80  * of the output matrix does not match the size of the input matrix.    
81  * If the input matrix is found to be singular (non-invertible), then the function returns    
82  * <code>ARM_MATH_SINGULAR</code>.  Otherwise, the function returns <code>ARM_MATH_SUCCESS</code>.    
83  */
84
85 arm_status arm_mat_inverse_f32(
86   const arm_matrix_instance_f32 * pSrc,
87   arm_matrix_instance_f32 * pDst)
88 {
89   float32_t *pIn = pSrc->pData;                  /* input data matrix pointer */
90   float32_t *pOut = pDst->pData;                 /* output data matrix pointer */
91   float32_t *pInT1, *pInT2;                      /* Temporary input data matrix pointer */
92   float32_t *pOutT1, *pOutT2;                    /* Temporary output data matrix pointer */
93   float32_t *pPivotRowIn, *pPRT_in, *pPivotRowDst, *pPRT_pDst;  /* Temporary input and output data matrix pointer */
94   uint32_t numRows = pSrc->numRows;              /* Number of rows in the matrix  */
95   uint32_t numCols = pSrc->numCols;              /* Number of Cols in the matrix  */
96
97 #ifndef ARM_MATH_CM0_FAMILY
98   float32_t maxC;                                /* maximum value in the column */
99
100   /* Run the below code for Cortex-M4 and Cortex-M3 */
101
102   float32_t Xchg, in = 0.0f, in1;                /* Temporary input values  */
103   uint32_t i, rowCnt, flag = 0u, j, loopCnt, k, l;      /* loop counters */
104   arm_status status;                             /* status of matrix inverse */
105
106 #ifdef ARM_MATH_MATRIX_CHECK
107
108
109   /* Check for matrix mismatch condition */
110   if((pSrc->numRows != pSrc->numCols) || (pDst->numRows != pDst->numCols)
111      || (pSrc->numRows != pDst->numRows))
112   {
113     /* Set status as ARM_MATH_SIZE_MISMATCH */
114     status = ARM_MATH_SIZE_MISMATCH;
115   }
116   else
117 #endif /*    #ifdef ARM_MATH_MATRIX_CHECK    */
118
119   {
120
121     /*--------------------------------------------------------------------------------------------------------------    
122      * Matrix Inverse can be solved using elementary row operations.    
123      *    
124      *    Gauss-Jordan Method:    
125      *    
126      *       1. First combine the identity matrix and the input matrix separated by a bar to form an    
127      *        augmented matrix as follows:    
128      *                        _                      _         _           _    
129      *                       |  a11  a12 | 1   0  |       |  X11 X12  |    
130      *                       |           |        |   =   |           |    
131      *                       |_ a21  a22 | 0   1 _|       |_ X21 X21 _|    
132      *    
133      *        2. In our implementation, pDst Matrix is used as identity matrix.    
134      *    
135      *        3. Begin with the first row. Let i = 1.    
136      *    
137      *        4. Check to see if the pivot for column i is the greatest of the column.    
138      *           The pivot is the element of the main diagonal that is on the current row.    
139      *           For instance, if working with row i, then the pivot element is aii.    
140      *           If the pivot is not the most significant of the columns, exchange that row with a row
141      *           below it that does contain the most significant value in column i. If the most
142      *         significant value of the column is zero, then an inverse to that matrix does not exist.
143      *           The most significant value of the column is the absolute maximum.
144      *    
145      *        5. Divide every element of row i by the pivot.    
146      *    
147      *        6. For every row below and  row i, replace that row with the sum of that row and    
148      *           a multiple of row i so that each new element in column i below row i is zero.    
149      *    
150      *        7. Move to the next row and column and repeat steps 2 through 5 until you have zeros    
151      *           for every element below and above the main diagonal.    
152      *    
153      *        8. Now an identical matrix is formed to the left of the bar(input matrix, pSrc).    
154      *           Therefore, the matrix to the right of the bar is our solution(pDst matrix, pDst).    
155      *----------------------------------------------------------------------------------------------------------------*/
156
157     /* Working pointer for destination matrix */
158     pOutT1 = pOut;
159
160     /* Loop over the number of rows */
161     rowCnt = numRows;
162
163     /* Making the destination matrix as identity matrix */
164     while(rowCnt > 0u)
165     {
166       /* Writing all zeroes in lower triangle of the destination matrix */
167       j = numRows - rowCnt;
168       while(j > 0u)
169       {
170         *pOutT1++ = 0.0f;
171         j--;
172       }
173
174       /* Writing all ones in the diagonal of the destination matrix */
175       *pOutT1++ = 1.0f;
176
177       /* Writing all zeroes in upper triangle of the destination matrix */
178       j = rowCnt - 1u;
179       while(j > 0u)
180       {
181         *pOutT1++ = 0.0f;
182         j--;
183       }
184
185       /* Decrement the loop counter */
186       rowCnt--;
187     }
188
189     /* Loop over the number of columns of the input matrix.    
190        All the elements in each column are processed by the row operations */
191     loopCnt = numCols;
192
193     /* Index modifier to navigate through the columns */
194     l = 0u;
195
196     while(loopCnt > 0u)
197     {
198       /* Check if the pivot element is zero..    
199        * If it is zero then interchange the row with non zero row below.    
200        * If there is no non zero element to replace in the rows below,    
201        * then the matrix is Singular. */
202
203       /* Working pointer for the input matrix that points    
204        * to the pivot element of the particular row  */
205       pInT1 = pIn + (l * numCols);
206
207       /* Working pointer for the destination matrix that points    
208        * to the pivot element of the particular row  */
209       pOutT1 = pOut + (l * numCols);
210
211       /* Temporary variable to hold the pivot value */
212       in = *pInT1;
213
214       /* Grab the most significant value from column l */
215       maxC = 0;
216       for (i = l; i < numRows; i++)
217       {
218         maxC = *pInT1 > 0 ? (*pInT1 > maxC ? *pInT1 : maxC) : (-*pInT1 > maxC ? -*pInT1 : maxC);
219         pInT1 += numCols;
220       }
221
222       /* Update the status if the matrix is singular */
223       if(maxC == 0.0f)
224       {
225         return ARM_MATH_SINGULAR;
226       }
227
228       /* Restore pInT1  */
229       pInT1 = pIn;
230
231       /* Destination pointer modifier */
232       k = 1u;
233       
234       /* Check if the pivot element is the most significant of the column */
235       if( (in > 0.0f ? in : -in) != maxC)
236       {
237         /* Loop over the number rows present below */
238         i = numRows - (l + 1u);
239
240         while(i > 0u)
241         {
242           /* Update the input and destination pointers */
243           pInT2 = pInT1 + (numCols * l);
244           pOutT2 = pOutT1 + (numCols * k);
245
246           /* Look for the most significant element to    
247            * replace in the rows below */
248           if((*pInT2 > 0.0f ? *pInT2: -*pInT2) == maxC)
249           {
250             /* Loop over number of columns    
251              * to the right of the pilot element */
252             j = numCols - l;
253
254             while(j > 0u)
255             {
256               /* Exchange the row elements of the input matrix */
257               Xchg = *pInT2;
258               *pInT2++ = *pInT1;
259               *pInT1++ = Xchg;
260
261               /* Decrement the loop counter */
262               j--;
263             }
264
265             /* Loop over number of columns of the destination matrix */
266             j = numCols;
267
268             while(j > 0u)
269             {
270               /* Exchange the row elements of the destination matrix */
271               Xchg = *pOutT2;
272               *pOutT2++ = *pOutT1;
273               *pOutT1++ = Xchg;
274
275               /* Decrement the loop counter */
276               j--;
277             }
278
279             /* Flag to indicate whether exchange is done or not */
280             flag = 1u;
281
282             /* Break after exchange is done */
283             break;
284           }
285
286           /* Update the destination pointer modifier */
287           k++;
288
289           /* Decrement the loop counter */
290           i--;
291         }
292       }
293
294       /* Update the status if the matrix is singular */
295       if((flag != 1u) && (in == 0.0f))
296       {
297         return ARM_MATH_SINGULAR;
298       }
299
300       /* Points to the pivot row of input and destination matrices */
301       pPivotRowIn = pIn + (l * numCols);
302       pPivotRowDst = pOut + (l * numCols);
303
304       /* Temporary pointers to the pivot row pointers */
305       pInT1 = pPivotRowIn;
306       pInT2 = pPivotRowDst;
307
308       /* Pivot element of the row */
309       in = *pPivotRowIn;
310
311       /* Loop over number of columns    
312        * to the right of the pilot element */
313       j = (numCols - l);
314
315       while(j > 0u)
316       {
317         /* Divide each element of the row of the input matrix    
318          * by the pivot element */
319         in1 = *pInT1;
320         *pInT1++ = in1 / in;
321
322         /* Decrement the loop counter */
323         j--;
324       }
325
326       /* Loop over number of columns of the destination matrix */
327       j = numCols;
328
329       while(j > 0u)
330       {
331         /* Divide each element of the row of the destination matrix    
332          * by the pivot element */
333         in1 = *pInT2;
334         *pInT2++ = in1 / in;
335
336         /* Decrement the loop counter */
337         j--;
338       }
339
340       /* Replace the rows with the sum of that row and a multiple of row i    
341        * so that each new element in column i above row i is zero.*/
342
343       /* Temporary pointers for input and destination matrices */
344       pInT1 = pIn;
345       pInT2 = pOut;
346
347       /* index used to check for pivot element */
348       i = 0u;
349
350       /* Loop over number of rows */
351       /*  to be replaced by the sum of that row and a multiple of row i */
352       k = numRows;
353
354       while(k > 0u)
355       {
356         /* Check for the pivot element */
357         if(i == l)
358         {
359           /* If the processing element is the pivot element,    
360              only the columns to the right are to be processed */
361           pInT1 += numCols - l;
362
363           pInT2 += numCols;
364         }
365         else
366         {
367           /* Element of the reference row */
368           in = *pInT1;
369
370           /* Working pointers for input and destination pivot rows */
371           pPRT_in = pPivotRowIn;
372           pPRT_pDst = pPivotRowDst;
373
374           /* Loop over the number of columns to the right of the pivot element,    
375              to replace the elements in the input matrix */
376           j = (numCols - l);
377
378           while(j > 0u)
379           {
380             /* Replace the element by the sum of that row    
381                and a multiple of the reference row  */
382             in1 = *pInT1;
383             *pInT1++ = in1 - (in * *pPRT_in++);
384
385             /* Decrement the loop counter */
386             j--;
387           }
388
389           /* Loop over the number of columns to    
390              replace the elements in the destination matrix */
391           j = numCols;
392
393           while(j > 0u)
394           {
395             /* Replace the element by the sum of that row    
396                and a multiple of the reference row  */
397             in1 = *pInT2;
398             *pInT2++ = in1 - (in * *pPRT_pDst++);
399
400             /* Decrement the loop counter */
401             j--;
402           }
403
404         }
405
406         /* Increment the temporary input pointer */
407         pInT1 = pInT1 + l;
408
409         /* Decrement the loop counter */
410         k--;
411
412         /* Increment the pivot index */
413         i++;
414       }
415
416       /* Increment the input pointer */
417       pIn++;
418
419       /* Decrement the loop counter */
420       loopCnt--;
421
422       /* Increment the index modifier */
423       l++;
424     }
425
426
427 #else
428
429   /* Run the below code for Cortex-M0 */
430
431   float32_t Xchg, in = 0.0f;                     /* Temporary input values  */
432   uint32_t i, rowCnt, flag = 0u, j, loopCnt, k, l;      /* loop counters */
433   arm_status status;                             /* status of matrix inverse */
434
435 #ifdef ARM_MATH_MATRIX_CHECK
436
437   /* Check for matrix mismatch condition */
438   if((pSrc->numRows != pSrc->numCols) || (pDst->numRows != pDst->numCols)
439      || (pSrc->numRows != pDst->numRows))
440   {
441     /* Set status as ARM_MATH_SIZE_MISMATCH */
442     status = ARM_MATH_SIZE_MISMATCH;
443   }
444   else
445 #endif /*      #ifdef ARM_MATH_MATRIX_CHECK    */
446   {
447
448     /*--------------------------------------------------------------------------------------------------------------       
449      * Matrix Inverse can be solved using elementary row operations.        
450      *        
451      *    Gauss-Jordan Method:       
452      *                
453      *       1. First combine the identity matrix and the input matrix separated by a bar to form an        
454      *        augmented matrix as follows:        
455      *                        _  _          _        _       _   _         _           _       
456      *                       |  |  a11  a12  | | | 1   0  |   |       |  X11 X12  |         
457      *                       |  |            | | |        |   |   =   |           |        
458      *                       |_ |_ a21  a22 _| | |_0   1 _|  _|       |_ X21 X21 _|       
459      *                              
460      *        2. In our implementation, pDst Matrix is used as identity matrix.    
461      *       
462      *        3. Begin with the first row. Let i = 1.       
463      *       
464      *        4. Check to see if the pivot for row i is zero.       
465      *           The pivot is the element of the main diagonal that is on the current row.       
466      *           For instance, if working with row i, then the pivot element is aii.       
467      *           If the pivot is zero, exchange that row with a row below it that does not        
468      *           contain a zero in column i. If this is not possible, then an inverse        
469      *           to that matrix does not exist.       
470      *           
471      *        5. Divide every element of row i by the pivot.       
472      *           
473      *        6. For every row below and  row i, replace that row with the sum of that row and        
474      *           a multiple of row i so that each new element in column i below row i is zero.       
475      *           
476      *        7. Move to the next row and column and repeat steps 2 through 5 until you have zeros       
477      *           for every element below and above the main diagonal.        
478      *                             
479      *        8. Now an identical matrix is formed to the left of the bar(input matrix, src).       
480      *           Therefore, the matrix to the right of the bar is our solution(dst matrix, dst).         
481      *----------------------------------------------------------------------------------------------------------------*/
482
483     /* Working pointer for destination matrix */
484     pOutT1 = pOut;
485
486     /* Loop over the number of rows */
487     rowCnt = numRows;
488
489     /* Making the destination matrix as identity matrix */
490     while(rowCnt > 0u)
491     {
492       /* Writing all zeroes in lower triangle of the destination matrix */
493       j = numRows - rowCnt;
494       while(j > 0u)
495       {
496         *pOutT1++ = 0.0f;
497         j--;
498       }
499
500       /* Writing all ones in the diagonal of the destination matrix */
501       *pOutT1++ = 1.0f;
502
503       /* Writing all zeroes in upper triangle of the destination matrix */
504       j = rowCnt - 1u;
505       while(j > 0u)
506       {
507         *pOutT1++ = 0.0f;
508         j--;
509       }
510
511       /* Decrement the loop counter */
512       rowCnt--;
513     }
514
515     /* Loop over the number of columns of the input matrix.     
516        All the elements in each column are processed by the row operations */
517     loopCnt = numCols;
518
519     /* Index modifier to navigate through the columns */
520     l = 0u;
521     //for(loopCnt = 0u; loopCnt < numCols; loopCnt++)   
522     while(loopCnt > 0u)
523     {
524       /* Check if the pivot element is zero..    
525        * If it is zero then interchange the row with non zero row below.   
526        * If there is no non zero element to replace in the rows below,   
527        * then the matrix is Singular. */
528
529       /* Working pointer for the input matrix that points     
530        * to the pivot element of the particular row  */
531       pInT1 = pIn + (l * numCols);
532
533       /* Working pointer for the destination matrix that points     
534        * to the pivot element of the particular row  */
535       pOutT1 = pOut + (l * numCols);
536
537       /* Temporary variable to hold the pivot value */
538       in = *pInT1;
539
540       /* Destination pointer modifier */
541       k = 1u;
542
543       /* Check if the pivot element is zero */
544       if(*pInT1 == 0.0f)
545       {
546         /* Loop over the number rows present below */
547         for (i = (l + 1u); i < numRows; i++)
548         {
549           /* Update the input and destination pointers */
550           pInT2 = pInT1 + (numCols * l);
551           pOutT2 = pOutT1 + (numCols * k);
552
553           /* Check if there is a non zero pivot element to     
554            * replace in the rows below */
555           if(*pInT2 != 0.0f)
556           {
557             /* Loop over number of columns     
558              * to the right of the pilot element */
559             for (j = 0u; j < (numCols - l); j++)
560             {
561               /* Exchange the row elements of the input matrix */
562               Xchg = *pInT2;
563               *pInT2++ = *pInT1;
564               *pInT1++ = Xchg;
565             }
566
567             for (j = 0u; j < numCols; j++)
568             {
569               Xchg = *pOutT2;
570               *pOutT2++ = *pOutT1;
571               *pOutT1++ = Xchg;
572             }
573
574             /* Flag to indicate whether exchange is done or not */
575             flag = 1u;
576
577             /* Break after exchange is done */
578             break;
579           }
580
581           /* Update the destination pointer modifier */
582           k++;
583         }
584       }
585
586       /* Update the status if the matrix is singular */
587       if((flag != 1u) && (in == 0.0f))
588       {
589         return ARM_MATH_SINGULAR;
590       }
591
592       /* Points to the pivot row of input and destination matrices */
593       pPivotRowIn = pIn + (l * numCols);
594       pPivotRowDst = pOut + (l * numCols);
595
596       /* Temporary pointers to the pivot row pointers */
597       pInT1 = pPivotRowIn;
598       pOutT1 = pPivotRowDst;
599
600       /* Pivot element of the row */
601       in = *(pIn + (l * numCols));
602
603       /* Loop over number of columns     
604        * to the right of the pilot element */
605       for (j = 0u; j < (numCols - l); j++)
606       {
607         /* Divide each element of the row of the input matrix     
608          * by the pivot element */
609         *pInT1 = *pInT1 / in;
610         pInT1++;
611       }
612       for (j = 0u; j < numCols; j++)
613       {
614         /* Divide each element of the row of the destination matrix     
615          * by the pivot element */
616         *pOutT1 = *pOutT1 / in;
617         pOutT1++;
618       }
619
620       /* Replace the rows with the sum of that row and a multiple of row i     
621        * so that each new element in column i above row i is zero.*/
622
623       /* Temporary pointers for input and destination matrices */
624       pInT1 = pIn;
625       pOutT1 = pOut;
626
627       for (i = 0u; i < numRows; i++)
628       {
629         /* Check for the pivot element */
630         if(i == l)
631         {
632           /* If the processing element is the pivot element,     
633              only the columns to the right are to be processed */
634           pInT1 += numCols - l;
635           pOutT1 += numCols;
636         }
637         else
638         {
639           /* Element of the reference row */
640           in = *pInT1;
641
642           /* Working pointers for input and destination pivot rows */
643           pPRT_in = pPivotRowIn;
644           pPRT_pDst = pPivotRowDst;
645
646           /* Loop over the number of columns to the right of the pivot element,     
647              to replace the elements in the input matrix */
648           for (j = 0u; j < (numCols - l); j++)
649           {
650             /* Replace the element by the sum of that row     
651                and a multiple of the reference row  */
652             *pInT1 = *pInT1 - (in * *pPRT_in++);
653             pInT1++;
654           }
655           /* Loop over the number of columns to     
656              replace the elements in the destination matrix */
657           for (j = 0u; j < numCols; j++)
658           {
659             /* Replace the element by the sum of that row     
660                and a multiple of the reference row  */
661             *pOutT1 = *pOutT1 - (in * *pPRT_pDst++);
662             pOutT1++;
663           }
664
665         }
666         /* Increment the temporary input pointer */
667         pInT1 = pInT1 + l;
668       }
669       /* Increment the input pointer */
670       pIn++;
671
672       /* Decrement the loop counter */
673       loopCnt--;
674       /* Increment the index modifier */
675       l++;
676     }
677
678
679 #endif /* #ifndef ARM_MATH_CM0_FAMILY */
680
681     /* Set status as ARM_MATH_SUCCESS */
682     status = ARM_MATH_SUCCESS;
683
684     if((flag != 1u) && (in == 0.0f))
685     {
686       pIn = pSrc->pData;
687       for (i = 0; i < numRows * numCols; i++)
688       {
689         if (pIn[i] != 0.0f)
690             break;
691       }
692       
693       if (i == numRows * numCols)
694         status = ARM_MATH_SINGULAR;
695     }
696   }
697   /* Return to application */
698   return (status);
699 }
700
701 /**    
702  * @} end of MatrixInv group    
703  */