Метод сопряжённых градиентов в С++

//----------метод сопряженных градиентов---------------
double *SoprGrad(double **A, int n){
    cout<<"Reshenie metodom sopryjennyh gradientov:"<<endl;
    int k=0;
    double *y=new double[n+1], *x=new double[n+1],
        *ax=new double[n+1], *ap=new double[n+1],
        *z=new double[n+1], *z1=new double[n+1],
        *p=new double[n+1];
    double a, b, nz;

    //X0==
    cout<<endl<<"Vvedite nachal'noe priblizhenie dlya matrici: "<<endl;
    for(int i=0;i<n;i++)
        cin>>x;

   
    double e;
    cout<<"Vvedite to4nost':";
    cin>>e;
   
//  cout<<"Nevyazka na4al'nogo priblijeniya ravna:"<<endl;
    ax=MatrixByVector(A,x,n);
    for (i=0; i<n; i++) z=A[n]-ax;
//  PrintV(z,n);
    if (Scal(z,z,n)!=0){
    for (i=0; i<n;i++)  p=z;
    nz=1000.;
    while(nz>e){
        ap=MatrixByVector(A,p,n);
        a=Scal(z,p,n)/Scal(z,ap,n);
        for (i=0; i<n; i++) {
            y=x+a*p;
            z1=z-a*ap;
        }
        nz=NormaV(z1,n);
        b=Scal(z1,ap,n)/Scal(p,ap,n);
        for (i=0; i<n; i++) {
            p=z1-b*p;
            z=z1;
            x=y;
        }
        k++;
    }
    cout<<"k="<<k<<endl;
    return y;}
    else {
        cout<<"k=="<<1<<endl;
        return x;}
}
double* MatrixByVector(double** H,double* V,int n){
    double* tmp=new double[n];
    double sum=0;
    for (int i=0;i<n;++i){
        for (int j=0;j<n;++j)
            sum+=V[j]*H[j];
        tmp=sum;
        sum=0;}
    return tmp;}
double NormaV(double *V, int n){
    double norma;
    double s=0;
    for(int i=0;i<n;i++)
        s+=V*V;
    norma=sqrt(s);
    return norma;}
double Scal(double *v1, double *v2, int n){
    double *v=new double[n+1];
    double s=0.;
    for (int i=0; i<n; i++){
        s=v1*v2+s;}
    return s;
}