D&C Zavolskov/Zemlyanskiy

#include <iostream>
#include <omp.h>

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include<math.h>
#include <fstream>
#include <unistd.h>
#include <sys/time.h>

using namespace std;

void foo (const char *in, const char *out)
{
    int t;
  printf ("%s -> %s: start\n", in, out);
  for( int i = 0; i < 1000000;i++){
    t +=1;
    t -=1;
  }

  /* Эмуляция того, что данная функция работает долго. На системах, отличных от linux,
   * данная функция может называться по другому */
  //sleep (1);

  printf ("%s -> %s: finish\n", in, out);
}

int copys(double ** t_or,double ** t_out,int count,int shift)
{
    for( int i = shift; i < count + shift; i++)
    {
        for( int j = shift; j < count + shift; j++)
        {
            t_out[i - shift][j-shift]=t_or[i][j];
        }
    }

    return 0;
    //из большой в маленькую.
}

double** create_array(int rows, int cols)
{
    double** array=new double*[rows];
    for (int i=0; i<rows; i++)
        array[i]=new double[cols];

    for( int i = 0 ; i < rows; i++){
        for( int j = 0; j < cols; j++)
            array[i][j] = 0;
    }
    return array;
}

void input_array(double ** array,int rows, int cols)
{
    srand(time(NULL));
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            array[i][j]=rand() % 300;
}

void input_array_s(double ** array,int rows)
{
    srand(time(NULL));
    array[0][0]=rand() % 200;
    array[1][0]=rand() % 100;
    array[0][1]=array[1][0];
    for (int i=1; i<rows - 1; i++){
        array[i][i]=array[i-1][i-1] - rand() % 200;
        array[i+1][i]=rand() % 100;
        array[i][i+1]=array[i+1][i];
    }
    array[rows-1][rows-1] = array[rows-2][rows-2] - rand() % 200;
}

void fresh_array(double ** array,int rows, int cols)
{
    srand(time(NULL));
    for (int i=0; i<rows; i++)
        for (int j=0; j<cols; j++)
            array[i][j]=0.0;
}

void print_array(double ** array, int rows, int cols)
{

    /*ofstream fout("cppstudio2.txt",ios::app);
    fout << "=";
    fout << rows;
    fout << "/";
    fout.close();*/
    for (int i=0; i<rows; i++)
    {
        for (int j=0; j<cols; j++)
        printf("%ff ", array[i][j]);
        printf("\n");
    }

}

void delete_array(double **array,int rows)
{
    for (int i=0; i<rows; i++) delete [] array[i];
    delete [] array;
}

void MultiplyWithOutAMP(double ** aMatrix,double ** bMatrix,double ** product, int x,int y,int z)
{
    for (int row = 0; row < x; row++)
    {
        for (int col = 0; col < z; col++)
        {
            for (int inner = 0; inner < y; inner++)
            {
                product[row][col] += aMatrix[row][inner] * bMatrix[inner][col];
            }
        }
    }
}

float f(float x, double ** v, double ** u, double bm, int n)
{
    double res = 0;
    for( int i = 0; i < n; i++){
        res += u[i][0]*u[i][0]/(v[i][i] - x);
    }
    return res*bm + 1;
}

float df (float x, double ** v, double ** u, double bm, int n)
{
    double res = 0;
    for( int i = 0; i < n; i++){
        res += u[i][0]*u[i][0]/((v[i][i] - x)*(v[i][i] - x));
    }
    return res*bm;
}

double newton( double x0,double err, int maxmitr, double ** v, double ** u, double bm, int n)
{
    int itr;
    double t = x0 + 0.135;
    double h,x1;
    //x0+=0.001;
    //printf("\nEnter x0, allowed error and maximum iterations\n");
    //scanf("%f %f %d", &x0, &allerr, &maxmitr);
    for (itr=1; itr<=maxmitr; itr++)
    {
        //printf("-");
        h=0.0001*f(x0,v,u,bm,n)/df(x0,v,u,bm,n);
        x1=x0-h;
        //printf(" At Iteration no. %3d, x = %9.6f\n", itr, x1);
        if (fabs(x1) > 100000){
            return t;
        }
        if (fabs(h) < err )
        {
            //printf("After %3d iterations, root = %8.6f\n", itr, x1);
            return x1;
        }
        x0=x1;
    }
    //printf(" The required solution does not converge or iterations are insufficient\n");
    return t;
}

/*

double newton( double x0,double err, int maxmitr, double ** v, double ** u, double bm, int n) {
  x0+=0.001;
  double x1  = x0 - f(x0,v,u,bm,n)/df(x0,v,u,bm,n);
  while (fabs(x1-x0)>err) {
    x0 = x1;
    x1 = x1 - f(x1,v,u,bm,n)/df(x1,v,u,bm,n);
  }
  return x1;
}
*/

