精华内容
下载资源
问答
  • view plaincopy [cpp]view plaincopy 测试程序如下: [cpp]view plaincopy (说明:下面代码分别包含了三个矩阵相乘的函数) [cpp]view plaincopy #include#include #include"complex.h" voidmatrix_mul( constdouble*...

    来至http://blog..net/gtatcs/article/details/8887082;

    下面代码参考TI的实现:

    [cpp]view

    plaincopy

    /*NAME*/

    /*DSPF_dp_mat_mul_cplx--Complexmatrixmultiplication*/

    /**/

    /*USAGE*/

    /**/

    /*ThisroutinehasthefollowingCprototype:*/

    /**/

    /*voidDSPF_dp_mat_mul_cplx(*/

    /*constdouble*x,*/

    /*intr1,*/

    /*intc1,*/

    /*constdouble*y,*/

    /*intc2,*/

    /*double*restrictr*/

    /*)*/

    /**/

    /*x[2*r1*c1]:Inputmatrixcontainingr1*c1complex*/

    /*floatingpointnumbershavingr1rowsandc1*/

    /*columnsofcomplexnumbers.*/

    /*r1:No.ofrowsinmatrixx.*/

    /*c1:No.ofcolumnsinmatrixx.*/

    /*Alsono.ofrowsinmatrixy.*/

    /*y[2*c1*c2]:Inputmatrixcontainingc1*c2complex*/

    /*floatingpointnumbershavingc1rowsandc2*/

    /*columnsofcomplexnumbers.*/

    /*c2:No.ofcolumnsinmatrixy.*/

    /*r[2*r1*c2]:Outputmatrixofc1*c2complexfloating*/

    /*pointnumbershavingc1rowsandc2columnsof*/

    /*complexnumbers.*/

    /**/

    /*Complexnumbersarestoredconsecutivelywith*/

    /*realvaluesstoredinevenpositionsand*/

    /*imaginaryvaluesinoddpositions.*/

    /**/

    /*DESCRIPTION*/

    /**/

    /*Thisfunctioncomputestheexpression"r=x*y"forthe*/

    /*matricesxandy.Thecolumnardimensionofxmustmatchtherow*/

    /*dimensionofy.Theresultingmatrixhasthesamenumberofrows*/

    /*asxandthesamenumberofcolumnsasy.*/

    /**/

    /*Eachelementofthematrixisassumedtobecomplexnumberswith*/

    /*Realvaluesarestoredinevenpositionsandimaginary*/

    /*valuesinoddpositions.*/

    /**/

    /*TECHNIQUES*/

    /**/

    /*1.Innermostloopisunrolledtwice.*/

    /*2.Outermostloopisexecutedinparallelwithinnnerloops.*/

    /**/

    /*ASSUMPTIONS*/

    /**/

    /*1.r1,r2>=1,c1shouldbeamultipleof2and>=2.*/

    /*2.xshouldbepaddedwith6words*/

    /**/

    /**/

    /*CCODE*/

    /**/

    voidDSPF_dp_mat_mul_cplx(constdouble*x,intr1,intc1,

    constdouble*y,intc2,double*r)

    {

    doublereal,imag;

    inti,j,k;

    for(i=0;i

    {

    for(j=0;j

    {

    real=0;

    imag=0;

    for(k=0;k

    {

    real+=(x[i*2*c1+2*k]*y[k*2*c2+2*j]

    -x[i*2*c1+2*k+1]*y[k*2*c2+2*j+1]);

    imag+=(x[i*2*c1+2*k]*y[k*2*c2+2*j+1]

    +x[i*2*c1+2*k+1]*y[k*2*c2+2*j]);

    }

    r[i*2*c2+2*j]=real;

    r[i*2*c2+2*j+1]=imag;

    }

    }

    }

    /*NOTES*/

    /**/

    /*1.Realvaluesarestoredinevenwordpositionsandimaginary*/

    /*valuesinoddpositions.*/

    /*2.Endian:ThiscodeisLITTLEENDIAN.*/

    /*3.Interruptibility:Thiscodeisinterrupttolerantbutnot*/

    /*interruptible.*/

    /**/

    /*CYCLES*/

    /**/

    /*8*r1*c1*c2'+18*(r1*c2)+40WHEREc2'=2*ceil(c2/2)*/

    /*Whenr1=3,c1=4,c2=4,cycles=640*/

    /*Whenr1=4,c1=4,c2=5,cycles=1040*/

    /**/

    /*CODESIZE*/

    /**/

    /*83***ytes*/

    /*------------------------------------------------------------------------*/

    /*Copyright(c)2004TexasInstruments,Incorporated.*/

    /*AllRightsReserved.*/

    /*========================================================================*/

    #ifndefDSPF_DP_MAT_MUL_CPLX_H_

    #defineDSPF_DP_MAT_MUL_CPLX_H_1

    voidDSPF_dp_mat_mul_cplx(

    constdouble*x,

    intr1,

    intc1,

    constdouble*y,

    intc2,

    double*restrictr

    );

    #endif

    /*========================================================================*/

    /*Endoffile:DSPF_dp_mat_mul_cplx.h*/

    /*------------------------------------------------------------------------*/

    [cpp]view

    plaincopy

    [cpp]view

    plaincopy

    测试程序如下:

    [cpp]view

    plaincopy

    (说明:下面代码分别包含了三个矩阵相乘的函数)

    [cpp]view

    plaincopy

    #include#include

    #include"complex.h"

    voidmatrix_mul(

    constdouble*x,

    intr1,

    intc1,

    constdouble*y,

    intc2,

    double*r

    )

    {

    inti,j,k;

    for(i=0;i

    {

    for(j=0;j

    {

    for(k=0;k

    {

    //r[i*r1+j]=r[i*r1+j]+x[i*r1+k]*y[k*c1+j];

    r[i*c2+j]=r[i*c2+j]+x[i*c1+k]*y[k*c2+j];

    }

    }

    }

    }

    voidmatrix_mul_cplx(

    constComplex*x,

    intr1,

    intc1,

    constComplex*y,

    intc2,

    Complex*r

    )

    {

    inti,j,k;

    Complextemp;

    temp.real=0;

    temp.image=0;

    for(i=0;i

    {

    for(j=0;j

    {

    for(k=0;k

    {

    //r[i*r1+j]=r[i*r1+j]+x[i*r1+k]*y[k*c1+j];

    temp=complex_mul(x[i*c1+k],y[k*c2+j]);

    r[i*c2+j]=complex_add(r[i*c2+j],temp);

    }

    }

    }

    }

    voidDSPF_dp_mat_mul_cplx(constdouble*x,intr1,intc1,

    constdouble*y,intc2,double*r)

    {

    doublereal,imag;

    inti,j,k;

    for(i=0;i

    {

    for(j=0;j

    {

    real=0;

    imag=0;

    for(k=0;k

    {

    real+=(x[i*2*c1+2*k]*y[k*2*c2+2*j]

    -x[i*2*c1+2*k+1]*y[k*2*c2+2*j+1]);

    imag+=(x[i*2*c1+2*k]*y[k*2*c2+2*j+1]

    +x[i*2*c1+2*k+1]*y[k*2*c2+2*j]);

    }

    r[i*2*c2+2*j]=real;

    r[i*2*c2+2*j+1]=imag;

    }

    }

    }

    intmain()

    {

    //matrix_multest

    /*

    doublex[9]={1,2,3,4,5,6,7,8,9};//矩阵按行存放

    intr1=3,c1=3,c2=1;

    doubley[3]={10,11,12};

    doubler[3]={0};

    inti,j;

    //printf("r:%lf\n%lf\n",r[0],r[1]);

    printf("r=:\n");

    for(i=0;i

    {

    for(j=0;j

    {

    //printf("%lf\t",r[i*r1+j]);

    printf("%lf\t",r[i*c2+j]);

    }

    printf("\n");

    }

    matrix_mul(x,r1,c1,y,c2,r);

    printf("x*y=:\n");

    for(i=0;i

    展开全文
  • 此函数计算输入复矩阵 A(2,2,m) 的累积矩阵乘法此函数应使用 Maltab Mex 函数进行编译。 米 = 1000; % 矩阵数A = rand(2,2,m) + 1i * rand(2,2,m); % 示例矩阵; B = cM(A); % B 是一个 2x2 矩阵 此函数比等效的 for...
  • 复数矩阵乘法C语言实现

    千次阅读 2013-05-08 15:13:22
    来至http://blog.csdn.net/gtatcs/article/details/8887082; 下面代码参考TI的实现: [cpp] view plaincopy "font-size:14px;">/* NAME

    来至http://blog.csdn.net/gtatcs/article/details/8887082;

    下面代码参考TI的实现:

    1. <span style="font-size:14px;">/* NAME                                                                    */  
    2. /*     DSPF_dp_mat_mul_cplx -- Complex matrix multiplication                    */  
    3. /*                                                                         */  
    4. /*  USAGE                                                                   */  
    5. /*                                                                          */  
    6. /*       This routine has the following C prototype:                        */  
    7. /*                                                                          */  
    8. /*       void DSPF_dp_mat_mul_cplx(                                              */  
    9. /*                              const double* x,                            */  
    10. /*                              int r1,                                     */  
    11. /*                              int c1,                                     */  
    12. /*                              const double* y,                            */  
    13. /*                              int c2,                                     */  
    14. /*                              double* restrict r                          */  
    15. /*                           )                                              */  
    16. /*                                                                          */  
    17. /*            x[2*r1*c1]:   Input matrix containing r1*c1 complex           */  
    18. /*                          floating point numbers having r1 rows and c1    */  
    19. /*                          columns of complex numbers.                     */  
    20. /*            r1        :   No. of rows in matrix x.                        */  
    21. /*            c1        :   No. of columns in matrix x.                     */  
    22. /*                          Also no. of rows in matrix y.                   */  
    23. /*            y[2*c1*c2]:   Input matrix containing c1*c2 complex           */  
    24. /*                          floating point numbers having c1 rows and c2    */  
    25. /*                          columns of complex numbers.                     */  
    26. /*            c2        :   No. of columns in matrix y.                     */  
    27. /*            r[2*r1*c2]:   Output matrix of c1*c2 complex floating         */  
    28. /*                          point numbers having c1 rows and c2 columns of  */  
    29. /*                          complex numbers.                                */  
    30. /*                                                                          */  
    31. /*                          Complex numbers are stored consecutively with   */  
    32. /*                          real values stored in even positions and        */  
    33. /*                          imaginary values in odd positions.              */  
    34. /*                                                                          */  
    35. /*  DESCRIPTION                                                             */  
    36. /*                                                                          */  
    37. /*        This function computes the expression "r = x * y" for the         */  
    38. /*        matrices x and y. The columnar dimension of x must match the row  */  
    39. /*        dimension of y. The resulting matrix has the same number of rows  */  
    40. /*        as x and the same number of columns as y.                         */  
    41. /*                                                                          */  
    42. /*        Each element of the matrix is assumed to be complex numbers with  */  
    43. /*        Real values are stored in even positions and imaginary            */  
    44. /*        values in odd positions.                                          */  
    45. /*                                                                          */  
    46. /*  TECHNIQUES                                                              */  
    47. /*                                                                          */  
    48. /*        1. Innermost loop is unrolled twice.                              */  
    49. /*        2. Outermost loop is executed in parallel with innner loops.      */  
    50. /*                                                                          */  
    51. /*  ASSUMPTIONS                                                             */  
    52. /*                                                                          */  
    53. /*        1. r1,r2>=1,c1 should be a multiple of 2 and >=2.                 */  
    54. /*        2. x should be padded with 6 words                                */  
    55. /*                                                                          */  
    56. /*                                                                          */  
    57. /*  C CODE                                                                  */  
    58. /*                                                                          */  
    59. void DSPF_dp_mat_mul_cplx(const double* x, int r1, int c1,   
    60.                           const double* y, int c2, double* r)  
    61. {                                                                   
    62.     double real, imag;                                              
    63.     int i, j, k;                                                    
    64.                                                         
    65.     for(i = 0; i < r1; i++)                                        
    66.     {                                                               
    67.         for(j = 0; j < c2; j++)                                       
    68.         {                                                             
    69.             real=0;                                                     
    70.             imag=0;                                                     
    71.                                                                  
    72.             for(k = 0; k < c1; k++)                                     
    73.             {                                                           
    74.                 real += (x[i*2*c1 + 2*k]*y[k*2*c2 + 2*j]                    
    75.                 -x[i*2*c1 + 2*k + 1] * y[k*2*c2 + 2*j + 1]);                
    76.                                                                      
    77.                 imag+=(x[i*2*c1 + 2*k] * y[k*2*c2 + 2*j + 1]               
    78.                 + x[i*2*c1 + 2*k + 1] * y[k*2*c2 + 2*j]);                
    79.             }                                                           
    80.             r[i*2*c2 + 2*j] = real;                                    
    81.             r[i*2*c2 + 2*j + 1] = imag;    
    82.         }  
    83.     }  
    84. }    
    85. /*  NOTES                                                                   */  
    86. /*                                                                          */  
    87. /*        1. Real values are stored in even word positions and imaginary    */  
    88. /*           values in odd positions.                                       */  
    89. /*        2. Endian: This code is LITTLE ENDIAN.                            */  
    90. /*        3. Interruptibility: This code is interrupt tolerant but not      */  
    91. /*           interruptible.                                                 */  
    92. /*                                                                          */  
    93. /*  CYCLES                                                                  */  
    94. /*                                                                          */  
    95. /*        8*r1*c1*c2'+ 18*(r1*c2)+40 WHERE c2'=2*ceil(c2/2)                 */  
    96. /*        When r1=3, c1=4, c2=4, cycles = 640                               */  
    97. /*        When r1=4, c1=4, c2=5, cycles = 1040                              */  
    98. /*                                                                          */  
    99. /*  CODESIZE                                                                */  
    100. /*                                                                          */  
    101. /*        832 bytes                                                         */  
    102. /* ------------------------------------------------------------------------ */  
    103. /*            Copyright (c) 2004 Texas Instruments, Incorporated.           */  
    104. /*                           All Rights Reserved.                           */  
    105. /* ======================================================================== */  
    106. #ifndef DSPF_DP_MAT_MUL_CPLX_H_  
    107. #define DSPF_DP_MAT_MUL_CPLX_H_ 1  
    108.   
    109. void DSPF_dp_mat_mul_cplx(  
    110.                        const double* x,  
    111.                        int r1,  
    112.                        int c1,  
    113.                        const double* y,  
    114.                        int c2,  
    115.                        double* restrict r  
    116.                     );  
    117.   
    118. #endif  
    119. /* ======================================================================== */  
    120. /*  End of file:  DSPF_dp_mat_mul_cplx.h                                         */  
    121. /* ------------------------------------------------------------------------ */</span>  
    1. <span style="font-size:14px;">  
    2. </span>  
    1. <span style="font-size:14px;">测试程序如下:</span>  
    1. <span style="font-size:14px;">(说明:下面代码分别包含了三个矩阵相乘的函数)</span>  
    1. <pre name="code" class="cpp"><span style="font-size:14px;">#include <stdio.h></span></pre><pre name="code" class="cpp"><span style="font-size:14px;">#include <stdlib.h>  
    2. #include "complex.h"  
    3.   
    4. void matrix_mul(  
    5.                      const double* x,  
    6.                      int r1,  
    7.                      int c1,  
    8.                      const double* y,  
    9.                      int c2,  
    10.                      double* r  
    11.                     )  
    12. {  
    13.     int i,j,k;  
    14.     for(i=0;i<r1;++i)  
    15.     {  
    16.         for(j=0;j<c2;++j)  
    17.         {  
    18.             for(k=0;k<c1;++k)  
    19.             {  
    20.                 //r[i*r1 + j] =r[i*r1 + j] + x[i*r1 + k] * y[k*c1 + j];  
    21.                 r[i*c2 + j] =r[i*c2 + j] + x[i*c1 + k] * y[k*c2 + j];  
    22.             }  
    23.         }  
    24.     }  
    25. }  
    26.   
    27. void matrix_mul_cplx(  
    28.                 const Complex* x,  
    29.                 int r1,  
    30.                 int c1,  
    31.                 const Complex* y,  
    32.                 int c2,  
    33.                 Complex* r  
    34.                 )  
    35. {  
    36.     int i,j,k;  
    37.     Complex temp;  
    38.     temp.real = 0;  
    39.     temp.image = 0;  
    40.     for(i=0;i<r1;++i)  
    41.     {  
    42.         for(j=0;j<c2;++j)  
    43.         {  
    44.             for(k=0;k<c1;++k)  
    45.             {  
    46.                 //r[i*r1 + j] =r[i*r1 + j] + x[i*r1 + k] * y[k*c1 + j];  
    47.                 temp = complex_mul(x[i*c1 + k], y[k*c2 + j]);  
    48.                 r[i*c2 + j] =complex_add(r[i*c2 + j], temp);  
    49.             }  
    50.         }  
    51.     }  
    52. }  
    53.   
    54. void DSPF_dp_mat_mul_cplx(const double* x, int r1, int c1,   
    55.                           const double* y, int c2, double* r)  
    56. {                                                                   
    57.     double real, imag;                                              
    58.     int i, j, k;                                                    
    59.                                                         
    60.     for(i = 0; i < r1; i++)                                        
    61.     {                                                               
    62.         for(j = 0; j < c2; j++)                                       
    63.         {                                                             
    64.             real=0;                                                     
    65.             imag=0;                                                     
    66.                                                                  
    67.             for(k = 0; k < c1; k++)                                     
    68.             {                                                           
    69.                 real += (x[i*2*c1 + 2*k]*y[k*2*c2 + 2*j]                    
    70.                 -x[i*2*c1 + 2*k + 1] * y[k*2*c2 + 2*j + 1]);                
    71.                                                                      
    72.                 imag+=(x[i*2*c1 + 2*k] * y[k*2*c2 + 2*j + 1]               
    73.                 + x[i*2*c1 + 2*k + 1] * y[k*2*c2 + 2*j]);                
    74.             }                                                           
    75.             r[i*2*c2 + 2*j] = real;                                    
    76.             r[i*2*c2 + 2*j + 1] = imag;    
    77.         }  
    78.     }  
    79. }       
    80.           
    81. int main()  
    82. {  
    83.     //matrix_mul test  
    84.     /* 
    85.     double x[9] = {1,2,3,4,5,6,7,8,9};  //矩阵按行存放 
    86.     int r1 = 3, c1 = 3,c2 = 1; 
    87.     double y[3] = {10,11,12}; 
    88.     double r[3] = {0}; 
    89.     int i,j; 
    90.     //printf("r  :  %lf\n%lf\n",r[0],r[1]); 
    91.  
    92.     printf("r = : \n"); 
    93.     for(i=0;i<r1;++i) 
    94.     { 
    95.         for(j=0;j<c2;++j) 
    96.         { 
    97.             //printf("%lf\t",r[i*r1 + j]); 
    98.             printf("%lf\t",r[i*c2 + j]); 
    99.         } 
    100.         printf("\n"); 
    101.     } 
    102.  
    103.     matrix_mul(x,r1,c1,y,c2,r); 
    104.     printf("x*y = : \n"); 
    105.     for(i=0;i<r1;++i) 
    106.     { 
    107.         for(j=0;j<c2;++j) 
    108.         { 
    109.             printf("%lf\t",r[i*c2 + j]); 
    110.         } 
    111.         printf("\n"); 
    112.     } 
    113. */  
    114.     /* 
    115.     //matrix_mul_cplx test 
    116.     Complex x[4] = {{1,2}, {3,4}, {5,6}, {7,8}};    //矩阵按行存放 
    117.     int r1 = 2, c1 = 2,c2 = 1; 
    118.     Complex y[2] = {{1,2}, {3,4}}; 
    119.     Complex r[2] = {{0,0}, {0,0}}; 
    120.     int i,j; 
    121.     //printf("r  :  %lf\n%lf\n",r[0],r[1]); 
    122.      
    123.     printf("r = : \n"); 
    124.     for(i=0;i<r1;++i) 
    125.     { 
    126.         for(j=0;j<c2;++j) 
    127.         { 
    128.             printf("%lf + j*%lf\t",r[i*c2 + j].real , r[i*c2 + j].image); 
    129.         } 
    130.         printf("\n"); 
    131.     } 
    132.      
    133.     matrix_mul_cplx(x,r1,c1,y,c2,r); 
    134.     printf("x*y = : \n"); 
    135.     for(i=0;i<r1;++i) 
    136.     { 
    137.         for(j=0;j<c2;++j) 
    138.         { 
    139.             printf("%lf + j*%lf\t",r[i*c2 + j].real , r[i*c2 + j].image); 
    140.         } 
    141.         printf("\n"); 
    142.     } 
    143.     */  
    144.     //DSPF_dp_mat_mul_cplx test  
    145.     Complex x[4] = {{1,2}, {3,4}, {5,6}, {7,8}};    //矩阵按行存放  
    146.     int r1 = 2, c1 = 2,c2 = 1;  
    147.     Complex y[2] = {{1,2}, {3,4}};  
    148.     Complex r[2] = {{0,0}, {0,0}};  
    149.     int i,j;  
    150.     //printf("r  :  %lf\n%lf\n",r[0],r[1]);  
    151.       
    152.     printf("r = : \n");  
    153.     for(i=0;i<r1;++i)  
    154.     {  
    155.         for(j=0;j<c2;++j)  
    156.         {  
    157.             printf("%lf + j*%lf\t",r[i*c2 + j].real , r[i*c2 + j].image);  
    158.         }  
    159.         printf("\n");  
    160.     }  
    161.       
    162.     DSPF_dp_mat_mul_cplx((double*)x, r1,c1, (double*)y, c2,  (double*)r);  
    163.   
    164.     printf("x*y = : \n");  
    165.     for(i=0;i<r1;++i)  
    166.     {  
    167.         for(j=0;j<c2;++j)  
    168.         {  
    169.             printf("%lf + j*%lf\t",r[i*c2 + j].real , r[i*c2 + j].image);  
    170.         }  
    171.         printf("\n");  
    172.     }  
    173.   
    174.     return 0;  
    175. }  
    176. </span></pre><pre name="code" class="cpp"></pre><span style="font-size:14px"><br>  
    177. </span><p></p>  
    178. <pre></pre>  
    179. <span style="font-size:14px"><br>  
    180. </span>  
    181. <p></p>  
    182. <p><span style="font-size:14px"><br>  
    183. </span></p>  
    184. <p><br>  
    185. </p>  
    展开全文
  • 接着上一篇博客CBF中for循环变矩阵乘法的思想(arrayfire)的续。 上一篇主要讲了算法思想的改变,但是只是测试了实数,没有测试复数的效果,实际项目中都是复数的运用,所以这次添上复数的代码及测试结果。 这次在...

    接着上一篇博客CBF中for循环变矩阵乘法的思想(arrayfire)的续。

    上一篇主要讲了算法思想的改变,但是只是测试了实数,没有测试复数的效果,实际项目中都是复数的运用,所以这次添上复数的代码及测试结果。

    这次在添加arrayfire的代码之前,先看看不用这个库的一个C++代码形式:

    for(i=0;i<360;i++) //角度搜索 (-90:0.5:89)
    	{
    		theta=theta+0.5;
    		for(j=0;j<NN;j++)    //产生复指数信号exp(-jay*2*pi*[0:N-1]*f0*d/c*sin(theta(i)*pi/180));
    		{
    			a_s[1][j+1].real=cos(-2*PI*j*f0*d/c*sin(theta*PI/180)); 
    			a_s[1][j+1].imag=sin(-2*PI*j*f0*d/c*sin(theta*PI/180));
    		}
    		//cout<<"as"<<endl;
    		//CMatrix_print(a_s,1,3);
    		CMatrix_multiply(a_s,1,NN,R1,NN,NN,temp_11);
    		CMatrix_transpose(a_s,1,NN,a_sT);
    		CMatrix_multiply(temp_11,1,NN,a_sT,NN,1,temp_22);
    		//cout<<endl;
    		//CMatrix_print(temp_22,1, 1);
    		temp_G=Cabs(temp_22[1][1])*Cabs(temp_22[1][1]);  //G_MVDR(k)=abs((w_MVDR*s).^2);
    		vec1.push_back(temp_G);
    		*(G_beam+i)=10*(log(temp_G)/log(10));          // G=10.*log10(G);
    	}
    然后是arrayfire的代码形式

    /*******************************************************
     * Copyright (c) 2016.6.19
     * 算法的改进:主要解决了两个for循环的问题
     *第一个for是角度的问题,循环360个角度或者更小来得到每个角度的
     *输出功率。
     *第二个for循环是通道的问题,目前是8通道,对于arrarfire来讲这个很好做
     *直接生成1行8列的矩阵
     * 核心在于第一版本需要用for循环角度,现在的改变是将角度做成360行1列的矩阵
     *直接与通道矩阵相乘,变成360行8列的新矩阵a_s,此时注意求功率谱时候a_s*R*a_sT
     *变成了360*360的矩阵,此时只有对角线上的360个点有用,对应360个循环,虽然在矩阵乘法
     *上量有所增加,但是从根本上解决了并行问题,初步简单(没有其他常数和复数)
     *测试 速度比之前提升至少100倍,并且arrafire具有很好的稳定性,顺便推广一下。
     **此处只是思想,代码指示简单测试**
     ********************************************************/
    
    #include <arrayfire.h>
    #include <cstdio>
    #include <cstdlib>
    #include<vector>
    #include"iostream"
    
    
    using namespace af;
    using namespace std;
    std:: vector<float> vecl;
    int main(int argc, char *argv[])
    {
        try {
    
    
            // Select a device and display arrayfire info
            int device = argc > 1 ? atoi(argv[1]) : 0;
            af::setDevice(device);
            af::info();
    		/*常规波束形成原型*/
    		/*
    		for i = 1:length(theta)
             a_s = exp(-jay*2*pi*[0:N-1]*f0*d/c*sin(theta(i)*pi/180));
               beam_single(i) =(a_s *(R1) *a_s');
            end*/
    		/*arrayfire改进版*/
    		int const NN=8;//8通道
    		int const MM=32768;//每通道的数据长度
    		double const PI=3.1415;
    		double const d=0.05;//阵元间距
    		double const f0=15000.0;//频率
    		double const c=1500.0;//声速
                    array DATAReal = randu(MM, NN);
    		array DATAImag = randu(MM, NN);
    		array DATA=af::complex(DATAReal,DATAImag);
    		array DATAT=transpose(DATA,true);
    		array R=matmul(DATAT,DATA);
    		//af_print(R);
    		array theta=seq(0,359);
    		//af_print(theta);
    		array N=transpose(seq(1,8));
    		//af_print(N);
    		timer::start();
    		//array a_sReal=cos(-2*PI*N*f0*d/c*sin(theta*PI)/180);
    		array a_s_N=(-2*PI*N*f0*d/c);
    		//af_print(a_s_N);
    		array a_s_theta=(sin(theta*PI/180));
    		//af_print(a_s_theta);
    		array a_sO=matmul(a_s_theta,a_s_N);//其余是常数,省略处理,。
    		//af_print(a_sO);
    		array a_sReal=cos(a_sO);
    		array a_sImag=sin(a_sO);
    		array a_s=af::complex(a_sReal,a_sImag);
    		array a_sT=transpose(a_s,true);
    		array beamtemp=matmul(a_s,R);
    		//int s=DATA.dims(0);
    		//int s1=DATA.dims(1);
    		//af_print( s);
    		//af_print(s1);
    		array beamDOA=matmul(beamtemp,a_sT);
    		array beamDOA_N=af::abs(diag(beamDOA));
    		//af_print(beamDOA_N);
    		printf("version 2 elapsed seconds: %g\n", timer::stop());
    		
    
    		timer::start();
    		/*arrayfire基础版*/
    		int thetaN=360;
    		array beamDOAT=randu(360,1);
    		//af_print(beamDOAT);
    		for(int i=0;i<thetaN;i++)
    		{
    			array a_sRealN=cos(-2*PI*N*f0*d/c*sin(i*PI/180));
    		    array a_simagN=sin(-2*PI*N*f0*d/c*sin(i*PI/180));
    			array a_sN=af::complex(a_sRealN,a_simagN );
    			//af_print(a_sN);
    			array a_sNT=transpose(a_sN,true);
    			array temp1=matmul(a_sN,R);
    			//af_print(temp1);
    			array temp2=matmul(temp1,a_sNT);
    			//af_print(temp2);
    			array temp_G=af::abs(temp2);
    			beamDOAT(i,0)=temp_G(0);
    		}
    		//af_print(beamDOAT);
    
    		printf("version 1elapsed seconds: %g\n", timer::stop());
    
        } catch (af::exception& e) {
    
            fprintf(stderr, "%s\n", e.what());
            throw;
        }
    	system("pause");
        return 0;
    }
    
    实际的运行结果如下所示:




    同等数据量下对比CUDA程序运行结果


    同等数据量下对比C++程序运行结果为0.19s。


    综合数据来看,arrayfire在处理复数情况下速度会大幅降低(实数时是0.001s,见另外一篇博客测试),但是从效果来看,不会比CPU或者CUDA(GPU)慢,并且很稳定,代码写起来简单,所以值得一试。


    展开全文
  • CUDA之矩阵乘法——复数

    千次阅读 2016-05-30 16:55:25
    做好矩阵乘法和转置之后本来开心得不行的! 准备上手做个最基本的波束形成了! 突然发现希尔伯特变换完以后需要进行各种复数的运算…所以临时补写了一个复数乘法… 学着学着好像有点感觉了~!还是蛮有意思的。...

    做好矩阵乘法和转置之后本来开心得不行的!
    准备上手做个最基本的波束形成了!
    突然发现希尔伯特变换完以后需要进行各种复数的运算…所以临时补写了一个复数乘法…
    学着学着好像有点感觉了~!还是蛮有意思的。当然前提是能调试成功。
    用一句傅小姐的名言鼓励一下“只要心甘情愿任何事情都会变得简单!”。

    代码

    __device__ float GetReal(const Matrix A, int row, int col) {
    
        return A.real[row * A.stride + col];
    }
    
    __device__ float GetImag(const Matrix A, int row, int col) {
        return A.imag[row * A.stride + col];
    }
    
    __device__ void SetElement(Matrix A, int row, int col, float valueR, float valueI) {
        A.real[row * A.stride + col] = valueR;
        A.imag[row * A.stride + col] = valueI;
    }
    
    __device__ Matrix GetSubMatrix(Matrix A, int row, int col) {
        Matrix Asub;
        Asub.width = BLOCK_SIZE;
        Asub.height = BLOCK_SIZE;
        Asub.stride = A.stride;
        Asub.real = &A.real[A.stride * BLOCK_SIZE * row+ BLOCK_SIZE * col];
        Asub.imag = &A.imag[A.stride * BLOCK_SIZE * row+ BLOCK_SIZE * col];
        return Asub;
    }
    
    __global__ void CMatMulKernel(Matrix A, Matrix B, Matrix C) {
        int blockRow = blockIdx.y;
        int blockCol = blockIdx.x;
        Matrix Csub = GetSubMatrix(C, blockRow, blockCol);
        float CvalueR = 0;
        float CvalueI = 0;
        int row = threadIdx.y;
        int col = threadIdx.x;
        for (int m = 0; m < (A.width / BLOCK_SIZE); ++m) {
            Matrix Asub = GetSubMatrix(A, blockRow, m);
            Matrix Bsub = GetSubMatrix(B, m, blockCol);
            __shared__ float AsR[BLOCK_SIZE][BLOCK_SIZE];
            __shared__ float AsI[BLOCK_SIZE][BLOCK_SIZE];
            __shared__ float BsR[BLOCK_SIZE][BLOCK_SIZE];
            __shared__ float BsI[BLOCK_SIZE][BLOCK_SIZE];
            AsR[row][col] = GetReal(Asub, row, col);
            AsI[row][col] = GetImag(Asub, row, col);
            BsR[row][col] = GetReal(Bsub, row, col);
            BsI[row][col] = GetImag(Bsub, row, col);
            __syncthreads();
            for (int e = 0; e < BLOCK_SIZE; ++e)
            {
                CvalueR += AsR[row][e] * BsR[e][col]-AsI[row][e]*BsI[e][col];
                CvalueI += AsR[row][e] * BsI[e][col]+AsI[row][e]*BsR[e][col];
            }
            __syncthreads();
        }
        SetElement(Csub, row, col, CvalueR,CvalueI);
    }
    
    void CMatMul(const Matrix A, const Matrix B, Matrix C) {
        Matrix d_A;
        d_A.width = d_A.stride = A.width; 
        d_A.height = A.height;
        size_t size = A.width * A.height * sizeof(float);
        cudaMalloc((void**)&d_A.real, size);
        cudaMalloc((void**)&d_A.imag, size);
        cudaMemcpy(d_A.real, A.real, size,
        cudaMemcpyHostToDevice);
        cudaMemcpy(d_A.imag, A.imag, size,
        cudaMemcpyHostToDevice);
    
        Matrix d_B;
        d_B.width = d_B.stride = B.width; 
        d_B.height = B.height;
        size = B.width * B.height * sizeof(float);
        cudaMalloc((void**)&d_B.real, size);
        cudaMalloc((void**)&d_B.imag, size);
        cudaMemcpy(d_B.real, B.real, size,
        cudaMemcpyHostToDevice);
        cudaMemcpy(d_B.imag, B.imag, size,
        cudaMemcpyHostToDevice);
    
        Matrix d_C;
        d_C.width = d_C.stride = C.width; 
        d_C.height = C.height;
        size = C.width * C.height * sizeof(float);
        cudaMalloc((void**)&d_C.real, size);
        cudaMalloc((void**)&d_C.imag, size);
    
        dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
        dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);
        CMatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);
        cudaMemcpy(C.real, d_C.real, size,
        cudaMemcpyDeviceToHost);
        cudaMemcpy(C.imag, d_C.imag, size,
        cudaMemcpyDeviceToHost);
    
        cudaFree(d_A.real);
        cudaFree(d_A.imag);
        cudaFree(d_B.real);
        cudaFree(d_B.imag);
        cudaFree(d_C.real);
        cudaFree(d_C.imag);
    }

    为什么这么久过去了我依然不能上传图片。

    展开全文
  • Intel MKL 求解大型稀疏复数矩阵 (C/C++)实例。 适用与求解大型对称或者非对称稀疏复数矩阵
  • 复数 矩阵.rar

    2019-07-15 10:59:47
    (1)完成矩阵加法,乘法运算. (2)完成复数的结构定义,并实现复数的赋值,加、减、输出等操作。
  • 矩阵乘法的一些应用

    万次阅读 多人点赞 2017-09-18 20:35:35
    矩阵乘法 例题 例题1斐波那契数列 problem 构造矩阵式 矩阵快速幂code 矩阵乘法code code 例题2JZOJsenior1240Fibonacci sequence problem analysis code 例题3NOI2012随机数生成器 problem 构造矩阵式 快速
  • 矩阵乘法与运用

    2017-08-14 20:31:54
    矩阵乘法:矩阵相乘最重要的方法是一般矩阵乘积。它只有在第一个矩阵的列数(column)和第二个矩阵的行数(row)相同时才有意义[1] 。一般单指矩阵乘积时,指的便是  一般矩阵乘积。一个m×n的矩阵就是m×n个数排成...
  • 文件为PYNQ-Z2板实现矩阵乘法加速所需文件,详细操作流程可见博客:https://blog.csdn.net/qq_42334072/article/details/106769534
  • C++实现复数矩阵求逆 matlab inv

    千次阅读 2020-07-06 11:10:57
    C++实现复数矩阵求逆 matlab inv一、引言二、原理2.1 实数矩阵求逆2.2 复数矩阵求逆三、代码四、测试 一、引言 之前偶尔一次有用到将matlab转为C++语言的需求,其中matlab有一个inv函数可以非常方便的求矩阵的逆,...
  • 实现复数矩阵的加减乘除

    千次阅读 2019-07-15 23:03:44
    (2)使用*实现两个复数矩阵乘法; (3)使用!实现复数矩阵的共轭转置操作; (4)使用<<与>>实现复数矩阵的输出与输入; (5)要求兼容实数矩阵运算 其实:;兼容矩阵运算也可以使用转换构造函数,...
  • 支持任意阶的复数矩阵求逆,a是复数矩阵的实属部分,b是虚数部分,c、d分别是输出实属和虚数部分
  • 矩阵乘法是线性代数中常见的问题之一, 许多数值计算问题都包含着矩阵乘法的计算
  • 复数乘法法则 设 x=a+bi,y=c+di (a、b、c、d∈R)是任意两个复数,那么它们的积(a+bi)(c+di)=(ac-bd)+(bc+ad)i.  其实就是把两个复数相乘,类似两个多项式相乘,展开得: ac+adi+bci+bdi^2,  由i^2=-1,所以...
  • 矩阵乘法 稀疏矩阵与向量相乘(向量中的所有初始项均为1),结果存储在同一向量中,以进行下一次迭代。 迭代次数应指定为参数参数。 具有稀疏矩阵的.txt文件也应作为参数提供 如果要查看结果,则应将第三个参数...
  • MKL库矩阵乘法

    千次阅读 2016-05-21 16:13:22
    MKL库矩阵乘法 Posted on 2008年09月3日by cici1020 我再次来学术一把。 我今天因为不想看3000多页的大文档,就在GOOGLE狂翻用Intel的Math Kernel Library的库如何作矩阵运算。最后还是看那个3000多页...
  • 复数矩阵的加/减/乘运算涉及到其复数元胞(cell)的相加减运算,由于complex.h头文件中只给出了复数乘法运算,故而复数的加减运算函数需要自己定义功能。 二、矩阵加/减法运算 即复数矩阵复数矩阵之间的加/减法...
  • 矩阵乘法回顾

    2020-09-29 20:54:56
    常用矩阵加速,和其它知识点结合
  • 复数矩阵求逆的 C 语言程序

    千次阅读 2019-05-20 11:38:28
    关于复数矩阵求逆,如果使用 MATLAB,就非常简单。我们先用一个 MATLAB 的例子来说明,等会要将 C 语言的程序和 MATLAB 的程序进行对比。 close all; clear all; clc; %定义矩阵a为复数矩阵 a = [[4+2*i,3+1*i,4+3*...
  • 复数矩阵相乘的扩展矩阵计算方法

    千次阅读 2018-07-04 21:44:18
    在诸多信号处理领域,我们所涉及到的信号都是复数信号,那么matlab如何处理复数矩阵呢? 假设矩阵AA\boldsymbol{A}和BB\boldsymbol{B}为复数矩阵,当我们运用matlab计算两个矩阵的乘积时(矩阵AA\boldsymbol{A}和BB\...
  • 复数矩阵行列式VB和C#计算程序 这是一个计算复数矩阵行列的VB和C#代码,需要的下载来看看哦
  • 矩阵-矩阵乘法 复数矩阵 矩阵符号 如果用表示所有实数的集合,那么我们用表示所有的实数矩阵组成的向量空间,即: 其中,大写字母(如)表示矩阵,带下标的小写字母(如)表示矩阵中的元素。除了用表示矩阵中第...
  • 二维矩阵乘法加法

    2020-03-20 20:08:58
    二维矩阵乘法和加法 引用方法 package 乘法; public class Jvzhencf { int data[][]; private int hang; private int lie; public int getHang() { return hang; } public void setHang(int hang) { this.hang = ...
  • 我再次来学术一把。我今天因为不想看3000多页...前人栽树后人乘凉……关键词:Intel MKL MKL库 矩阵乘积 矩阵乘法 样例 例子 代码 标程 sample code这样以后人家都能搜到这里了:DMKL的大文档里,有各种matrix*matri...
  • cublas 矩阵乘法

    2019-12-11 11:13:10
    cublas是cuda用来解决线性代数的问题的一个函数库,而且对于矩阵运算来说,其效率比大部分人自己写核函数高不少,只是cublas不同于C++,是列优先存储,因此参数一不小心设的不对,结果大不相同,所以记录下来;...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 6,858
精华内容 2,743
关键字:

复数矩阵乘法