Некорректная BGK

double ATA_norm(sparse_matrix_t& A, const matrix_descr descr, double* vec, int size) {
    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, vec, 0, vec);
    double res = cblas_ddot(size, vec, 1, vec, 1);
    res = sqrt(res);
    return res;
}

void Bidiagonalization_ATA(
    int nrows,
    sparse_matrix_t& A,
    double** p,
    double** q,
    int size_q
)
{
    double** z = new double* [size_q];
    for (int i = 0; i < size_q; i++) z[i] = new double[nrows];

    double** w = new double* [size_q];
    for (int i = 0; i < size_q; i++) w[i] = new double[nrows];

    double* ATAp = new double[nrows];
    double* ATAq = new double[nrows];
    double* tmp = new double[nrows];
    double* alpha = new double[nrows];
    double* beta = new double[nrows];

    ofstream outf("ab_check.txt");
    matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL,
        SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };
    //q[0] = (1, 0, 0, 0, 0, 0, 0, 0, 0)
    q[0][0] = 1;
    for (int i = 1; i < nrows; i++) q[0][i] = 0;
    //set w_0 = ATA * q_0
    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[0], 0, tmp);
    print(tmp, nrows, "Aq");
    mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, 1., A, descr, tmp, 0, w[0]);
    print(w[0], nrows, "ATAq");
    //alpha_0 = ||w_0||
    alpha[0] = ATA_norm(A, descr, w[0], nrows);
    //alpha[0] = e_norm(w[0], nrows);
    cout << "alpha = " << alpha[0];
    //p_0 = alpha_0^-1 * w_0
    cblas_daxpby(nrows, 1.0 / alpha[0], w[0], 1, 0, p[0], 1);

    for (int i = 0; i < size_q - 1; i++)
    {
        //z_i = ATA * p_i - alpha_i * q_i
        mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE,
            1., A, descr, p[i], 0, tmp);
        mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE,
            1., A, descr, tmp, 0, ATAp);
        cblas_daxpy(nrows, -1 * alpha[i], q[i], 1, ATAp, 1);
        copy(ATAp, ATAp + nrows, z[i]);
        //beta_i = ||z_i||
        beta[i] = ATA_norm(A, descr, z[i], nrows);
        //beta[i] = e_norm(z[i], nrows);
        //q_i+1 = beta_i^-1 * z_i
        cblas_daxpby(nrows, 1.0 / beta[i], z[i], 1, 0, q[i + 1], 1);


        //w_i = ATA * q_i - beta_i-1 * p_i-1
        mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[i + 1], 0, tmp);
        mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, 1., A, descr, tmp, 0, ATAq);
        cblas_daxpy(nrows, -1 * beta[i], p[i], 1, ATAq, 1);
        copy(ATAq, ATAq + nrows, w[i + 1]);
        //alpha_i = ||w_i||
        alpha[i + 1] = ATA_norm(A, descr, w[i + 1], nrows);
        //alpha[i + 1] = e_norm(w[i + 1], nrows);
        //outf << "alpha_" << i+1 << " = " << setprecision(3) << alpha[i+1] << "\t";
        //p_i = alpha_i^-1 * w_i
        cblas_daxpby(nrows, 1.0 / alpha[i + 1], w[i + 1], 1, 0, p[i + 1], 1);
    }
    outf << "\t";
    for (int i = 0; i < size_q; i++) outf << i << "\t";
    outf << "\nalpha\t";
    for (int i = 0; i < size_q; i++) outf << setprecision(2) << alpha[i] << "\t";

    outf << "\nbeta \t";
    for (int i = 0; i < size_q - 1; i++) outf << setprecision(2) << beta[i] << "\t";

    outf.close();

    delete[] w;
    delete[] z;
    delete[] ATAp;
    delete[] ATAq;
    delete[] alpha;
    delete[] beta;
}