针对大型稀疏矩阵的求解，需要提前对矩阵进行预处理

/** * 预处理文件: data.txt 第1: int nSize(矩阵中非零个数) *N个int值表示每行中第一个数值的偏移 nSize个列 nSize个double类型的值 *v1_b.txt 下面N个double类型的值 *v1_x0.txt下面N个double类型的值 * **/ #include <stdio.h> #include <stdlib.h> #include <math.h> #define NX 360 #define NY 180 #define NZ 38 #define N NX*NY*NZ double *AV,*b,*x0,*MV; int *AI,*AJ; int nSize; int InitAll() { FILE *fp1=fopen("data.txt","rb"); if(fp1==NULL){ printf("can't open read file\n"); exit(0); } int size; fread(&size,4,1,fp1); AV=(double *)malloc(size* sizeof(double)); AI=(int*)malloc(N*sizeof(int)); AJ=(int*)malloc(size* sizeof(int)); fread(AI,4,N,fp1); fread(AJ,4,size,fp1); fread(AV,8,size,fp1); fclose(fp1); FILE *fp2=fopen("v1_b.txt","rb"); if(fp2==NULL){ printf("can't open read file\n"); exit(0); } b=(double *)malloc(N* sizeof(double)); fread(b,8,N,fp2); fclose(fp2); FILE *fp3=fopen("v1_x0.txt","rb"); if(fp3==NULL){ printf("can't open read file\n"); exit(0); } x0=(double *)malloc(N* sizeof(double)); fread(x0,8,N,fp3); fclose(fp3); return size; } void getAmVector(double *result,double *data) //矩阵向量相乘 { int ind=0; int i,j; int size; double tmp; for(i=0;i<N-1;i++) //N-1 { ind=AI[i]; size=AI[i+1]-AI[i]; tmp=0.0; for(j=0;j<size;j++){ tmp=tmp+AV[ind+j]*data[AJ[ind+j]]; } result[i]=tmp; } size=nSize-AI[N-1]; tmp=0.0; ind=AI[N-1]; for(j=0;j<size;j++){ tmp=tmp+AV[ind+j]*data[AJ[ind+j]]; } result[N-1]=tmp; } double checkSum(double *data) { int ind=0; int i,j; int size; double tmp; double sum=0.0; for(i=0;i<N-1;i++) //N-1行 { ind=AI[i]; size=AI[i+1]-AI[i]; tmp=0.0; for(j=0;j<size;j++){ tmp=tmp+AV[ind+j]*data[AJ[ind+j]]; } sum=sum+(tmp-b[i])*(tmp-b[i]); } size=nSize-AI[N-1]; //printf("%d\n",AI[N-1]); tmp=0.0; ind=AI[N-1]; for(j=0;j<size;j++){ tmp=tmp+AV[ind+j]*data[AJ[ind+j]]; } sum=sum+(tmp-b[N-1])*(tmp-b[N-1]); return sum; } double VectorDianJi(double *dataA,double *dataB) { int i; double sum=0.0; for(i=0;i<N;i++){ sum=sum+dataA[i]*dataB[i]; } return sum; } void JacobiCal(double *Rm,double *R) //矩阵Rm=M^(-1)*R Jacobi迭代 { int ind=0; int i,j,size; double tmp,m_rst; for(i=0;i<N-1;i++){ ind=AI[i]; tmp=0.0; m_rst=0.0; size=AI[i+1]-AI[i]; for(j=0;j<size;j++){ if(AJ[ind+j]==i){ //对角线上的元素 tmp=AV[ind+j]; }else{ m_rst=m_rst+AV[ind+j]*R[AJ[ind+j]]; } } Rm[i]= (R[i]-m_rst)/tmp; } size=nSize-AI[N-1]; ind=AI[N-1]; tmp=0.0; m_rst=0.0; for(j=0;j<size;j++){ if(AJ[ind+j]==(N-1)){ //对角线上的元素 tmp=AV[ind+j]; }else{ m_rst=m_rst+AV[ind+j]*R[AJ[ind+j]]; } } Rm[N-1]=(R[N-1]-m_rst)/tmp; } void calGCR() { int k=5; int i,j; double *R=(double*)malloc(N*sizeof(double)); double *Rm=(double*)malloc(N*sizeof(double)); double *p=(double*)malloc(N*sizeof(double)); double *Ap=(double*)malloc(N* sizeof(double)); double *ARm=(double*)malloc(N* sizeof(double)); double *Apj=(double*)malloc(N* sizeof(double)); double *pj[5]; double betaij[5]; for(i=0;i<5;i++){ pj[i]=(double*)malloc(N*sizeof(double)); betaij[i]=0; } getAmVector(R,x0); for(i=0;i<N;i++){ //get the R=b-Ax0 R[i]=b[i]-R[i]; } JacobiCal(Rm,R); //预处理 for(j=0;j<N;j++){ p[j]=Rm[j]; } double ai,checkVal; int ind=0,m=0; i=1; do{ printf("Jacobi times------>%d\n",i); getAmVector(Ap,p); //get Api-1 向量 for(j=0;j<N;j++){ pj[(i-1)%k][j]=p[j]; //存储5个子空间向量 } ai=VectorDianJi(R,Ap)/VectorDianJi(Ap,Ap); //printf("ai---->%.30lf\n",ai); for(j=0;j<N;j++){ //xi=x(i-1)+a(i-1)*p(i-1) x0[j]=x0[j]+ai*p[j]; R[j]=R[j]-ai*Ap[j]; } JacobiCal(Rm,R); //预处理 ind=0; for(j=((i-1)/k)*k;j<i;j++){ getAmVector(ARm,Rm); getAmVector(Apj,pj[ind]); betaij[ind]=-VectorDianJi(ARm,Apj)/VectorDianJi(Apj,Apj); ind++; } for(j=0;j<N;j++){ p[j]=Rm[j]; } ind=0; for(j=((i-1)/k)*k;j<i;j++){ for(m=0;m<N;m++){ p[m]=p[m]+betaij[ind]*pj[ind][m]; } ind++; } checkVal=checkSum(x0); printf("checkSum---->%.30lf\n",checkVal); i++; }while(checkVal>0.0000000001); free(R);free(Rm);free(p);free(Ap);free(ARm); free(Apj); for(i=0;i<5;i++){ free(pj[i]); } } double checkData() { int ind=0; int i,j; int size; double tmp; double sum=0.0; for(i=0;i<N-1;i++) //N-1行 { ind=AI[i]; size=AI[i+1]-AI[i]; tmp=0.0; for(j=0;j<size;j++){ tmp=tmp+AV[ind+j]*x0[AJ[ind+j]]; } sum=sum+(tmp-b[i])*(tmp-b[i]); } size=nSize-AI[N-1]; //printf("%d\n",AI[N-1]); tmp=0.0; ind=AI[N-1]; for(j=0;j<size;j++){ tmp=tmp+AV[ind+j]*x0[AJ[ind+j]]; } sum=sum+(tmp-b[N-1])*(tmp-b[N-1]); return sum; } int main(int argc, char **argv) { AI=NULL; AJ=NULL; AV=NULL; b=NULL; x0=NULL; printf(" start Init all data\n"); nSize=InitAll(); printf(" end Init all data----N: %d nSize: %d\n",N,nSize); // double checkValue=checkData(); // printf("checksum ---->%.30lf\n",checkValue); calGCR(); free(AI); free(AJ);free(AV);free(b);free(x0); free(MV); return 0; }