Advertisement
Divinty2

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. }
Advertisement
RAW Paste Data Copied
Advertisement