]> git.donarmstrong.com Git - qmk_firmware.git/blob - tool/mbed/mbed-sdk/libraries/dsp/cmsis_dsp/MatrixFunctions/arm_mat_inverse_f32.c
Squashed 'tmk_core/' changes from 7967731..b9e0ea0
[qmk_firmware.git] / tool / mbed / mbed-sdk / libraries / dsp / cmsis_dsp / MatrixFunctions / arm_mat_inverse_f32.c
1 /* ----------------------------------------------------------------------    
2 * Copyright (C) 2010-2013 ARM Limited. All rights reserved.    
3 *    
4 * $Date:        1. March 2013 
5 * $Revision:    V1.4.1
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 till 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 *pInT3, *pInT4;                      /* 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 coluimns, 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 absolut 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     pInT2 = 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         *pInT2++ = 0.0f;
171         j--;
172       }
173
174       /* Writing all ones in the diagonal of the destination matrix */
175       *pInT2++ = 1.0f;
176
177       /* Writing all zeroes in upper triangle of the destination matrix */
178       j = rowCnt - 1u;
179       while(j > 0u)
180       {
181         *pInT2++ = 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       pInT3 = pOut + (l * numCols);
210
211       /* Temporary variable to hold the pivot value */
212       in = *pInT1;
213
214       /* Destination pointer modifier */
215       k = 1u;
216
217      /* Grab the most significant value from column l */
218       maxC = 0;
219       for (i = 0; i < numRows; i++)
220       {
221         maxC = *pInT1 > 0 ? (*pInT1 > maxC ? *pInT1 : maxC) : (-*pInT1 > maxC ? -*pInT1 : maxC);
222         pInT1 += numCols;
223       }
224
225       /* Update the status if the matrix is singular */
226       if(maxC == 0.0f)
227       {
228         status = ARM_MATH_SINGULAR;
229         break;
230       }
231
232       /* Restore pInT1  */
233       pInT1 -= numRows * numCols;
234       
235       /* Check if the pivot element is the most significant of the column */
236       if( (in > 0.0f ? in : -in) != maxC)
237       {
238         /* Loop over the number rows present below */
239         i = numRows - (l + 1u);
240
241         while(i > 0u)
242         {
243           /* Update the input and destination pointers */
244           pInT2 = pInT1 + (numCols * l);
245           pInT4 = pInT3 + (numCols * k);
246
247           /* Look for the most significant element to    
248            * replace in the rows below */
249           if((*pInT2 > 0.0f ? *pInT2: -*pInT2) == maxC)
250           {
251             /* Loop over number of columns    
252              * to the right of the pilot element */
253             j = numCols - l;
254
255             while(j > 0u)
256             {
257               /* Exchange the row elements of the input matrix */
258               Xchg = *pInT2;
259               *pInT2++ = *pInT1;
260               *pInT1++ = Xchg;
261
262               /* Decrement the loop counter */
263               j--;
264             }
265
266             /* Loop over number of columns of the destination matrix */
267             j = numCols;
268
269             while(j > 0u)
270             {
271               /* Exchange the row elements of the destination matrix */
272               Xchg = *pInT4;
273               *pInT4++ = *pInT3;
274               *pInT3++ = Xchg;
275
276               /* Decrement the loop counter */
277               j--;
278             }
279
280             /* Flag to indicate whether exchange is done or not */
281             flag = 1u;
282
283             /* Break after exchange is done */
284             break;
285           }
286
287           /* Update the destination pointer modifier */
288           k++;
289
290           /* Decrement the loop counter */
291           i--;
292         }
293       }
294
295       /* Update the status if the matrix is singular */
296       if((flag != 1u) && (in == 0.0f))
297       {
298         status = ARM_MATH_SINGULAR;
299
300         break;
301       }
302
303       /* Points to the pivot row of input and destination matrices */
304       pPivotRowIn = pIn + (l * numCols);
305       pPivotRowDst = pOut + (l * numCols);
306
307       /* Temporary pointers to the pivot row pointers */
308       pInT1 = pPivotRowIn;
309       pInT2 = pPivotRowDst;
310
311       /* Pivot element of the row */
312       in = *pPivotRowIn;
313
314       /* Loop over number of columns    
315        * to the right of the pilot element */
316       j = (numCols - l);
317
318       while(j > 0u)
319       {
320         /* Divide each element of the row of the input matrix    
321          * by the pivot element */
322         in1 = *pInT1;
323         *pInT1++ = in1 / in;
324
325         /* Decrement the loop counter */
326         j--;
327       }
328
329       /* Loop over number of columns of the destination matrix */
330       j = numCols;
331
332       while(j > 0u)
333       {
334         /* Divide each element of the row of the destination matrix    
335          * by the pivot element */
336         in1 = *pInT2;
337         *pInT2++ = in1 / in;
338
339         /* Decrement the loop counter */
340         j--;
341       }
342
343       /* Replace the rows with the sum of that row and a multiple of row i    
344        * so that each new element in column i above row i is zero.*/
345
346       /* Temporary pointers for input and destination matrices */
347       pInT1 = pIn;
348       pInT2 = pOut;
349
350       /* index used to check for pivot element */
351       i = 0u;
352
353       /* Loop over number of rows */
354       /*  to be replaced by the sum of that row and a multiple of row i */
355       k = numRows;
356
357       while(k > 0u)
358       {
359         /* Check for the pivot element */
360         if(i == l)
361         {
362           /* If the processing element is the pivot element,    
363              only the columns to the right are to be processed */
364           pInT1 += numCols - l;
365
366           pInT2 += numCols;
367         }
368         else
369         {
370           /* Element of the reference row */
371           in = *pInT1;
372
373           /* Working pointers for input and destination pivot rows */
374           pPRT_in = pPivotRowIn;
375           pPRT_pDst = pPivotRowDst;
376
377           /* Loop over the number of columns to the right of the pivot element,    
378              to replace the elements in the input matrix */
379           j = (numCols - l);
380
381           while(j > 0u)
382           {
383             /* Replace the element by the sum of that row    
384                and a multiple of the reference row  */
385             in1 = *pInT1;
386             *pInT1++ = in1 - (in * *pPRT_in++);
387
388             /* Decrement the loop counter */
389             j--;
390           }
391
392           /* Loop over the number of columns to    
393              replace the elements in the destination matrix */
394           j = numCols;
395
396           while(j > 0u)
397           {
398             /* Replace the element by the sum of that row    
399                and a multiple of the reference row  */
400             in1 = *pInT2;
401             *pInT2++ = in1 - (in * *pPRT_pDst++);
402
403             /* Decrement the loop counter */
404             j--;
405           }
406
407         }
408
409         /* Increment the temporary input pointer */
410         pInT1 = pInT1 + l;
411
412         /* Decrement the loop counter */
413         k--;
414
415         /* Increment the pivot index */
416         i++;
417       }
418
419       /* Increment the input pointer */
420       pIn++;
421
422       /* Decrement the loop counter */
423       loopCnt--;
424
425       /* Increment the index modifier */
426       l++;
427     }
428
429
430 #else
431
432   /* Run the below code for Cortex-M0 */
433
434   float32_t Xchg, in = 0.0f;                     /* Temporary input values  */
435   uint32_t i, rowCnt, flag = 0u, j, loopCnt, k, l;      /* loop counters */
436   arm_status status;                             /* status of matrix inverse */
437
438 #ifdef ARM_MATH_MATRIX_CHECK
439
440   /* Check for matrix mismatch condition */
441   if((pSrc->numRows != pSrc->numCols) || (pDst->numRows != pDst->numCols)
442      || (pSrc->numRows != pDst->numRows))
443   {
444     /* Set status as ARM_MATH_SIZE_MISMATCH */
445     status = ARM_MATH_SIZE_MISMATCH;
446   }
447   else
448 #endif /*      #ifdef ARM_MATH_MATRIX_CHECK    */
449   {
450
451     /*--------------------------------------------------------------------------------------------------------------       
452          * Matrix Inverse can be solved using elementary row operations.        
453          *        
454          *      Gauss-Jordan Method:       
455          *                     
456          *         1. First combine the identity matrix and the input matrix separated by a bar to form an        
457          *        augmented matrix as follows:        
458          *                                      _  _          _     _      _   _         _             _       
459          *                                         |  |  a11  a12  | | | 1   0  |   |       |  X11 X12  |         
460          *                                         |  |            | | |        |   |   =   |           |        
461          *                                         |_ |_ a21  a22 _| | |_0   1 _|  _|       |_ X21 X21 _|       
462          *                                                
463          *              2. In our implementation, pDst Matrix is used as identity matrix.    
464          *       
465          *              3. Begin with the first row. Let i = 1.       
466          *       
467          *          4. Check to see if the pivot for row i is zero.       
468          *                 The pivot is the element of the main diagonal that is on the current row.       
469          *                 For instance, if working with row i, then the pivot element is aii.       
470          *                 If the pivot is zero, exchange that row with a row below it that does not        
471          *                 contain a zero in column i. If this is not possible, then an inverse        
472          *                 to that matrix does not exist.       
473          *             
474          *          5. Divide every element of row i by the pivot.       
475          *             
476          *          6. For every row below and  row i, replace that row with the sum of that row and        
477          *                 a multiple of row i so that each new element in column i below row i is zero.       
478          *             
479          *          7. Move to the next row and column and repeat steps 2 through 5 until you have zeros       
480          *                 for every element below and above the main diagonal.        
481          *                                        
482          *              8. Now an identical matrix is formed to the left of the bar(input matrix, src).       
483          *                 Therefore, the matrix to the right of the bar is our solution(dst matrix, dst).         
484          *----------------------------------------------------------------------------------------------------------------*/
485
486     /* Working pointer for destination matrix */
487     pInT2 = pOut;
488
489     /* Loop over the number of rows */
490     rowCnt = numRows;
491
492     /* Making the destination matrix as identity matrix */
493     while(rowCnt > 0u)
494     {
495       /* Writing all zeroes in lower triangle of the destination matrix */
496       j = numRows - rowCnt;
497       while(j > 0u)
498       {
499         *pInT2++ = 0.0f;
500         j--;
501       }
502
503       /* Writing all ones in the diagonal of the destination matrix */
504       *pInT2++ = 1.0f;
505
506       /* Writing all zeroes in upper triangle of the destination matrix */
507       j = rowCnt - 1u;
508       while(j > 0u)
509       {
510         *pInT2++ = 0.0f;
511         j--;
512       }
513
514       /* Decrement the loop counter */
515       rowCnt--;
516     }
517
518     /* Loop over the number of columns of the input matrix.     
519        All the elements in each column are processed by the row operations */
520     loopCnt = numCols;
521
522     /* Index modifier to navigate through the columns */
523     l = 0u;
524     //for(loopCnt = 0u; loopCnt < numCols; loopCnt++)   
525     while(loopCnt > 0u)
526     {
527       /* Check if the pivot element is zero..    
528        * If it is zero then interchange the row with non zero row below.   
529        * If there is no non zero element to replace in the rows below,   
530        * then the matrix is Singular. */
531
532       /* Working pointer for the input matrix that points     
533        * to the pivot element of the particular row  */
534       pInT1 = pIn + (l * numCols);
535
536       /* Working pointer for the destination matrix that points     
537        * to the pivot element of the particular row  */
538       pInT3 = pOut + (l * numCols);
539
540       /* Temporary variable to hold the pivot value */
541       in = *pInT1;
542
543       /* Destination pointer modifier */
544       k = 1u;
545
546       /* Check if the pivot element is zero */
547       if(*pInT1 == 0.0f)
548       {
549         /* Loop over the number rows present below */
550         for (i = (l + 1u); i < numRows; i++)
551         {
552           /* Update the input and destination pointers */
553           pInT2 = pInT1 + (numCols * l);
554           pInT4 = pInT3 + (numCols * k);
555
556           /* Check if there is a non zero pivot element to     
557            * replace in the rows below */
558           if(*pInT2 != 0.0f)
559           {
560             /* Loop over number of columns     
561              * to the right of the pilot element */
562             for (j = 0u; j < (numCols - l); j++)
563             {
564               /* Exchange the row elements of the input matrix */
565               Xchg = *pInT2;
566               *pInT2++ = *pInT1;
567               *pInT1++ = Xchg;
568             }
569
570             for (j = 0u; j < numCols; j++)
571             {
572               Xchg = *pInT4;
573               *pInT4++ = *pInT3;
574               *pInT3++ = Xchg;
575             }
576
577             /* Flag to indicate whether exchange is done or not */
578             flag = 1u;
579
580             /* Break after exchange is done */
581             break;
582           }
583
584           /* Update the destination pointer modifier */
585           k++;
586         }
587       }
588
589       /* Update the status if the matrix is singular */
590       if((flag != 1u) && (in == 0.0f))
591       {
592         status = ARM_MATH_SINGULAR;
593
594         break;
595       }
596
597       /* Points to the pivot row of input and destination matrices */
598       pPivotRowIn = pIn + (l * numCols);
599       pPivotRowDst = pOut + (l * numCols);
600
601       /* Temporary pointers to the pivot row pointers */
602       pInT1 = pPivotRowIn;
603       pInT2 = pPivotRowDst;
604
605       /* Pivot element of the row */
606       in = *(pIn + (l * numCols));
607
608       /* Loop over number of columns     
609        * to the right of the pilot element */
610       for (j = 0u; j < (numCols - l); j++)
611       {
612         /* Divide each element of the row of the input matrix     
613          * by the pivot element */
614         *pInT1 = *pInT1 / in;
615         pInT1++;
616       }
617       for (j = 0u; j < numCols; j++)
618       {
619         /* Divide each element of the row of the destination matrix     
620          * by the pivot element */
621         *pInT2 = *pInT2 / in;
622         pInT2++;
623       }
624
625       /* Replace the rows with the sum of that row and a multiple of row i     
626        * so that each new element in column i above row i is zero.*/
627
628       /* Temporary pointers for input and destination matrices */
629       pInT1 = pIn;
630       pInT2 = pOut;
631
632       for (i = 0u; i < numRows; i++)
633       {
634         /* Check for the pivot element */
635         if(i == l)
636         {
637           /* If the processing element is the pivot element,     
638              only the columns to the right are to be processed */
639           pInT1 += numCols - l;
640           pInT2 += numCols;
641         }
642         else
643         {
644           /* Element of the reference row */
645           in = *pInT1;
646
647           /* Working pointers for input and destination pivot rows */
648           pPRT_in = pPivotRowIn;
649           pPRT_pDst = pPivotRowDst;
650
651           /* Loop over the number of columns to the right of the pivot element,     
652              to replace the elements in the input matrix */
653           for (j = 0u; j < (numCols - l); j++)
654           {
655             /* Replace the element by the sum of that row     
656                and a multiple of the reference row  */
657             *pInT1 = *pInT1 - (in * *pPRT_in++);
658             pInT1++;
659           }
660           /* Loop over the number of columns to     
661              replace the elements in the destination matrix */
662           for (j = 0u; j < numCols; j++)
663           {
664             /* Replace the element by the sum of that row     
665                and a multiple of the reference row  */
666             *pInT2 = *pInT2 - (in * *pPRT_pDst++);
667             pInT2++;
668           }
669
670         }
671         /* Increment the temporary input pointer */
672         pInT1 = pInT1 + l;
673       }
674       /* Increment the input pointer */
675       pIn++;
676
677       /* Decrement the loop counter */
678       loopCnt--;
679       /* Increment the index modifier */
680       l++;
681     }
682
683
684 #endif /* #ifndef ARM_MATH_CM0_FAMILY */
685
686     /* Set status as ARM_MATH_SUCCESS */
687     status = ARM_MATH_SUCCESS;
688
689     if((flag != 1u) && (in == 0.0f))
690     {
691       status = ARM_MATH_SINGULAR;
692     }
693   }
694   /* Return to application */
695   return (status);
696 }
697
698 /**    
699  * @} end of MatrixInv group    
700  */