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