#include "stdafx.h"
#include <cmath>
#include <mpi.h>
#include <fstream>
#include <iostream>



long double xZeroFunction()
{
    return 0.0;
}
long double xOneFunction(long double t)
{
    return pow(2.0, 3.0 / 4.0)*exp(4.0*t);
}
long double CorrectXTFunction(long double x, long double t)
{
    return exp(4.0*t)*pow(2.0*x, 3.0 / 4.0);
}
long double StepXFunction(long double wiPlus1, long double wiMinus1, long double wi, long double tau, long double h)
{
    return wi + tau * (4.0*wi + ((3.0*wi*((wiPlus1 - 2.0*wi + wiMinus1) / pow(h, 2.0)) + pow((wiPlus1 - wi) / h, 2.0)) / (3.0*pow(wi, 2.0 / 3.0))));
}
long double ZeroZFunction(long double x)
{
    return pow(2.0*x, 3.0 / 4.0);
}



int main(int argc, char* argv[])
{
    const int numberOfStepsX = 20;
    const int numberOfStepsT = 2000;

    const long double stepX = 1.0 / numberOfStepsX;
    const long double stepT = 1.0 / numberOfStepsT;

    MPI_Init(&argc, &argv);

    int myid, numprocs;
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);

    long double matrixOfCorrectAnswers[numberOfStepsT][numberOfStepsX] = { 0 };
    long double matrixOfApproximateAnswers[numberOfStepsT][numberOfStepsX] = { 0 };

    long double currentX = 0.0;
    long double currentT = 0.0;
    double t1, t2;
    if (myid == 0)
    {
        for (int i = 0; i < numberOfStepsT; i++)
        {
            currentX = 0.0;
            for (int j = 0; j < numberOfStepsX; j++)
            {
                matrixOfCorrectAnswers[i][j] = CorrectXTFunction(currentX, currentT);
                currentX += stepX;
            }
            currentT += stepT;
        }
        // Algorithm starts
        t1 = MPI_Wtime();
        currentX = 0.0;
        for (int i = 0; i < numberOfStepsX; i++)
        {
            matrixOfApproximateAnswers[0][i] = ZeroZFunction(currentX);
            currentX += stepX;
        }

        currentT = 0.0;
        for (int i = 0; i < numberOfStepsT; i++)
        {
            matrixOfApproximateAnswers[i][0] = xZeroFunction();
            matrixOfApproximateAnswers[i][numberOfStepsX - 1] = xOneFunction(currentT);
            currentT += stepT;
        }
    }
    for (int i = 1; i < numberOfStepsT; i++)
    {
        MPI_Bcast(matrixOfApproximateAnswers[i - 1], numberOfStepsX, MPI_LONG_DOUBLE, 0, MPI_COMM_WORLD);
        for (int j = myid + 1; j < numberOfStepsX - 1; j += numprocs)
        {
            matrixOfApproximateAnswers[i][j] = StepXFunction(matrixOfApproximateAnswers[i - 1][j + 1],
                matrixOfApproximateAnswers[i - 1][j - 1], matrixOfApproximateAnswers[i - 1][j], stepT, stepX);
        }
        if (myid == 0)
        {
            for (int k = 1; k < numprocs; k++)
            {
                long double temp[numberOfStepsX];
                MPI_Recv(temp, numberOfStepsX, MPI_LONG_DOUBLE, k, MPI_ANY_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                for (int j = k; j < numberOfStepsX - 1; j += numprocs)
                {
                    matrixOfApproximateAnswers[i][j] = temp[j];
                }
            }
        }
        else
        {
            MPI_Send(&matrixOfApproximateAnswers[i], numberOfStepsX, MPI_LONG_DOUBLE, 0, 0, MPI_COMM_WORLD);
        }
    }

    if (myid == 0)
    {
        // Algorithm finishes
        t2 = MPI_Wtime();
        printf("Elapsed time is %f\n", t2 - t1);
        // Saving results to file
        std::ofstream fout1("/Users/Nazar Shevchuk/Kursova/results.txt");
        std::ofstream fout2("/Users/Nazar Shevchuk/Kursova/results_yavno.txt");
        for (int i = 0; i < numberOfStepsT; i++)
        {
            for (int j = 0; j < numberOfStepsX; j++)
            {
                //printf("%3.20Lf - %3.20Lf\n", matrixOfCorrectAnswers[i][j], matrixOfApproximateAnswers[i][j]);
                fout1 << matrixOfCorrectAnswers[i][j] << '\t';
                fout2 << matrixOfApproximateAnswers[i][j] << '\t';
            }
            fout1 << std::endl;
            fout2 << std::endl;
        }
        fout1.close();
        fout2.close();
    }
    MPI_Finalize();
    System ("pause");
    return 0;
}