void eugenVal(double ** u,double ** vt,double ** v, double bm,int n){
    //print_array(u,2,1);
    //printf("\n");
    //print_array(vt,2,1);
    //printf("\n");
    //printf("+");
    for(int i = 0;i < n; i++){
        v[i][i] = newton(vt[i][i] + 0.1,0.001,80,vt,u,bm,n);
    }
    //v[n-1][n-1] = newton(vt[n-1][n-1] +0.01,0.0001,80,vt,u,bm,n);
    return;
}
void uf(double ** q1, double ** q2, double ** u,int n){


    for( int i = 0; i < n/2; i++){
        //printf("%d - %d = %ff\n", n/2 - 1, i,u[i][0]);
        u[i][0]=q1[n/2 - 1][i];
    }

    for( int i = n/2; i < n; i++){
        u[i][0] = q2[0][i - n/2];
    }

}

int vf(double ** v1, double ** v2, double ** v,int n){


    for( int i = 0; i < n/2; i++){
        v[i][i]=v1[i][i];
    }

    for( int i = n/2; i < n; i++){
        v[i][i] = v2[i-n/2][i-n/2];
    }
    return 0;
}

void trans(double ** u, double ** ut, int n){
    for(int i = 0; i < n; i++)
        ut[0][i] = u[i][0];
    return;
}

void di(double ** d,double ** q1,double ** q2,double ** vt,double ** u,double bm, int n){
    double ** ut;
    ut = create_array(1,n);
    trans(u,ut,n);
    MultiplyWithOutAMP(u,ut,d,n,n,n);
    for(int i = 0; i < n; i++){
        d[i][i] += vt[i][i];
    }
    delete_array(ut,n);
    return;
}

void eugenVecTemp(double ** qt,double ** vt,double ** v,double ** u,int n){
    double ** temp = create_array(n,n);
    double ** temp2 = create_array(n,1);
    for(int i = 0; i < n; i++){
        fresh_array(temp,n,n);
        fresh_array(temp2,n,n);
        for(int j = 0; j < n ; j++){
            temp[j][j] = 1/(vt[j][j] - v[i][i]);
        }
            MultiplyWithOutAMP(temp,u,temp2,n,n,1);
        for(int k = 0; k < n ; k++){
            qt[k][i]=temp2[k][0];
        }
    }
    delete_array(temp,n);
    delete_array(temp2,n);
    return;
}

void eugenVec(double ** q,double **  qt,double ** q1,double ** q2, int n){
    double ** temp = create_array(n,n);

    for( int i = 0; i < n/2; i++){
        temp[i][i]=q1[i][i];
    }

    for( int i = n/2; i < n; i++){
        temp[i][i] = q2[i-n/2][i-n/2];
    }

    MultiplyWithOutAMP(temp,qt,q,n,n,n);
    delete_array(temp,n);
    return;
}

int proc(double ** t,double ** q, double ** v, int n)
{

    if( n == 1)
    {
        **q = 1;
        **v = **t;
        return 1;
    }
    else
    {

        double bm = t[n/2][n/2 - 1];

        double ** t1;
        t1 = create_array(n/2,n/2);

        double ** q1;
        q1 = create_array(n/2,n/2);

        double ** v1;
        v1 = create_array(n/2,n/2);

        copys(t,t1,n/2,0);
        t1[n/2-1][n/2-1] -= bm;

        double ** t2;
        t2 = create_array(n/2,n/2);

        double ** q2;
        q2 = create_array(n/2,n/2);


        double ** v2;
        v2 = create_array(n/2,n/2);


        copys(t,t2,n/2,n/2);

        t2[0][0] -= bm;

        omp_set_num_threads(1);

        #pragma omp parallel sections
            {
            #pragma omp section
                  {
                    proc(t1,q1,v1,n/2);
                  }
            #pragma omp section
                  {
                    proc(t2,q2,v2,n/2);
                  }

              }

        double ** u = create_array(n,1);
        uf(q1,q2,u,n);
        double ** vt = create_array(n,n);
        vf(v1,v2,vt,n);

        eugenVal(u,vt,v,bm,n);


        double ** qt = create_array(n,n);

        eugenVecTemp(qt,vt,v,u,n);

        eugenVec(q, qt,q1,q2,n);

        delete_array(q1,n/2);
        delete_array(q2,n/2);
        delete_array(v1,n/2);
        delete_array(v2,n/2);
        delete_array(t1,n/2);
        delete_array(t2,n/2);

        return 0;
    }
}


int main()
{
    double ** t;
    struct timeval start, end;
    gettimeofday(&start, NULL);
    int cout = 8;
    t = create_array(cout,cout);
    input_array_s(t,cout);


    print_array(t,cout,cout);
    double ** v;
    v = create_array(cout,cout);


    double ** q;
    q = create_array(cout,cout);


    proc(t,q,v,cout);
gettimeofday(&end, NULL);
    print_array(v,cout,cout);
    printf("\n-----------\n");
    print_array(q,cout,cout);

double delta = ((end.tv_sec  - start.tv_sec) * 1000000u +
         end.tv_usec - start.tv_usec) / 1.e6;
    printf("\n cout - %d, proc - 1, time - %ff",cout,delta);


    return 0;
}