Golub-Kahan Bidiagonalization

void Bidiagonalization_G_K(int nrows, sparse_matrix_t& A, double** q)
{
    double* Aq = new double[nrows];
    double* w = new double[nrows];
    double* z = new double[nrows];
    double* p = new double[nrows];
    double* Ap = new double[nrows];

    double alpha = 0.0;
    double beta = 0.0;

    matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };

    q[0][0] = 1.0;
    for (int i = 1; i < nrows; i++) q[0][i] = 0.;

    mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[0], 0, w); //Q[0] = q_1

    alpha = eucl_norm(w, nrows);
    cblas_daxpby(nrows, 1.0 / alpha, w, 1, 0, p, 1);  //p_1 = alpha^-1 * w

    for (int i = 2; i < 10; i++)
    {
        mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, p, 0, Ap);
        cblas_daxpy(nrows, -1 * alpha, q[i-2], 1, Ap, 1);
        copy(Ap, Ap + nrows, z);

        beta = eucl_norm(z, nrows);

        cblas_daxpby(nrows, 1.0 / beta, z, 1, 0, q[i-1], 1);  //q_i+1 = beta_i^-1 * z_i

        mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[i-1], 0, Aq);
        cblas_daxpy(nrows, -1 * beta, p, 1, Aq, 1);
        copy(Aq, Aq + nrows, w);  //w

        alpha = eucl_norm(w, nrows);

        cblas_daxpby(nrows, 1.0 / alpha, w, 1, 0, p, 1);  //p_i = alpha_i^-1 * w_i
    }
    delete[] Aq;
    delete[] w;
    delete[] z;
    delete[] p;
    delete[] Ap;
}