C语言中复数运算及调用blas,lapack中复数函数进行科学计算

C语言中常用的数据类型主要int, float ,double ,char 等，但在科学运算中复数扮演着重要角色。这里讲下C语言中的复数运算以及如何调用blas，lapack库中的复数函数来进行科学计算．

１．C语言中的复数运算。

C语言中若要用的复数，需要包含头文件complex.h，下面看看一些基本的例子

#include <stdio.h>
#include"complex.h"

int main()
{
    complex a, b, c, d, f;
    a = 1 + 2*I;
    b = 3 + 4*I;
    c = a+b;
    d = a*b;
    f = a/b;
    printf("complex a plus b is :%f, %f\n", creal(c), cimag(c));
    printf("complex a multiply b is :%f, %f\n", creal(d), cimag(d));
    printf("complex a divided b is :%f, %f\n", creal(f), cimag(f));

}

说明：

     (1).creal(c), cimag(c)表示分别取复数c的实部和虚部；

     (2).complex 默认是complex double 型数据结构，即实部和虚部都是double 型，当然也可以定义complex int型；

　 (3).虚数单位是用I来表示

生成可执行文件：

gcc -o test_complex test_complex1

得到的结果如下：

complex a plus b is :4.000000, 6.000000
complex a multiply b is :-5.000000, 10.000000
complex a / b is :0.440000, 0.080000

２．下面讲下如何利用调用blas的复函数来进行科学运算。

在blas，lapack中有三种表示复数的形式：利用C语言中的complex数据类型、结构体以及实部和复部交替存储。下面举例说明：

(1).利用C语言中complex数据类型

#include <stdio.h>
#include"stdlib.h"
#include"complex.h"

complex zdotu_(int *, complex *, int *,complex *, int *);    

int main()
{
    complex *bb, *cc, al;

    int i, j,  N, ONE=1;
    N =2;

    bb = malloc(N *sizeof(complex));
    cc = malloc(N *sizeof(complex));
        bb[0] = 1+I; bb[1]=1+I;
        cc[0] = 1+I; cc[1]=1+I;
        al= zdotu_(&N, cc, &ONE, bb, &ONE);
    printf("al:%f, %f\n", creal(al), cimag(al));

    free(bb);
    free(cc);
    return 0;
}

编译并生成可执行文件

gcc -o test_zdot.o -c test_zdot.c
gcc -o test_zdot test_zdot.c -lblas -lgfortran

执行：

./test_zdot

得到结果：

al:0.000000, 4.000000

(2).利用blas中定义的结构体CBLAS_TEST_ZOMPLEX

#include <stdio.h>
#include"stdlib.h"
#include"cblas_test.h"

#define zomplex CBLAS_TEST_ZOMPLEX 

zomplex zdotu_(int *, zomplex *, int *,zomplex *, int *);    
int main()
{

    zomplex *b, *c , nal;
    int i, N, ONE=1;
    N =2;
    b = malloc(N *sizeof(zomplex));
    c = malloc(N *sizeof(zomplex));

    for(i=0; i< N; i++)
    {
        (*(b+i)).real = 1;
        (*(b+i)).imag = 1;
        (*(c+i)).real = 1;
        (*(c+i)).imag = 1;
    }

    　nal= zdotu_(&N, b, &ONE, c, &ONE);
    printf("nal:%f, %f\n", nal.real, nal.imag);

    free(b);
    free(c);
    return 0;
}

说明：a. blas中在头文件cblas_test.h定义了结构体CBLAS_TEST_ZOMPLEX，其成员变量分别为：real, imag.

           b.由于结构体CBLAS_TEST_ZOMPLEX的书写太长，例子中通过#define zomplex CBLAS_TEST_ZOMPLEX ，用zomplex来取代。

编译和(1)中一样，链接好blas, gfortran库即可，得到结果如下：

nal:0.000000, 4.000000

(３).利用实部和虚部交叉存储

在blas中定义了结构体，其在物理存储空间中，其实部和虚部都是交叉存储的，即存第一个数的实部，再存第一个数的虚部，第二个数的实部，第二个数的虚部...,那么也可以定义数组，交替存储实部和虚部。例子如下：

#include <stdio.h>
#include"stdlib.h"
#include"cblas_test.h"

#define zomplex CBLAS_TEST_ZOMPLEX 
zomplex zdotu_(int *, double *, int *, double *, int *);    
int main()
{

    zomplex al;
   　double x[4], y[4];
    int i, j, N, ONE=1;
    N =2;

    x[0] = 1.0; x[1] = 1.0; x[2] = 1.0; x[3] = 1.0;
    y[0] = 1.0; y[1] = 1.0; y[2] = 1.0; y[3] = 1.0;

    　al= zdotu_(&N, x,  &ONE, y, &ONE);
    printf("al:%f, %f\n", al.real, al.imag);

    return 0;
}

注意：函数声明如下:

zomplex zdotu_(int *, double *, int *, double *, int *); 

这里函数的形参数据类型为：double，而不是zomplex型了！！！

编译并执行，得到得到结果如下：

nal:0.000000, 4.000000

　　根据自己的选用那种方式就需要声明哪种函数，因此在调用blas编程，需要先声明函数，以免出现意想不到的结果（从这三个例子中也可以看出，函数名一样，如果不声明函数的形参的数据类型，那得到的结果可能会有问题）。

   调用lapack和blas类似，复数运算也是这三种方式，根据自己的喜好选择。