Jacobi_Alex

#include <iostream>
#include <cstdlib>
#include <stdlib.h>
#include <cstdlib>
#include <fstream>
#include <cmath>
#include <stdio.h>
#include <iomanip>
#include <omp.h>

using namespace std;

int  sign(double a)
{
    if (a > 0) return 1;
    else return -1;
}

void PrintArr(double **arr, int n)
{
    for(int i=0; i<n; i++)
    {
        for(int j=0; j<n; j++)
        if (arr[i][j] < 10) cout<<setw(13)<< arr[i][j]<<"  ";
        else cout<< setw(13)<<arr[i][j]<<" ";
        cout<<endl;
    }
}

double **GenArr(int n)
{
    double **arr = new double*[n];
    for (int i=0; i < n; i++)
        arr[i] = new double[n];

    for (int i=0; i < n; i++)
        for (int j=i; j < n; j++)
            arr[i][j] = arr[j][i] = rand()%100;

    return arr;
}

double **CopyArr(double **arr, int n)
{
    double **arr1=new double*[n];
    for (int i=0; i < n; i++)
        arr1[i] = new double[n];

    for (int i=0; i < n; i++)
        for (int j=i; j < n; j++)
            arr1[i][j] = arr1[j][i] = arr[i][j];


    return arr1;
}

int sizef(int size)
{
    if (size == 200) return 300;
    else if (size == 300) return 500;
         else return size+500;
}

int treadnumf(int treadnum)
{
    if (treadnum == 32) return 48;
    else return treadnum*2;
}

int main(int argc, char **argv)
{
    int n=10;
    int iternum=0;
    int treadnum=1;
    double eps=0.001, c=0, s=0, r=0, t=0, tmp3=0, tmp4=0, tmp5=0, norm=0, stTime=0, fnTime=0, time=0;
    double **arr;
    //int k1=0,k2=0,k3=0,k4=0;
    ofstream fout("out.txt");

    fout<<endl<<setw(13)<<"size"<<setw(13)<<"treads"<<setw(13)<<"iteration"<<setw(13)<<"time"<<endl;

    for (n=200; n <= 2000; n=sizef(n))
    {
        /*Создание массива*/
        arr = GenArr(n);

        /*Вывод массива*/
        //PrintArr(arr,n);

        for (treadnum=1; treadnum <= 48; treadnum=treadnumf(treadnum))
        {
            double **arr1=CopyArr(arr,n);
            stTime = omp_get_wtime();
            iternum = 0;

            /*Подсчёт нормы*/
            norm=0;
            for(int i=0; i < n; i++)
                for(int j=i + 1; j < n; j++)
                    norm = norm + arr1[i][j] * arr1[i][j];

            /*Цикл итераций алгоритма*/
            while (1)
            {
                /*Проход по элементам зануления*/
                for(int ik=0; ik < n; ik++)
                for(int jk=ik + 1; jk < n; jk++)
                {
                    iternum++;
                    /*Условие выхода*/
                    if (norm <= eps) goto stop;

                    /*Вычисление параметров*/
                    if (arr1[ik][ik] == arr1[jk][jk]) c = s = sqrt(2) / 2;
                    else
                    {
                        r = (arr1[jk][jk] - arr1[ik][ik]) / (2 * arr1[ik][jk]);
                        t = sign(r) / (abs(r) + sqrt(1 + r * r));
                        c = 1 / (sqrt(1 + t * t));
                        s = c * t;
                    }

                    /*Пересчёт элементов*/
                    tmp3 = arr1[ik][ik];
                    tmp4 = arr1[jk][jk];
                    tmp5 = arr1[ik][jk];
                        omp_set_num_threads(treadnum);
                        #pragma omp parallel for
                        for(int i=0; i < n; i++)
                        {
                            double tmp1 = arr1[ik][i];
                            double tmp2 = arr1[jk][i];
                            arr1[i][ik] = arr1[ik][i] = c * tmp1 - s * tmp2;
                            arr1[i][jk] = arr1[jk][i] = s * tmp1 + c * tmp2;
                            //if (omp_get_thread_num() == 0) k1++;
                            //if (omp_get_thread_num() == 1) k2++;
                            //if (omp_get_thread_num() == 2) k3++;
                            //if (omp_get_thread_num() == 3) k4++;
                        }
                    arr1[ik][ik] = c * c * tmp3 - 2 * c * s * tmp5 + s * s * tmp4;
                    arr1[jk][jk] = s * s * tmp3 + 2 * c * s * tmp5 + c * c * tmp4;
                    arr1[ik][jk] = arr1[jk][ik] = 0;


                    /*Измеенение нормы*/
                    norm = norm - tmp5 * tmp5;
                }
            }

            stop:

            fnTime = omp_get_wtime();
            time = fnTime - stTime;
            fout<<endl<<setw(13)<<n<<setw(13)<<treadnum<<setw(13)<<iternum<<setw(13)<<time<<endl;
            //cout<<"k1="<<k1<<"  k2="<<k2<<"  k3="<<k3<<"  k4="<<k4<<endl;
            //PrintArr(arr,n);
            //PrintArr(arr1,n);

            /*Освобождение памяти*/
            for (int i=0; i < n; i++)
                delete[] arr1[i];
            delete[] arr1;
        }
        for (int i=0; i < n; i++)
            delete[] arr[i];
        delete[] arr;
    }

    fout.close();
    system ("pause");
    return 0;
}