openmp

#include <fstream>
#include <cstdlib>
#include <ctime>
#include <omp.h>
#include <mpi.h>
#include <math.h>
using namespace std;
void multpar(double **a, double *v,double *c, int n,int k);
void matrix(double **a, int n);
void vecter( double *v,int n);
void  coutt(double *a, int n,int k);
double scalarproduct(double *a,double *b, int n,int k);
void  multiply(double *a, double b, double *c, int n,int k);
void  sum(double *a, double*b, double *c,char p, int n,int k);
void multpar(double **a, double *v,double *c, int n,int k);
void ort(double *w,double **v,int k, int num,int procnum);
void mainn(int num, double *time, int k);
int processor=48;
int main(int argc, char* argv[])
{  int k=10,n=1000, l=0;


  double time[7];
  ofstream fout("result.txt");
    //for(int k=700;k<=700;k=k+100)
    for(n=1000;n<=20000;n=n+1000)
  {l=0;
   mainn(n,time,k );
   for(int j=1;j<=64;j=2*j)
    fout<<"processor_number="<<j<<"   size="<<n<<"    time="<<time[l++]<<endl;
    fout<<endl;
  }
     fout.close();
}

void mainn(int num, double *time , int k)
{
     //num;количество итераци
    unsigned int start_time,end_time;
    double **v = new double* [k];
    double *w = new double [num];
    double *s = new double [num];
    double *b = new double [k];
    double *c = new double [num];
    double *a = new double [k];
    double **matr = new double* [num];
    for (int count = 0; count < k; count++)
            v[count] = new double [num];
    for (int count = 0; count < num; count++)
       {
        matr[count] = new double [num];
       }
    b[0]=0;
    vecter(v[0],num);
    multiply(v[0],1/sqrt(scalarproduct(v[0],v[0],num,2)),v[0],num,2);
    matrix(matr,num);
    int l=0;

     for (int procnum =1;  procnum<=64 ;procnum=2*procnum)
    { if(procnum==64) procnum=48;
        start_time=clock();
    for (int i = 0; i < k-1;i++)
    {
     multpar(matr,v[i],w,num,procnum);
      a[i]=scalarproduct(w,v[i],num,procnum);


     multiply(v[i],a[i],c,num,procnum);
     sum(w,c,w,'-',num,procnum);
    if (i!=0)
     {
         multiply(v[i-1],b[i],c,num,procnum);
         sum(w,c,w,'-',num,procnum);
     }

       ort(w,v,i,num,procnum);
       ort(w,v,i,num,procnum);

     b[i+1]=sqrt(scalarproduct(w,w,num,procnum));
     if (b[i+1]==0)
         {
       //  cout<<"преждевременно"<<i<<endl;
         k=i;
         break;
         }
     multiply(w,1/b[i+1],v[i+1],num,procnum);
    }
    end_time=clock();
    time[l++]=float(end_time-start_time)/CLOCKS_PER_SEC;
    }
       /*  for (int i = 0; i< k;i++)
            {
          for (int j = 0; j < k;j++)
          {
             float f=scalarproduct(v[i],v[j],num,2);
             if(f<0.001) f=0;
             cout<<f<<"  ";
         }
         cout<<endl;
         }*/
    delete  []w;
    delete  []s ;
    delete  []b ;
    delete  []c ;
    delete  []a;
    for (int count = 0; count < k; count++)
            delete[] v[count] ;
    for (int count = 0; count < num; count++)
       {
        delete [] matr[count];
       }
      delete []v;
      delete []matr;

}
void multpar(double **a, double *v,double *c, int n,int k){
      int j;
      omp_set_num_threads(k);
    #pragma omp parallel for private (j)
    for (int i = 0; i < n; i++)
    {
       c[i] = 0;
       for ( j = 0; j < n; j++)
       {
        c[i] +=a[i][j]* v[j];
       }
     }
}

void matrix(double **a, int n)
{
    srand(time(NULL));
    int j;int i;
    for ( i = 0; i < n; i++)
    for ( j = i; j < n;j++)
         {
             a[i][j]=rand()%10;
             a[j][i]=a[i][j];
             }
}
void vecter( double *v,int n)
{
    int i;
    srand(time(NULL));
    for ( i = 0; i < n; i++)
    {
         v[i]= rand()%10;
    }
}
double scalarproduct(double *a,double *b, int n,int k)
{
     omp_set_num_threads(k);
    double s=0;int j=0;
   #pragma omp parallel for private (j) reduction(+:s)
    for( j=0;j<n; j++)
       {
        s=s+a[j]*b[j];
       }

    return s;
}

void  multiply(double *a, double b, double *c, int n,int k)
{    omp_set_num_threads(k);
     #pragma omp parallel for
    for( int j=0;j<n; j++)
        c[j]=a[j]*b;
}
void  sum(double *a, double *b, double *c,char p, int n,int k)
{       omp_set_num_threads(k);
    if(p=='-')
    #pragma omp parallel for
    for( int j=0;j<n; j++)
        c[j]=a[j]-b[j];
    else
    #pragma omp parallel for
    for( int j=0;j<n; j++)
        c[j]=a[j]+b[j];
}
void  coutt(double *a, int n)
{    for( int i=0;i<n; i++)
   //cout<<a[i]<<" ";
   //cout<<endl;
   ;
}


       void ort(double *w,double **v,int k, int num,int procnum)
       {
           double  c[num], s[num];
       for(int j = 0; j<num;j++)
       {
           s[j]=0;
       }

           for (int j = 0; j<k;j++)
      {
          multiply(v[j],scalarproduct(w,v[j],num,procnum),c,num,procnum);
          sum(s,c,s,'+',num,procnum);
      }
        sum(w,s,w,'-',num,procnum);
     }