Умножение матриц по Винограду

#include "stdafx.h"
#include <iostream>//i/o
#include <conio.h> //getch
#include <cmath>   //fabs, abs
using namespace std;
 
void ShowVector(int n, double * vec);
void PryamoiHod(int n, double **a, double *b);
void ObratniHod(int n, double **a, double *b, double *x);
 
int main()
{
        int i,j,n;
        double **a, *b, *x;
        do
        {
                std::cout<<"Enter NUM of equations: ";
                std::cin>>n;
                //Выделяем память под матрицу А и векторы В и Х
                a = new double *[n];
                b = new double  [n];
                x = new double  [n];
                for(i = 0; i < n; i++)
                {
                        a[i] = new double[n];
                        //Ввод a
                        for(j = 0; j < n; j++)
                        {
                                std::cout<<"a["<<i + 1<<"]["<<j + 1<<"] = ";
                                std::cin>>a[i][j];
                        }
                }
                //Ввод b
                for(i = 0; i < n; i++)
                {
                        std::cout<<"b["<<i + 1<<"] = ";
                        std::cin>>b[i];
                }
                
                std::cout<<"\tSee input\r\n";
                std::cout<<"Matrix A:\r\n";
                for(i = 0; i < n; i++)
                        ShowVector(n, a[i]);
                std::cout<<"Vector B:\r\n";
                ShowVector(n, b);
                
                std::cout<<"\tSolving on Gauss method\r\n";
                PryamoiHod(n, a, b);
                std::cout<<"Forvard Gauss course\r\n";//Прямой ход
                std::cout<<"Matrix A:\r\n";
                for(i = 0; i < n; i++)
                        ShowVector(n, a[i]);
                std::cout<<"Vector B:\r\n";
                ShowVector(n, b);
 
                ObratniHod(n, a, b, x);
                std::cout<<"Back Gauss course\r\n";//Обратный ход
                std::cout<<"Matrix A:\r\n";
                for(i = 0; i < n; i++)
                        ShowVector(n, a[i]);
                std::cout<<"Vector B:\r\n";
                ShowVector(n, b);
 
                std::cout<<"Results :\r\n";
                ShowVector(n, x);
 
                std::cout<<"Press Y for new input\r\n";
                //Чистим память
                delete [] a;
                delete [] b;
                delete [] x;
        }
        while(toupper(getch()) == 'Y');
        return 0;
}
 
void ShowVector(int n, double * vec)
{
        for(int i = 0; i < n; i++)
                std::cout<<vec[i]<<" ";
        std::cout<<endl;
}
 
void PryamoiHod(int n, double **a, double *b)
{
        double v;
        for(int k = 0,i,j,im; k < n - 1; k++)
        {
                im = k;
                for(i = k + 1; i < n; i++)
                {
                        if(abs(a[im][k]) < abs(a[i][k]))
                        {
                                im = i;
                        }
                }
                if(im != k)
                {
                        for(j = 0; j < n; j++)
                        {
                                v                = a[im][j];
                                a[im][j] = a[k][j];
                                a[k][j]  = v;
                        }
                        v     = b[im];
                        b[im] = b[k];
                        b[k]  = v;
                }
                for(i = k + 1; i < n; i++)
                {
                        v               = 1.0*a[i][k]/a[k][k];
                        a[i][k] = 0;
                        b[i]    = b[i] - v*b[k];
                                                if(v != 0)
                        for(j = k + 1; j < n; j++)
                        {
                                a[i][j] = a[i][j] - v*a[k][j];
                        }
                }
        }
}
 
void ObratniHod(int n, double **a, double *b, double *x)
{
        double s = 0;
        x[n - 1] = 1.0*b[n - 1]/a[n - 1][n - 1];
        for(int i = n - 2, j; 0 <= i; i--)
        {
                s = 0;
                for(j = i + 1; j < n; j++)
                {
                        s = s+a[i][j]*x[j];
                }
                x[i] = 1.0*(b[i] - s)/a[i][i];
        }
}