qr.rar_openmp_openmp qr_qr_qr openmp_qr分解

#include "omp.h" #include "stdio.h" #include "stdlib.h" #include <sys/time.h> #define a(x,y) a[x*M+y] #define q(x,y) q[x*M+y] #define A(x,y) A[x*M+y] #define Q(x,y) Q[x*M+y] #define R(x,y) R[x*M+y] #define Mark(x,y) Mark[x*M+y] int bcomplite; float temp; float *A; float *R; float *Q; int *Mark; double starttime; double time1; double time2; int p; double omp_time() { static int sec = -1; struct timeval tv; gettimeofday(&tv, (void *)0); if (sec < 0) sec = tv.tv_sec; return (tv.tv_sec - sec) + 1.0e-6*tv.tv_usec; } int main() { int M,N,m; int i,j,k,myid; float *ai,*qi,*aj,*qj; float c,s,sp; float *a,*q; FILE *fdA; p=1; omp_set_num_threads(p); starttime = omp_time(); printf("time1=%f\n",starttime); fdA=fopen("dataIn.txt","r"); fscanf(fdA,"%d %d", &M, &N); if(M != N) { puts("The input is error!"); exit(0); } A=(float*)malloc(sizeof(float)*M*M); Q=(float*)malloc(sizeof(float)*M*M); R=(float*)malloc(sizeof(float)*M*M); for(i = 0; i < M; i ++) { for(j = 0; j < M; j ++) fscanf(fdA, "%f", A+i*M+j); } fclose(fdA); #pragma omp parallel for for(i=0;i<M;i++) { #pragma omp parallel for for(j=0;j<M;j++) if (i==j) Q(i,j)=1.0; else Q(i,j)=0.0; } m=M/p; if (M%p!=0) m++; Mark=(int*)malloc(sizeof(int)*p*M);//在p线程上的整体i行是否执行 qi=(float*)malloc(sizeof(float)*M); qj=(float*)malloc(sizeof(float)*M); aj=(float*)malloc(sizeof(float)*M); ai=(float*)malloc(sizeof(float)*M); a=(float*)malloc(sizeof(float)*M*M); q=(float*)malloc(sizeof(float)*M*M); #pragma omp parallel for for(i=0;i<M;i++) { #pragma omp parallel for for(j=0;j<M;j++) { a(i,j)=A(i,j); q(i,j)=Q(i,j); } } #pragma omp parallel for for(i=0;i<p;i++) { #pragma omp parallel for for(j=0;j<M;j++) { Mark(i,j)=0; } } time1 = omp_time(); #pragma omp parallel default(shared) shared(Mark,m) private(bcomplite,i,j,k,c,s,sp,myid,temp) firstprivate(ai,aj,qi,qj) { myid=omp_get_thread_num(); if (myid==0) { for(j=0;j<m-1;j++) { for(i=j+1;i<m;i++) { sp=sqrt(a(j,j)*a(j,j)+a(i,j)*a(i,j)); c=a(j,j)/sp; s=a(i,j)/sp; for(k=0;k<M;k++) { aj[k]=c*a(j,k)+s*a(i,k); qj[k]=c*q(j,k)+s*q(i,k); ai[k]=-s*a(j,k)+c*a(i,k); qi[k]=-s*q(j,k)+c*q(i,k); } for(k=0;k<M;k++) { a(j,k)=aj[k]; q(j,k)=qj[k]; a(i,k)=ai[k]; q(i,k)=qi[k]; } } /* i */ Mark(myid,j)=1;//myid线程上的i行已经处理完成，下一个线程可以处理i行 } Mark(myid,m-1)=1; /* for j */ } /* my_rank==0 */ else /* my_rank!=0 */ { if (myid!=(p-1)) { // for(k=0;k<myid-1;k++)//上一线程完成整体的某行 for(j=0;j<(myid*m);j++) { //判断上一线程完成整体的某行 bcomplite=0; while(bcomplite==0) { if(Mark((myid-1),j)==1) bcomplite=1; } for(i=(myid*m);i<(myid*m+m);i++) { sp=sqrt(a(j,j)*a(j,j)+a(i,j)*a(i,j)); c=a(j,j)/sp; s=a(i,j)/sp; for(k=0;k<M;k++) { aj[k]=c*a(j,k)+s*a(i,k); qj[k]=c*q(j,k)+s*q(i,k); ai[k]=-s*a(j,k)+c*a(i,k); qi[k]=-s*q(j,k)+c*q(i,k); } for(k=0;k<M;k++) { a(j,k)=aj[k]; q(j,k)=qj[k]; a(i,k)=ai[k]; q(i,k)=qi[k]; } } Mark(myid,j)=1; } /* for j */ for(j=(myid*m);j<(myid*m+m-1);j++) { for(i=(j+1);i<m;i++) { sp=sqrt(a(j,(myid*m+j))*a(j,(myid*m+j))+a(i,(myid*m+j))*a(i,(myid*m+j))); c=a(j,(myid*m+j))/sp; s=a(i,(myid*m+j))/sp; for(k=0;k<M;k++) { aj[k]=c*a(j,k)+s*a(i,k); qj[k]=c*q(j,k)+s*q(i,k); ai[k]=-s*a(j,k)+c*a(i,k); qi[k]=-s*q(j,k)+c*q(i,k); } for(k=0;k<M;k++) { a(j,k)=aj[k]; q(j,k)=qj[k]; a(i,k)=ai[k]; q(i,k)=qi[k]; } } Mark(myid,j)=1; } /* for j */ Mark(myid,(myid*m+m-1))=1; } /* my_rank !=groupsize -1 */ if (myid==(p-1)) { for(j=0;j<(myid*m);j++) { //判断上一线程完成整体的某行 bcomplite=0; while(bcomplite==0) { if(Mark((myid-1),j)==1) bcomplite=1; } for(i=(myid*m);i<(myid*m+m);i++) { sp=sqrt(a(j,j)*a(j,j)+a(i,j)*a(i,j)); c=a(j,j)/sp; s=a(i,j)/sp; for(k=0;k<M;k++) { aj[k]=c*a(j,k)+s*a(i,k); qj[k]=c*q(j,k)+s*q(i,k); ai[k]=-s*a(j,k)+c*a(i,k); qi[k]=-s*q(j,k)+c*q(i,k); } for(k=0;k<M;k++) { a(j,k)=aj[k]; q(j,k)=qj[k]; a(i,k)=ai[k]; q(i,k)=qi[k]; } Mark(myid,j)=1; } /* for i */ for(k=0;k<M;k++) { Q(j,k)=q(j,k); R(j,k)=a(j,k); } }