复数矩阵计算行列式

来源:互联网 发布:flash游戏制作软件 编辑:程序博客网 时间:2024/04/30 03:03
项目上需要对复矩阵的行列式计算,根据计算一般矩阵行列式的代码改成了复矩阵行列式计算。
#include<iostream>#include<stdio.h>#include<stdlib.h>#include<time.h>#include<fstream>#include <iomanip>using namespace std;#define ROW 25#define COL 25typedef struct  {float Real;float Image;}Complex;Complex add(Complex a,Complex b){Complex c;c.Real=a.Real+b.Real;c.Image=a.Image+b.Image;return c;}Complex sub(Complex a,Complex b){Complex c;c.Real=a.Real-b.Real;c.Image=a.Image-b.Image;return c;}Complex Mul(Complex a,Complex b){Complex c;c.Real=a.Real*b.Real-a.Image*b.Image;c.Image=a.Real*b.Image+b.Real*a.Image;return c;}Complex Div(Complex a,Complex b){Complex c;c.Real=(a.Real*b.Real+a.Image*b.Image)/(b.Real*b.Real+b.Image*b.Image);c.Image=(a.Image*b.Real-a.Real*b.Image)/(b.Real*b.Real+b.Image*b.Image);return c;}Complex matrix_det(const Complex *mat, int n){    int i,j,k,is,js,l,v;Complex det;float flag,pivot,tmp;    Complex *cpmat;Complex tempValue;    if(mat == NULL)                            /* 检查输入的指针是否为空*/    {      printf("matrix pointer is NULL.\n");    }    cpmat = (Complex*)malloc(n*n*sizeof(Complex));  /* 将输入矩阵的内容拷贝一份,以免破坏*/    for(i=0; i<n*n; i++)      cpmat[i] = mat[i];det.Real = 1.0;                                 /* 设置行列式值初置*/ det.Image=0.0;flag= 1.0;                                /* 这个变量原来记录行列式值的符号*/    for(k=0; k<n-1; k++){           /* 最多进行n-1次消去*/         pivot = 0.0;                             /* 选择主元*/for(i=k; i<n; i++){for(j=k; j<n; j++){tmp = cpmat[i*n+j].Real*cpmat[i*n+j].Real+cpmat[i*n+j].Image*cpmat[i*n+j].Image;if(tmp > pivot){pivot = tmp;is = i;js = j;}}}if(pivot < 1e-5){                          /* 如果找到的主元小于eps,则认为是0。*/det.Real = 0.0;                             /*此时行列式值也是0。*/det.Image=0.0;return(det);} if(is != k){                              /* 判断是否需要行交换*/flag = -flag;                          /* 行交换一次,行列式值变号*/for(j=k; j<n; j++){                     /* 进行行交换*/l = k*n + j;v = is*n + j;tempValue = cpmat[l];cpmat[l] = cpmat[v];cpmat[v] = tempValue;}      }      if(js != k){                              /* 判断是否需要列交换*/flag = -flag;                          /* 列交换一次,行列式值变号*/for(i=k; i<n; i++){                     /* 进行列交换*/l = i*n + k;v = i*n + js;tempValue = cpmat[v];cpmat[v] = cpmat[l];cpmat[l] = tempValue;}     }     for(i=k+1; i<n; i++){                    /* 进行消去*/tempValue=Div(cpmat[i*n+k],cpmat[k*n+k]);   /* 记录下此值,减少除法的次数*/   for(j=k+1; j<n; j++) {                 /* 消去*/Complex aa;aa=Mul(tempValue,cpmat[k*n+j]);cpmat[i*n+j]=sub(cpmat[i*n+j],aa);} }      det = Mul(det,cpmat[k*n+k]);           /*更新det的值*/    }det = Mul(det,cpmat[k*n+k]); /* 最终更新det的值*/det.Real=flag*det.Real;det.Image=flag*det.Image;    free(cpmat);           return(det);}int main(){Complex *inputdata;Complex result;inputdata=(Complex*)malloc(sizeof(Complex)*ROW*COL);ofstream myfile("example.txt");srand(time(0));for(int i=0;i<ROW*COL;i++){inputdata[i].Real=(float)rand()/10000;inputdata[i].Image=(float)rand()/10000;}//cout<<"[";myfile<<"b=[";for(int i=0;i<ROW;i++){for(int j=0;j<COL;j++){//cout<<inputdata[i*COL+j].Real<<"+"<<inputdata[i*COL+j].Image<<"i";myfile<<inputdata[i*COL+j].Real<<"+"<<inputdata[i*COL+j].Image<<"i";if(j<COL-1){//cout<<",";myfile<<",";}}if(i<COL-1){//cout<<";";myfile<<";";}}//cout<<"]"<<endl;myfile<<"]"<<endl;;result=matrix_det(inputdata,ROW);//cout<<"Result:"<<endl;//cout<<result.Real<<"+"<<result.Image<<"i"<<endl;myfile<<setiosflags(ios::showpoint)<<result.Real<<"+"<<result.Image<<"i"<<endl;myfile<<flush;myfile.close();return 0;}
计算结果与matlab计算对比,没有问题
0 0
原创粉丝点击