Conjugate Gradients

//Метод сопряжённых градиентов с использованием библиотеки MKL
void CG(
    int nrows,                  //число строк в матрице
    sparse_matrix_t& A,         //квадратная с.п.о. матрица в формате csr
    double* u,                  //вектор решения
    double* f,                  //правая часть
    double eps,                 //заданная погрешность
    int itermax,                //максимальное число итераций
    int* number_of_operations   //суммарное число выполненных итераций
)
{
    double* Ap = new double[nrows];
    double* p = new double[nrows];
    double* r = new double[nrows];

    int iter = 0;
    double alpha = 0.0;
    double beta = 0.0;
    double rr = 0;

    matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };

    copy(f, f + nrows, r);
    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, -1., A, descr, u, 1, r);

    eps *= eps * cblas_ddot(nrows, f, 1, f, 1);
    copy(r, r + nrows, p);
    rr = cblas_ddot(nrows, r, 1, r, 1);
    while (rr > eps && iter < itermax)
    {
        mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, p, 0, Ap);

        alpha = rr / cblas_ddot(nrows, Ap, 1, p, 1);

        cblas_daxpy(nrows, alpha, p, 1, u, 1);
        cblas_daxpy(nrows, -alpha, Ap, 1, r, 1);

        beta = rr;
        rr = cblas_ddot(nrows, r, 1, r, 1);
        beta = rr / beta;

        cblas_daxpby(nrows, 1, r, 1, beta, p, 1);

        iter++;
        *number_of_operations = iter;
    }
    delete[] Ap;
    delete[] p;
    delete[] r;
}