Subversion Repositories DashDisplay

Rev

Blame | Last modification | View Log | Download | RSS feed

  1. /* ----------------------------------------------------------------------      
  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_cmplx_mult_q31.c      
  9. *      
  10. * Description:  Floating-point matrix multiplication.      
  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. #include "arm_math.h"
  41.  
  42. /**    
  43.  * @ingroup groupMatrix    
  44.  */
  45.  
  46. /**      
  47.  * @addtogroup CmplxMatrixMult      
  48.  * @{      
  49.  */
  50.  
  51. /**      
  52.  * @brief Q31 Complex matrix multiplication      
  53.  * @param[in]       *pSrcA points to the first input complex matrix structure      
  54.  * @param[in]       *pSrcB points to the second input complex matrix structure      
  55.  * @param[out]      *pDst points to output complex matrix structure      
  56.  * @return              The function returns either      
  57.  * <code>ARM_MATH_SIZE_MISMATCH</code> or <code>ARM_MATH_SUCCESS</code> based on the outcome of size checking.      
  58.  *      
  59.  * @details      
  60.  * <b>Scaling and Overflow Behavior:</b>      
  61.  *      
  62.  * \par      
  63.  * The function is implemented using an internal 64-bit accumulator.      
  64.  * The accumulator has a 2.62 format and maintains full precision of the intermediate      
  65.  * multiplication results but provides only a single guard bit. There is no saturation      
  66.  * on intermediate additions. Thus, if the accumulator overflows it wraps around and      
  67.  * distorts the result. The input signals should be scaled down to avoid intermediate      
  68.  * overflows. The input is thus scaled down by log2(numColsA) bits      
  69.  * to avoid overflows, as a total of numColsA additions are performed internally.      
  70.  * The 2.62 accumulator is right shifted by 31 bits and saturated to 1.31 format to yield the final result.      
  71.  *      
  72.  *      
  73.  */
  74.  
  75. arm_status arm_mat_cmplx_mult_q31(
  76.   const arm_matrix_instance_q31 * pSrcA,
  77.   const arm_matrix_instance_q31 * pSrcB,
  78.   arm_matrix_instance_q31 * pDst)
  79. {
  80.   q31_t *pIn1 = pSrcA->pData;                    /* input data matrix pointer A */
  81.   q31_t *pIn2 = pSrcB->pData;                    /* input data matrix pointer B */
  82.   q31_t *pInA = pSrcA->pData;                    /* input data matrix pointer A  */
  83.   q31_t *pOut = pDst->pData;                     /* output data matrix pointer */
  84.   q31_t *px;                                     /* Temporary output data matrix pointer */
  85.   uint16_t numRowsA = pSrcA->numRows;            /* number of rows of input matrix A */
  86.   uint16_t numColsB = pSrcB->numCols;            /* number of columns of input matrix B */
  87.   uint16_t numColsA = pSrcA->numCols;            /* number of columns of input matrix A */
  88.   q63_t sumReal1, sumImag1;                      /* accumulator */
  89.   q31_t a0, b0, c0, d0;
  90.   q31_t a1, b1, c1, d1;
  91.  
  92.  
  93.   /* Run the below code for Cortex-M4 and Cortex-M3 */
  94.  
  95.   uint16_t col, i = 0u, j, row = numRowsA, colCnt;      /* loop counters */
  96.   arm_status status;                             /* status of matrix multiplication */
  97.  
  98. #ifdef ARM_MATH_MATRIX_CHECK
  99.  
  100.  
  101.   /* Check for matrix mismatch condition */
  102.   if((pSrcA->numCols != pSrcB->numRows) ||
  103.      (pSrcA->numRows != pDst->numRows) || (pSrcB->numCols != pDst->numCols))
  104.   {
  105.  
  106.     /* Set status as ARM_MATH_SIZE_MISMATCH */
  107.     status = ARM_MATH_SIZE_MISMATCH;
  108.   }
  109.   else
  110. #endif /*      #ifdef ARM_MATH_MATRIX_CHECK    */
  111.  
  112.   {
  113.     /* The following loop performs the dot-product of each row in pSrcA with each column in pSrcB */
  114.     /* row loop */
  115.     do
  116.     {
  117.       /* Output pointer is set to starting address of the row being processed */
  118.       px = pOut + 2 * i;
  119.  
  120.       /* For every row wise process, the column loop counter is to be initiated */
  121.       col = numColsB;
  122.  
  123.       /* For every row wise process, the pIn2 pointer is set    
  124.        ** to the starting address of the pSrcB data */
  125.       pIn2 = pSrcB->pData;
  126.  
  127.       j = 0u;
  128.  
  129.       /* column loop */
  130.       do
  131.       {
  132.         /* Set the variable sum, that acts as accumulator, to zero */
  133.         sumReal1 = 0.0;
  134.         sumImag1 = 0.0;
  135.  
  136.         /* Initiate the pointer pIn1 to point to the starting address of the column being processed */
  137.         pIn1 = pInA;
  138.  
  139.         /* Apply loop unrolling and compute 4 MACs simultaneously. */
  140.         colCnt = numColsA >> 2;
  141.  
  142.         /* matrix multiplication        */
  143.         while(colCnt > 0u)
  144.         {
  145.  
  146.           /* Reading real part of complex matrix A */
  147.           a0 = *pIn1;
  148.  
  149.           /* Reading real part of complex matrix B */
  150.           c0 = *pIn2;
  151.  
  152.           /* Reading imaginary part of complex matrix A */
  153.           b0 = *(pIn1 + 1u);
  154.  
  155.           /* Reading imaginary part of complex matrix B */
  156.           d0 = *(pIn2 + 1u);
  157.  
  158.           /* Multiply and Accumlates */
  159.           sumReal1 += (q63_t) a0 *c0;
  160.           sumImag1 += (q63_t) b0 *c0;
  161.  
  162.           /* update pointers */
  163.           pIn1 += 2u;
  164.           pIn2 += 2 * numColsB;
  165.  
  166.           /* Multiply and Accumlates */
  167.           sumReal1 -= (q63_t) b0 *d0;
  168.           sumImag1 += (q63_t) a0 *d0;
  169.  
  170.           /* c(m,n) = a(1,1)*b(1,1) + a(1,2) * b(2,1) + .... + a(m,p)*b(p,n) */
  171.  
  172.           /* read real and imag values from pSrcA and pSrcB buffer */
  173.           a1 = *pIn1;
  174.           c1 = *pIn2;
  175.           b1 = *(pIn1 + 1u);
  176.           d1 = *(pIn2 + 1u);
  177.  
  178.           /* Multiply and Accumlates */
  179.           sumReal1 += (q63_t) a1 *c1;
  180.           sumImag1 += (q63_t) b1 *c1;
  181.  
  182.           /* update pointers */
  183.           pIn1 += 2u;
  184.           pIn2 += 2 * numColsB;
  185.  
  186.           /* Multiply and Accumlates */
  187.           sumReal1 -= (q63_t) b1 *d1;
  188.           sumImag1 += (q63_t) a1 *d1;
  189.  
  190.           a0 = *pIn1;
  191.           c0 = *pIn2;
  192.  
  193.           b0 = *(pIn1 + 1u);
  194.           d0 = *(pIn2 + 1u);
  195.  
  196.           /* Multiply and Accumlates */
  197.           sumReal1 += (q63_t) a0 *c0;
  198.           sumImag1 += (q63_t) b0 *c0;
  199.  
  200.           /* update pointers */
  201.           pIn1 += 2u;
  202.           pIn2 += 2 * numColsB;
  203.  
  204.           /* Multiply and Accumlates */
  205.           sumReal1 -= (q63_t) b0 *d0;
  206.           sumImag1 += (q63_t) a0 *d0;
  207.  
  208.           /* c(m,n) = a(1,1)*b(1,1) + a(1,2) * b(2,1) + .... + a(m,p)*b(p,n) */
  209.  
  210.           a1 = *pIn1;
  211.           c1 = *pIn2;
  212.  
  213.           b1 = *(pIn1 + 1u);
  214.           d1 = *(pIn2 + 1u);
  215.  
  216.           /* Multiply and Accumlates */
  217.           sumReal1 += (q63_t) a1 *c1;
  218.           sumImag1 += (q63_t) b1 *c1;
  219.  
  220.           /* update pointers */
  221.           pIn1 += 2u;
  222.           pIn2 += 2 * numColsB;
  223.  
  224.           /* Multiply and Accumlates */
  225.           sumReal1 -= (q63_t) b1 *d1;
  226.           sumImag1 += (q63_t) a1 *d1;
  227.  
  228.           /* Decrement the loop count */
  229.           colCnt--;
  230.         }
  231.  
  232.         /* If the columns of pSrcA is not a multiple of 4, compute any remaining MACs here.    
  233.          ** No loop unrolling is used. */
  234.         colCnt = numColsA % 0x4u;
  235.  
  236.         while(colCnt > 0u)
  237.         {
  238.           /* c(m,n) = a(1,1)*b(1,1) + a(1,2) * b(2,1) + .... + a(m,p)*b(p,n) */
  239.           a1 = *pIn1;
  240.           c1 = *pIn2;
  241.  
  242.           b1 = *(pIn1 + 1u);
  243.           d1 = *(pIn2 + 1u);
  244.  
  245.           /* Multiply and Accumlates */
  246.           sumReal1 += (q63_t) a1 *c1;
  247.           sumImag1 += (q63_t) b1 *c1;
  248.  
  249.           /* update pointers */
  250.           pIn1 += 2u;
  251.           pIn2 += 2 * numColsB;
  252.  
  253.           /* Multiply and Accumlates */
  254.           sumReal1 -= (q63_t) b1 *d1;
  255.           sumImag1 += (q63_t) a1 *d1;
  256.  
  257.           /* Decrement the loop counter */
  258.           colCnt--;
  259.         }
  260.  
  261.         /* Store the result in the destination buffer */
  262.         *px++ = (q31_t) clip_q63_to_q31(sumReal1 >> 31);
  263.         *px++ = (q31_t) clip_q63_to_q31(sumImag1 >> 31);
  264.        
  265.         /* Update the pointer pIn2 to point to the  starting address of the next column */
  266.         j++;
  267.         pIn2 = pSrcB->pData + 2u * j;
  268.  
  269.         /* Decrement the column loop counter */
  270.         col--;
  271.  
  272.       } while(col > 0u);
  273.  
  274.       /* Update the pointer pInA to point to the  starting address of the next row */
  275.       i = i + numColsB;
  276.       pInA = pInA + 2 * numColsA;
  277.  
  278.       /* Decrement the row loop counter */
  279.       row--;
  280.  
  281.     } while(row > 0u);
  282.  
  283.     /* Set status as ARM_MATH_SUCCESS */
  284.     status = ARM_MATH_SUCCESS;
  285.   }
  286.  
  287.   /* Return to application */
  288.   return (status);
  289. }
  290.  
  291. /**    
  292.  * @} end of MatrixMult group    
  293.  */
  294.