API tools faq
paste
Advertisement
Divinty2

MME, Golub-Kahan Bidiagonalization

Divinty2
Sep 27th, 2022
812
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 11.34 KB | Source Code | 0 0
raw download clone embed print report
  1. //Нормы
  2. double eucl_norm(double* vec,
  3.                  int size) {
  4.     double res = cblas_ddot(size, vec, 1, vec, 1);
  5.     res = sqrt(res);
  6.     return res;
  7. }
  8.  
  9. double ATA_norm(sparse_matrix_t& A,
  10.                 const matrix_descr descr,
  11.                 double* vec,
  12.                 int size) {
  13.     mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, vec, 0, vec);
  14.     double res = cblas_ddot(size, vec, 1, vec, 1);
  15.     res = sqrt(res);
  16.     return res;
  17. }
  18.  
  19. //Метод минимальных ошибок (недоделанный)
  20. void MME(int nrows,               //число строк в матрице
  21.     sparse_matrix_t& A,           //матрица в формате CSR
  22.     double* u,                    //начальное приближение
  23.     double* f,                    //правая часть
  24.     double eps,                   //точность вычислений
  25.     int itermax,                  //максимальное число итераций
  26.     int* number_of_operations,    //число выполненных итераций
  27.     double* u_exact,              //точное решение
  28.     double** P
  29. )
  30. {
  31.     ofstream outf("MME.txt");
  32.     matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };
  33.  
  34.     double* Ap = new double[nrows];
  35.     //double* p = new double[nrows];
  36.     double* r = new double[nrows];
  37.     double* v_0 = new double[nrows];
  38.     double* v = new double[nrows];
  39.  
  40.     int iter = 0; double alpha = 0.0; double beta = 0.0; double rho = 0.0; double sigma = 0.0;
  41.     double sigma_1 = 0; double L = 0; double delta_0 = 0;   double rr = 0;
  42.     double Phi = 0; double vv = 0;
  43.  
  44.     copy(f, f + nrows, r);
  45.     //r_0 = f - Au_0
  46.     mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, -1., A, descr, u, 1, r);
  47.     rr = cblas_ddot(nrows, r, 1, r, 1);
  48.     eps *= eps * cblas_ddot(nrows, f, 1, f, 1);
  49.  
  50.     //Ap = Ap_0
  51.     mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, P[0], 0, Ap);
  52.  
  53.     //rho_0 = (p_0, p_0)
  54.     rho = cblas_ddot(nrows, P[0], 1, P[0], 1);
  55.     //alpha_0 = rr / rho_0
  56.     alpha = rr / rho;
  57.  
  58.     //delta_0 = ||u_exact - u_0||
  59.     copy(u_exact, u_exact + nrows, v);
  60.     cblas_daxpy(nrows, -1, u, 1, v, 1);
  61.     vv = cblas_ddot(nrows, v, 1, v, 1);
  62.     copy(v, v + nrows, v_0);
  63.     delta_0 = eucl_norm(v, nrows);
  64.     Phi = pow(delta_0, 2);  //Phi_0 = delta_0^2
  65.  
  66.     outf << "Iter: " << iter << "\n" << "u: ";
  67.     for (int i = 0; i < nrows; i++) { outf << u[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  68.     outf << "\n" << "v: ";
  69.     for (int i = 0; i < nrows; i++) { outf << v[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  70.     outf << "\n" << "r: ";
  71.     for (int i = 0; i < nrows; i++) { outf << r[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  72.     outf << "\n" << "p: ";
  73.     for (int i = 0; i < nrows; i++) { outf << P[iter][i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  74.     outf << "\n" << "Ap: ";
  75.     for (int i = 0; i < nrows; i++) { outf << Ap[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  76.     outf << "\nalpha = " << alpha << ", rho = " << rho
  77.         << ", (r_n, r_n) = " << rr << ", L = " << L
  78.         << ", delta^2 = " << vv << ", Phi = " << Phi << endl;
  79.  
  80.     while (rr > eps && iter < itermax-1)
  81.     {
  82.         ////Phi_n = delta_0^2 - summ_(k = 0, ..., n-1) L_k
  83.         sigma_1 = cblas_ddot(nrows, v_0, 1, P[iter], 1);
  84.         L = pow(sigma_1, 2) / rho;
  85.         Phi -= L;
  86.        
  87.         //u_n = u_n-1 + alpha_n-1 * p_n-1
  88.         cblas_daxpy(nrows, alpha, P[iter], 1, u, 1);
  89.         //delta_n = ||u_exact - u_n||
  90.         copy(u_exact, u_exact + nrows, v);
  91.         cblas_daxpy(nrows, -1, u, 1, v, 1);
  92.         vv = cblas_ddot(nrows, v, 1, v, 1);
  93.  
  94.         //r_n = r_n-1 - alpha_n-1 * Ap_n-1
  95.         cblas_daxpy(nrows, -1 * alpha, Ap, 1, r, 1);
  96.         rr = cblas_ddot(nrows, r, 1, r, 1);
  97.  
  98.         //sigma_n = (r_n, p_n-1)
  99.         sigma = cblas_ddot(nrows, r, 1, P[iter], 1);
  100.         //beta_n = -sigma_n / rho_n-1
  101.         beta = -1 * sigma / rho;
  102.  
  103.         ////p_n = r_n + beta_n * p_n-1
  104.         //cblas_daxpby(nrows, 1, r, 1, beta, P[iter], 1);
  105.         //mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, P[iter], 0, Ap);
  106.        
  107.         //rho_n = (p_n, p_n)
  108.         rho = cblas_ddot(nrows, P[iter], 1, P[iter], 1);
  109.         //alpha_n = -alpha_n-1 * sigma_n / rho_n
  110.         alpha = -1 * alpha * sigma / rho;
  111.  
  112.         iter++;
  113.         *number_of_operations = iter;
  114.  
  115.         outf << "Iter: " << iter << "\n" << "u: ";
  116.         for (int i = 0; i < nrows; i++) { outf << u[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  117.         outf << "\n" << "v: ";
  118.         for (int i = 0; i < nrows; i++) { outf << v[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  119.         outf << "\n" << "r: ";
  120.         for (int i = 0; i < nrows; i++) { outf << r[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  121.         outf << "\n" << "p: ";
  122.         for (int i = 0; i < nrows; i++) { outf << P[iter][i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  123.         outf << "\n" << "Ap: ";
  124.         for (int i = 0; i < nrows; i++) { outf << Ap[i] << ", "; if ((i + 1) % 7 == 0) outf << "\n"; }
  125.         outf << "\nalpha = " << alpha << ", beta = " << beta << ", rho = " << rho << ", sigma = " << sigma
  126.             << ", (r_n, r_n) = " << rr << ", L = " << L << ", delta^2 = " << vv
  127.             << ", Phi = " << Phi << endl;
  128.     }
  129.     outf.close();
  130.     delete[] r;
  131.     //delete[] p;
  132.     delete[] Ap;
  133.     delete[] v_0;
  134.     delete[] v;
  135. }
  136.  
  137.  
  138. //Бидиагонилизация Голуба-Кахана
  139.  
  140. void Bidiagonalization_G_K(
  141.     int nrows,
  142.     sparse_matrix_t& A,
  143.     double** p,
  144.     double** q,
  145.     int size_q,
  146.     double* r_0
  147. )
  148. {
  149.     double** z = new double* [size_q];
  150.     for (int i = 0; i < size_q; i++) z[i] = new double[nrows];
  151.  
  152.     double** w = new double* [size_q];
  153.     for (int i = 0; i < size_q; i++) w[i] = new double[nrows];
  154.  
  155.     double* Ap = new double[nrows];
  156.     double* Aq = new double[nrows];
  157.     double* tmp = new double[nrows];
  158.     double* alpha = new double[nrows];
  159.     double* beta = new double[nrows];
  160.  
  161.     ofstream outf("ab_check.txt");
  162.     matrix_descr descr = {SPARSE_MATRIX_TYPE_GENERAL,
  163.         SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT};
  164.     //q[0] = (1, 0, 0, 0, 0, 0, 0, 0, 0)
  165.     //q[0][0] = 10;
  166.     //for (int i = 1; i < nrows; i++) q[0][i] = 0;
  167.  
  168.     //q_0 = Ar_0
  169.     mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, r_0, 0, q[0]);
  170.     //Нормирование q_0
  171.     double norm_q_0 = eucl_norm(q[0], nrows);
  172.     cblas_daxpby(nrows, 1.0 / norm_q_0, q[0], 1, 0, tmp, 1);
  173.     copy(tmp, tmp + nrows, q[0]);
  174.  
  175.     //set w_0 = Aq_0
  176.     mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE,
  177.         1., A, descr, q[0], 0, w[0]);  
  178.     //alpha_0 = ||w_0||
  179.     alpha[0] = eucl_norm(w[0], nrows);
  180.     //p_0 = alpha_0^-1 * w_0
  181.     cblas_daxpby(nrows, 1.0 / alpha[0], w[0], 1, 0, p[0], 1);  
  182.  
  183.     for (int i = 0; i < size_q - 1; i++)
  184.     {
  185.         //z_i = Ap_i - alpha_i * q_i
  186.         mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE,
  187.             1., A, descr, p[i], 0, Ap);
  188.         cblas_daxpy(nrows, -1 * alpha[i], q[i], 1, Ap, 1);
  189.         copy(Ap, Ap + nrows, z[i]);  
  190.         //beta_i = ||z_i||
  191.         beta[i] = eucl_norm(z[i], nrows);
  192.         //q_i+1 = beta_i^-1 * z_i
  193.         cblas_daxpby(nrows, 1.0 / beta[i], z[i], 1, 0, q[i + 1], 1);
  194.  
  195.  
  196.         //w_i = Aq_i - beta_i-1 * p_i-1
  197.         mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE,
  198.             1., A, descr, q[i+1], 0, Aq);
  199.         cblas_daxpy(nrows, -1 * beta[i], p[i], 1, Aq, 1);
  200.         copy(Aq, Aq + nrows, w[i+1]);  
  201.         //alpha_i = ||w_i||
  202.         alpha[i+1] = eucl_norm(w[i+1], nrows);
  203.         //outf << "alpha_" << i+1 << " = " << setprecision(3) << alpha[i+1] << "\t";
  204.         //p_i = alpha_i^-1 * w_i
  205.         cblas_daxpby(nrows, 1.0 / alpha[i+1], w[i+1], 1, 0, p[i+1], 1);  
  206.     }
  207.     outf << "\t";
  208.     for (int i = 0; i < size_q; i++) outf << i << "\t";
  209.     outf << "\nalpha\t";
  210.     for (int i = 0; i < size_q; i++) outf << setprecision(2) << alpha[i] << "\t";
  211.  
  212.     outf << "\nbeta \t";
  213.     for (int i = 0; i < size_q - 1; i++) outf << setprecision(2) << beta[i] << "\t";
  214.  
  215.     /*for (int n = 0; n < size_q - 1; n++) {
  216.         cout << "z_" << n << " = ";
  217.         for (int i = 0; i < nrows; i++)
  218.             cout << z[n][i] << ", ";
  219.         cout << "\n";
  220.     }*/
  221.  
  222.     outf.close();
  223.  
  224.     delete[] w;
  225.     delete[] z;
  226.     delete[] Ap;
  227.     delete[] Aq;
  228.     delete[] tmp;
  229.     delete[] alpha;
  230.     delete[] beta;
  231. }
  232.  
  233.  
  234. void Bidiagonalization_ATA(
  235.     int nrows,
  236.     sparse_matrix_t& A,
  237.     double** p,
  238.     double** q,
  239.     int size_q
  240. )
  241. {
  242.     double** z = new double* [size_q];
  243.     for (int i = 0; i < size_q; i++) z[i] = new double[nrows];
  244.  
  245.     double** w = new double* [size_q];
  246.     for (int i = 0; i < size_q; i++) w[i] = new double[nrows];
  247.  
  248.     double* ATAp = new double[nrows];
  249.     double* ATAq = new double[nrows];
  250.     double* tmp = new double[nrows];
  251.     double* alpha = new double[nrows];
  252.     double* beta = new double[nrows];
  253.  
  254.     ofstream outf("ab_check.txt");
  255.     matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL,
  256.         SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };
  257.     //q[0] = (1, 0, 0, 0, 0, 0, 0, 0, 0)
  258.     q[0][0] = 1;
  259.     for (int i = 1; i < nrows; i++) q[0][i] = 0;
  260.     //set w_0 = ATA * q_0
  261.     mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[0], 0, tmp);
  262.     print(tmp, nrows, "Aq");
  263.     mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, 1., A, descr, tmp, 0, w[0]);
  264.     print(w[0], nrows, "ATAq");
  265.     //alpha_0 = ||w_0||
  266.     alpha[0] = ATA_norm(A, descr, w[0], nrows);
  267.     //alpha[0] = e_norm(w[0], nrows);
  268.     cout << "alpha = " << alpha[0];
  269.     //p_0 = alpha_0^-1 * w_0
  270.     cblas_daxpby(nrows, 1.0 / alpha[0], w[0], 1, 0, p[0], 1);
  271.  
  272.     for (int i = 0; i < size_q - 1; i++)
  273.     {
  274.         //z_i = ATA * p_i - alpha_i * q_i
  275.         mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE,
  276.             1., A, descr, p[i], 0, tmp);
  277.         mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE,
  278.             1., A, descr, tmp, 0, ATAp);
  279.         cblas_daxpy(nrows, -1 * alpha[i], q[i], 1, ATAp, 1);
  280.         copy(ATAp, ATAp + nrows, z[i]);
  281.         //beta_i = ||z_i||
  282.         beta[i] = ATA_norm(A, descr, z[i], nrows);
  283.         //beta[i] = e_norm(z[i], nrows);
  284.         //q_i+1 = beta_i^-1 * z_i
  285.         cblas_daxpby(nrows, 1.0 / beta[i], z[i], 1, 0, q[i + 1], 1);
  286.  
  287.         //w_i = ATA * q_i - beta_i-1 * p_i-1
  288.         mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1., A, descr, q[i + 1], 0, tmp);
  289.         mkl_sparse_d_mv(SPARSE_OPERATION_TRANSPOSE, 1., A, descr, tmp, 0, ATAq);
  290.         cblas_daxpy(nrows, -1 * beta[i], p[i], 1, ATAq, 1);
  291.         copy(ATAq, ATAq + nrows, w[i + 1]);
  292.         //alpha_i = ||w_i||
  293.         alpha[i + 1] = ATA_norm(A, descr, w[i + 1], nrows);
  294.         //alpha[i + 1] = e_norm(w[i + 1], nrows);
  295.         //outf << "alpha_" << i+1 << " = " << setprecision(3) << alpha[i+1] << "\t";
  296.         //p_i = alpha_i^-1 * w_i
  297.         cblas_daxpby(nrows, 1.0 / alpha[i + 1], w[i + 1], 1, 0, p[i + 1], 1);
  298.     }
  299.     outf << "\t";
  300.     for (int i = 0; i < size_q; i++) outf << i << "\t";
  301.     outf << "\nalpha\t";
  302.     for (int i = 0; i < size_q; i++) outf << setprecision(2) << alpha[i] << "\t";
  303.  
  304.     outf << "\nbeta \t";
  305.     for (int i = 0; i < size_q - 1; i++) outf << setprecision(2) << beta[i] << "\t";
  306.  
  307.     outf.close();
  308.  
  309.     delete[] w;
  310.     delete[] z;
  311.     delete[] ATAp;
  312.     delete[] ATAq;
  313.     delete[] alpha;
  314.     delete[] beta;
  315. }
  316.  
  317.  
  318. //Проверка ортогонализации
  319. void Check_orthogonalization(int nrows,
  320.     sparse_matrix_t& A,
  321.     double** v,
  322.     int size_v,
  323.     double* r_0
  324. )
  325. {
  326.     //cout << "q_0 = r_0 = ";
  327.     //for (int i = 0; i < nrows; i++)
  328.     //  cout << v[0][i] << ", ";
  329.     //cout << "\n";
  330.     //cout << "Norm_q_0 = " << eucl_norm(v[0], nrows) << "\n";
  331.  
  332.     ofstream outf("Check_ort_matr.txt");
  333.     matrix_descr descr = { SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_UPPER, SPARSE_DIAG_UNIT };
  334.  
  335.     for (int n = 0; n < size_v; n++) {
  336.         outf << "v_" << n << " = ";
  337.         for (int i = 0; i < nrows; i++)
  338.             outf << v[n][i] << ", ";
  339.         outf << "\n";
  340.     }
  341.  
  342.     double vv = 0;
  343.  
  344.     outf << "(v_i, v_j)" << "\t";
  345.     for (int i = 0; i < size_v; i++) outf << i << "\t";
  346.     outf << "\n";
  347.     for (int i = 0; i < size_v; i++) {
  348.         outf << i << "\t";
  349.         for (int j = 0; j < size_v; j++) {
  350.             vv = cblas_ddot(nrows, v[i], 1, v[j], 1);
  351.             outf << setprecision(1) << vv << "\t";
  352.         }
  353.         outf << "\n";
  354.     }
  355.     outf.close();
  356. }
Tags: C++
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!