CR_AT

void CR_AT(int nrows,
    sparse_matrix_t& A,
    double* u,
    double* f,
    double eps,
    int itermax,
    int* number_of_operations)
{
    double* Ap = new double[nrows];
    double* Atr = new double[nrows]; // A^T*r^n
    double* p = new double[nrows];
    double* r = new double[nrows];

    int iter = 0;
    double alpha = 0.0;
    double beta = 0.0;
    double rr = 0;       // (r^n, r^n)
    double Atr_p = 0;    // (A^T*r^n, p^n)
    double Atr_Atr = 0;  // (A^T*r^n, A^T*r^n)

    matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };

    eps *= eps * cblas_ddot(nrows, f, 1, f, 1);
    copy(f, f + nrows, r);
    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, -1., A, descr, u, 1, r);
    rr = cblas_ddot(nrows, r, 1, r, 1);

    mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, 1., A, descr, r, 0, Atr);
    copy(Atr, Atr + nrows, p);
    Atr_p = cblas_ddot(nrows, Atr, 1, p, 1);
    Atr_Atr = cblas_ddot(nrows, Atr, 1, Atr, 1);


    while (rr > eps && iter < itermax) {
        mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, p, 0, Ap);
        alpha = Atr_p / cblas_ddot(nrows, Ap, 1, Ap, 1);

        iter++;

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

        rr = cblas_ddot(nrows, r, 1, r, 1);
        mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, 1., A, descr, r, 0, Atr);

        beta = Atr_Atr;
        Atr_Atr = cblas_ddot(nrows, Atr, 1, Atr, 1);
        beta = Atr_Atr / beta;

        cblas_daxpby(nrows, 1, Atr, 1, beta, p, 1);
        Atr_p = cblas_ddot(nrows, Atr, 1, p, 1);

        *number_of_operations = iter;
    }

    delete[] Ap;
    delete[] Atr;
    delete[] p;
    delete[] r;
}