# Golub-Kahan Bidiagonalization

May 24th, 2022
927
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. void Bidiagonalization_G_K(int nrows, sparse_matrix_t& A, double** q)
2. {
3.     double* Aq = new double[nrows];
4.     double* w = new double[nrows];
5.     double* z = new double[nrows];
6.     double* p = new double[nrows];
7.     double* Ap = new double[nrows];
8.
9.     double alpha = 0.0;
10.     double beta = 0.0;
11.
12.     matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };
13.
14.     q[0][0] = 1.0;
15.     for (int i = 1; i < nrows; i++) q[0][i] = 0.;
16.
17.     mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[0], 0, w); //Q[0] = q_1
18.
19.     alpha = eucl_norm(w, nrows);
20.     cblas_daxpby(nrows, 1.0 / alpha, w, 1, 0, p, 1);  //p_1 = alpha^-1 * w
21.
22.     for (int i = 2; i < 10; i++)
23.     {
24.         mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, p, 0, Ap);
25.         cblas_daxpy(nrows, -1 * alpha, q[i-2], 1, Ap, 1);
26.         copy(Ap, Ap + nrows, z);
27.
28.         beta = eucl_norm(z, nrows);
29.
30.         cblas_daxpby(nrows, 1.0 / beta, z, 1, 0, q[i-1], 1);  //q_i+1 = beta_i^-1 * z_i
31.
32.         mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[i-1], 0, Aq);
33.         cblas_daxpy(nrows, -1 * beta, p, 1, Aq, 1);
34.         copy(Aq, Aq + nrows, w);  //w
35.
36.         alpha = eucl_norm(w, nrows);
37.
38.         cblas_daxpby(nrows, 1.0 / alpha, w, 1, 0, p, 1);  //p_i = alpha_i^-1 * w_i
39.     }
40.     delete[] Aq;
41.     delete[] w;
42.     delete[] z;
43.     delete[] p;
44.     delete[] Ap;
45. }