复数的乘积c语言,复数矩阵乘法C语言实现

下面代码参考TI的实现：/* NAME */

/* DSPF_dp_mat_mul_cplx -- Complex matrix multiplication */

/* */

/* USAGE */

/* */

/* This routine has the following C prototype: */

/* */

/* void DSPF_dp_mat_mul_cplx( */

/* const double* x, */

/* int r1, */

/* int c1, */

/* const double* y, */

/* int c2, */

/* double* restrict r */

/* ) */

/* */

/* x[2*r1*c1]: Input matrix containing r1*c1 complex */

/* floating point numbers having r1 rows and c1 */

/* columns of complex numbers. */

/* r1 : No. of rows in matrix x. */

/* c1 : No. of columns in matrix x. */

/* Also no. of rows in matrix y. */

/* y[2*c1*c2]: Input matrix containing c1*c2 complex */

/* floating point numbers having c1 rows and c2 */

/* columns of complex numbers. */

/* c2 : No. of columns in matrix y. */

/* r[2*r1*c2]: Output matrix of c1*c2 complex floating */

/* point numbers having c1 rows and c2 columns of */

/* complex numbers. */

/* */

/* Complex numbers are stored consecutively with */

/* real values stored in even positions and */

/* imaginary values in odd positions. */

/* */

/* DESCRIPTION */

/* */

/* This function computes the expression "r = x * y" for the */

/* matrices x and y. The columnar dimension of x must match the row */

/* dimension of y. The resulting matrix has the same number of rows */

/* as x and the same number of columns as y. */

/* */

/* Each element of the matrix is assumed to be complex numbers with */

/* Real values are stored in even positions and imaginary */

/* values in odd positions. */

/* */

/* TECHNIQUES */

/* */

/* 1. Innermost loop is unrolled twice. */

/* 2. Outermost loop is executed in parallel with innner loops. */

/* */

/* ASSUMPTIONS */

/* */

/* 1. r1,r2>=1,c1 should be a multiple of 2 and >=2. */

/* 2. x should be padded with 6 words */

/* */

/* */

/* C CODE */

/* */

void DSPF_dp_mat_mul_cplx(const double* x, int r1, int c1,

const double* y, int c2, double* r)

{

double real, imag;

int i, j, k;

for(i = 0; i < r1; i++)

{

for(j = 0; j < c2; j++)

{

real=0;

imag=0;

for(k = 0; k < c1; 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. Real values are stored in even word positions and imaginary */

/* values in odd positions. */

/* 2. Endian: This code is LITTLE ENDIAN. */

/* 3. Interruptibility: This code is interrupt tolerant but not */

/* interruptible. */

/* */

/* CYCLES */

/* */

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

/* When r1=3, c1=4, c2=4, cycles = 640 */

/* When r1=4, c1=4, c2=5, cycles = 1040 */

/* */

/* CODESIZE */

/* */

/* 832 bytes */

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

/* Copyright (c) 2004 Texas Instruments, Incorporated. */

/* All Rights Reserved. */

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

#ifndef DSPF_DP_MAT_MUL_CPLX_H_

#define DSPF_DP_MAT_MUL_CPLX_H_ 1

void DSPF_dp_mat_mul_cplx(

const double* x,

int r1,

int c1,

const double* y,

int c2,

double* restrict r

);

#endif

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

/* End of file: DSPF_dp_mat_mul_cplx.h */

/* ------------------------------------------------------------------------ */测试程序如下：(说明：下面代码分别包含了三个矩阵相乘的函数)#include #include

#include "complex.h"

void matrix_mul(

const double* x,

int r1,

int c1,

const double* y,

int c2,

double* r

)

{

int i,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];

}

}

}

}

void matrix_mul_cplx(

const Complex* x,

int r1,

int c1,

const Complex* y,

int c2,

Complex* r

)

{

int i,j,k;

Complex temp;

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);

}

}

}

}

void DSPF_dp_mat_mul_cplx(const double* x, int r1, int c1,

const double* y, int c2, double* r)

{

double real, imag;

int i, j, k;

for(i = 0; i < r1; i++)

{

for(j = 0; j < c2; j++)

{

real=0;

imag=0;

for(k = 0; k < c1; 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;

}

}

}

int main()

{

//matrix_mul test

/*

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

int r1 = 3, c1 = 3,c2 = 1;

double y[3] = {10,11,12};

double r[3] = {0};

int i,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

{

for(j=0;j

{

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

}

printf("\n");

}

*/

/*

//matrix_mul_cplx test

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

int r1 = 2, c1 = 2,c2 = 1;

Complex y[2] = {{1,2}, {3,4}};

Complex r[2] = {{0,0}, {0,0}};

int i,j;

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

printf("r = : \n");

for(i=0;i

{

for(j=0;j

{

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

}

printf("\n");

}

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

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

for(i=0;i

{

for(j=0;j

{

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

}

printf("\n");

}

*/

//DSPF_dp_mat_mul_cplx test

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

int r1 = 2, c1 = 2,c2 = 1;

Complex y[2] = {{1,2}, {3,4}};

Complex r[2] = {{0,0}, {0,0}};

int i,j;

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

printf("r = : \n");

for(i=0;i

{

for(j=0;j

{

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

}

printf("\n");

}

DSPF_dp_mat_mul_cplx((double*)x, r1,c1, (double*)y, c2, (double*)r);

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

for(i=0;i

{

for(j=0;j

{

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

}

printf("\n");

}

return 0;

}