Zinovev/Pronin

#include <cstdio>
#include <iostream>
#include <fstream>
#include <math.h>
#include <ctime>
#include <chrono>
#include <omp.h>
#include <iostream>

using namespace std;

const double EPS = 0.001;

//Матричные функции
double **CreateTestMatrix(int size);
double **CopyMatrix(double** src, int size);
void MatrixEcho(double** d, int size);
double ToDiag_opt(double** d, int size, double eps);
void ToDiag_no_opt(double** d, int size, double eps);
void ToDiag_parallel(double** d, int size, double eps);
double Norm(double** d, int size);
void FindMax(int &j, int&k, double**d, int size);
void Rotate (double**d, int size, int j, int k, double c, double s);

//скалярные функции
int signum(double x)
{
    if (x>0)
        return 1;
    if (x<0)
        return -1;
    else
        return 0;
}

double abs(double x)
{
    if (x>=0)
        return x;
    else
        return -x;
}

int thread_nums;

int main(int argc, char **argv)
{

    srand(time(0));
    ofstream out;
    out.open("result.txt");

    double **d;
    int size=200;
    out<<"name\tperf\tsize\tthreads\titers\ttime"<<endl;
    for (size=400; size<=1600; size=size+200)
    {
        //cout<<size<<endl;
        double iters=0;
        d=CreateTestMatrix(size);
        for (thread_nums=1; thread_nums<=48; thread_nums++)
        {
            if ((thread_nums==1)||(thread_nums%4 == 0))
            {
                //cout<<thread_nums<<endl;


                double** d1=CopyMatrix(d,size);
                omp_set_num_threads(thread_nums);
                double start = omp_get_wtime();
                iters=ToDiag_opt(d1,size,EPS);
                double end = omp_get_wtime();
                double sum_time=end-start;
                double R=iters*(double)(size*size*size*6) / (sum_time);
                out<<"\t"<<R<<"\t"<<size<<"\t"<<thread_nums<<"\t"<<iters<<"\t"<<sum_time<< endl;
                for (int i=0;i<size;i++)
                {
                    delete[] d1[i];
                }
                delete[] d1;

            }
        }


    }
    out.close();
    return 0;
}

double** CreateTestMatrix(int size)
{
    double **d = new double*[size];
    for (int k = 0; k < size; k++)
    {
        d[k] = new double[size];
    }
    int i,j;

    for (i = 0; i < size; i++){
        double res[size];
    //  #pragma omp paralell for private(i,size)
        for (j = i; j < size; j++)
            {
                //d[i][j] = Random(-100, 100);

                res[j] = (double)(rand()) / RAND_MAX;
                res[j]=res[j]*(100 - (-100)) + (-100);

            }

            for (j = i; j < size; j++)
            {
                d[i][j]=res[j];
            }

    }




    for (i = 0; i < size; i++)
    {
        for (j = 0; j <=i; j++)
        {
            d[i][j] = d[j][i];
        }

    }

    return d;
}

double **CopyMatrix (double** src, int size)
{
    double **d = new double*[size];
    for (int k = 0; k < size; k++)
    {
        d[k] = new double[size];
    }
    int i,j;
    for (i=0;i<size;i++)
    {
        for (j=0;j<size;j++)
        {
            d[i][j]=src[i][j];
        }
    }
    return d;

}


void MatrixEcho(double** d, int size=0)
{
    for (int i=0;i<size;i++)
    {
        for (int j=0;j<size;j++)
        {
            if (abs(d[i][j])>0.0001)
                cout<<d[i][j]<<" ";
            else
                cout<<"0"<<" ";
        }
        cout<<endl;
    }
}



void ToDiag_no_opt(double **d, int size, double eps)
{
    double norm= Norm(d, size);
    long iters=0;
    while (norm>eps)
    {
        iters++;
        int j=0,k=1;
        //Нашли наибольший элемент
        FindMax(j, k, d,size);
        norm=sqrt(norm*norm-d[j][k]*d[j][k]);
        if (abs(d[j][k])>0.0001)
        {
            double tau,t,c,s;
            if(d[j][j]!=d[k][k])
            {
                tau=(d[j][j] - d[k][k])/(2*d[j][k]);
                t=signum(tau)/(abs(tau) + sqrt(1+tau*tau));
                c=1/sqrt(1+t*t);
                s=c*t;
            }
            else
            {
                c=1/(sqrt(2));
                s=c;
            }
            Rotate(d,size, j,k, c,s);

        }



    }
    cout<<iters<<" "<<endl;
}

double ToDiag_opt(double **d, int size, double eps)
{
    double norm= Norm(d, size);
    double iters=0;
    while (norm>eps)
    {
        for (int j=0;j<size;j++)
        {
            for (int k=j+1;k<size;k++)
            {
            iters++;


            norm=sqrt(norm*norm-d[j][k]*d[j][k]);


            if (abs(d[j][k])>0.0001)
            {
                double tau,t,c,s;
                if(d[j][j]!=d[k][k])
                {
                    tau=(d[j][j] - d[k][k])/(2*d[j][k]);
                    t=signum(tau)/(abs(tau) + sqrt(1+tau*tau));
                    c=1/sqrt(1+t*t);
                    s=c*t;
                }
                else
                {
                    c=1/(sqrt(2));
                    s=c;
                }

                Rotate(d,size, j,k, c,s);
                if (norm<eps)
                break;

            }
            }
            if (norm < eps)
                break;
        }



    }
    //cout<<iters<<" "<<endl;
    return iters;
}

void Rotate (double**d, int size, int j, int k, double c, double s)
{
    #pragma omp parallel
    {

    #pragma omp for
    /*for (int q=0;q<size/thread_nums;q++ )
    {*/
        for (int i=0;i<size;i++)
        {
            /*int cur_index=q*thread_nums+i;
            if (cur_index>=size)
                break;*/
            double d1=d[j][i];
            double d2=d[k][i];
            double c1=c,s1=s;
            int j1=j, k1=k;
            d[j1][i]=d1*c1+d2*s1;
            d[k1][i]=-s1*d1+c1*d2;
            if ((i!=j1)&&(i!=k1))
            {
                d[i][j1]=d[j1][i];
                d[i][k1]=d[k1][i];
            }
        //cout<<i<<" ";
        }
    }
    //}
    //cout<<endl;
    double d1=d[j][j];
    double d2=d[k][j];
    d[j][j]=d1*c+d2*s;
    d[k][j]=0;

    d1=d[j][k];
    d2=d[k][k];
    d[j][k]=0;
    d[k][k]=-s*d1+c*d2;
}

double Norm(double** d ,int size)
{
    double result=0;
    for (int i=0;i<size;i++)
    {
        for (int j=i+1;j<size;j++)
        {
            result+=d[i][j]*d[i][j];
        }
    }
    result=sqrt(result);
    return result;

}

void FindMax(int &j, int&k, double**d, int size)
{
    double max=d[j][k];
    for (int p=0;p<size;p++)
    {
        for (int q=p+1;q<size;q++)
        {
            if (abs(d[p][q])>max)
            {
                max=abs(d[p][q]);
                j=p;
                k=q;
            }
        }
    }
}