Метод наискорейшего спуска

#include <cmath>
#include <stdio.h>
 
const int N=4;
void main() 
{
    double Q[N][N] ={{500.1, 11.2, 11.1, -13.1}, {-3.3, 500.1, 30.1, -20.1}, {7.5, 1.3, 500.1, 10}, {1.7, 7.5, -1.8, 500.1}};
    double p[N]={1.3, 1.1, 20, 1.1};
    double x0[N]={0, 0, 0, 0};
    double x1[N]={0, 0, 0, 0};
    double r[N]={0, 0, 0, 0};
    double r1[N]={0, 0, 0, 0};
    double s, s1=0;
    int k=0;
 
    do
    {
        for(int i=0; i<4; i++)
        {
            r[i]=p[i];
            for (int j=0; j<4; j++)
                r[i]-=Q[i][j]*x0[j];
        }
        
        s=0;
        for(int i=0; i<4; i++)
            s+=r[i]*r[i];
 
        for(int i=0; i<4; i++)
        {
            r1[i]=0;
            for(int j=0; j<4; j++)
                r1[i]+=Q[i][j]*r[j];
        }
 
        s1=0;
        for(int i=0; i<4; i++)
            s1+=r[i]*r1[i];
        
        s/=s1;
 
        for(int i=0; i<4; i++)
            x1[i]+=s*r[i];
        
        s=0;
        
        for(int i=0; i<4; i++)
        {
            s+=(x0[i]-x1[i])*(x0[i]-x1[i]);
            x0[i]=x1[i];
        }
        k++;
    }
    while(k<50000 && sqrt(s)>0.000001);
    printf("x1[0]:%4f\nx1[1]:%4f\nx1[2]:%4f\nx1[3]:%4f\n", x1[0],x1[1],x1[2],x1[3]);
    printf("k:%d\n",k);
}