AlgoJacobi

Jacobi Algorithm

Nov 21st, 2016
82
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 5.24 KB | None | 0 0
  1. #include <mpi.h>
  2. #include <iostream>
  3. #include "math.h"
  4. #include <iostream>
  5. #include <fstream>
  6. #include <cstdlib>
  7. #include <ctime>
  8.  
  9. using namespace std;
  10.  
  11. /*
  12. int read_data(double** a) {
  13. ifstream f_inp("matrix.txt");
  14. int n = 0;
  15. f_inp >> n;
  16. for (int j = 0; j < n; j++){
  17. for (int k = 0; k < n; k++){
  18. f_inp >> a[j][k];
  19. }
  20. }
  21. f_inp.close();
  22. return n;
  23. }
  24. */
  25.  
  26.  
  27. int main(int argc, char *argv[]) {
  28.  
  29.     int rank;
  30.     int size;
  31.     int rc;
  32.  
  33.     if (rc = MPI_Init(&argc, &argv))
  34.     {
  35.         cout << "ошибка запуска, выполнение остановлено " << endl;
  36.         MPI_Abort(MPI_COMM_WORLD, rc);
  37.     }
  38.  
  39.     MPI_Comm_size(MPI_COMM_WORLD, &size);
  40.     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  41.  
  42.     if (rank == 0)
  43.         cout << "number of processes: " << size << endl;
  44.  
  45.     if (size < 2){
  46.         MPI_Abort(MPI_COMM_WORLD, MPI_ERR_OTHER);
  47.     }
  48.  
  49.     for (int count = 0; count < 10; count++){
  50.         int n = 30 + 30 * count;
  51.         double **a;
  52.         double *d;
  53.  
  54.         d = new double[n];
  55.         a = new double*[n];
  56.  
  57.         clock_t t1, t2;
  58.  
  59.         if (rank == 0){
  60.  
  61.             for (int i = 0; i < n; i++)
  62.                 a[i] = new double[n];
  63.  
  64.             for (int i = 0; i < n; i++)
  65.             for (int j = 0; j < n; j++)
  66.                 a[i][j] = rand() % 101 - 50;
  67.             /*
  68.             int num = read_data(a);
  69.             if (num == 0){
  70.             cout << "read_data error" << endl;
  71.             system("pause");
  72.             return -1;
  73.             }
  74.             */
  75.  
  76.             t1 = clock();
  77.         }
  78. //////////////////////////////////////////////////////////
  79. //          Jacobi
  80. //////////////////////////////////////////////////////////
  81.  
  82.         double *b;
  83.         double *z;
  84.         b = new double[n];
  85.         z = new double[n];
  86.  
  87.         if (rank == 0)
  88.         for (int i = 0; i < n; ++i){
  89.             z[i] = 0.;
  90.             b[i] = d[i] = a[i][i];
  91.         }
  92.  
  93.         for (int i = 0; i < 50; ++i)
  94.         {
  95.             double sm = 0.;
  96.             double tresh;
  97.             int p, q;
  98.             bool BR = false;
  99.  
  100.             if (rank == 0){
  101.                 for (p = 0; p < n - 1; ++p)
  102.                 {
  103.                     for (q = p + 1; q < n; ++q)
  104.                         sm += fabs(a[p][q]);
  105.                 }
  106.                 if (sm == 0)
  107.                     BR = true;
  108.             }
  109.            
  110.             MPI_Bcast(&BR, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
  111.             if (BR)
  112.                 break;
  113.  
  114.             if (rank == 0)
  115.                 tresh = i < 3 ? 0.2 * sm / (n*n) : 0.;
  116.  
  117.             MPI_Bcast(&tresh, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  118.  
  119.             for (p = 0; p < n - 1; ++p)
  120.             {
  121.                 for (q = p + 1; q < n; ++q)
  122.                 {
  123.                     double apqdqdp[3];
  124.                     if (rank == 0){
  125.                         apqdqdp[0] = a[p][q];
  126.                         apqdqdp[1] = d[q];
  127.                         apqdqdp[2] = d[p];
  128.                     }
  129.                     MPI_Bcast(apqdqdp, 3, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  130.  
  131.                     double g = 1e12 * fabs(apqdqdp[0]);
  132.                     if (i >= 3 && fabs(apqdqdp[1]) > g && fabs(apqdqdp[1]) > g){
  133.                         apqdqdp[0] = 0.;
  134.                         if (rank == 0)
  135.                             a[p][q] = 0.;
  136.                     }
  137.                     else
  138.                     if (fabs(apqdqdp[0]) > tresh)
  139.                     {
  140.                         const double theta = 0.5 * (apqdqdp[1] - apqdqdp[2]) / apqdqdp[0];
  141.                         double t = 1. / (fabs(theta) + sqrt(1. + theta*theta));
  142.                         if (theta < 0)
  143.                             t = -t;
  144.                         const double c = 1. / sqrt(1. + t*t);
  145.                         const double s = t * c;
  146.                         const double tau = s / (1. + c);
  147.                         const double h = t * apqdqdp[0];
  148.                         if (rank == 0){
  149.                             z[p] -= h;
  150.                             z[q] += h;
  151.                             d[p] -= h;
  152.                             d[q] += h;
  153.                             a[p][q] = 0.;
  154.                         }
  155.  
  156.                         if (rank == 0){
  157.                             double gh[2];
  158.                             double ajpq[2];
  159.  
  160.                             for (int j = 0; j < p; j++){
  161.                                 gh[0] = a[j][p];
  162.                                 gh[1] = a[j][q];
  163.                                 MPI_Send(gh, 2, MPI_DOUBLE, j % (size - 1) + 1, j, MPI_COMM_WORLD);
  164.                             }
  165.                             for (int j = p + 1; j < q; j++){
  166.                                 gh[0] = a[p][j];
  167.                                 gh[1] = a[j][q];
  168.                                 MPI_Send(gh, 2, MPI_DOUBLE, j % (size - 1) + 1, j, MPI_COMM_WORLD);
  169.                             }
  170.                             for (int j = q + 1; j < n; j++){
  171.                                 gh[0] = a[p][j];
  172.                                 gh[1] = a[q][j];
  173.                                 MPI_Send(gh, 2, MPI_DOUBLE, j % (size - 1) + 1, j, MPI_COMM_WORLD);
  174.                             }
  175.  
  176.                             MPI_Status status;
  177.                             for (int k = 0; k < n-2; k++){
  178.                                 MPI_Recv(ajpq, 2, MPI_DOUBLE, MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
  179.                                 int jtag = status.MPI_TAG;
  180.                                 if (jtag < p){
  181.                                     a[jtag][p] = ajpq[0];
  182.                                     a[jtag][q] = ajpq[1];
  183.                                 }
  184.                                 else if (jtag > p && jtag < q){
  185.                                     a[p][jtag] = ajpq[0];
  186.                                     a[jtag][q] = ajpq[1];
  187.                                 }
  188.                                 else if (jtag > q && jtag < n){
  189.                                     a[p][jtag] = ajpq[0];
  190.                                     a[q][jtag] = ajpq[1];
  191.                                 }
  192.                             }
  193.                         }
  194.                         else
  195.                         {
  196.                             double gh[2], ajpq[2];
  197.  
  198.                             for (int j = rank - 1; j < n; j += size - 1)
  199.                             {
  200.                                 if (j != p && j != q){
  201.                                     MPI_Recv(gh, 2, MPI_DOUBLE, 0, j, MPI_COMM_WORLD, NULL);
  202.                                     ajpq[0] = gh[0] - s * (gh[1] + gh[0] * tau);
  203.                                     ajpq[1] = gh[1] + s * (gh[0] - gh[1] * tau);
  204.                                     MPI_Send(ajpq, 2, MPI_DOUBLE, 0, j, MPI_COMM_WORLD);
  205.                                 }
  206.                             }
  207.                         }
  208.                     }
  209.                 }
  210.             }
  211.  
  212.  
  213.             if (rank == 0)
  214.                 for (int p = 0; p < n; ++p){
  215.                     d[p] = (b[p] += z[p]);
  216.                     z[p] = 0.;
  217.                 }
  218.         }
  219.  
  220.         free(z);
  221.         free(b);
  222.  
  223. //////////////////////////////////////////////////////////////////////
  224.  
  225.         if (rank == 0){
  226.             t2 = clock();
  227.  
  228.             double seconds = ((double)t2 - t1) / ((double)CLOCKS_PER_SEC);
  229.  
  230.             FILE* f_outp = fopen("outp.txt", "a");
  231.             fprintf(f_outp, "* %f\n", seconds);
  232.             cout << '*' << seconds << endl;
  233.             fclose(f_outp);
  234.  
  235.             /* вывод результата: */
  236.             //for (int i = 0; i < num; i++)
  237.             //cout << d[i] << endl;
  238.         }
  239.         MPI_Barrier(MPI_COMM_WORLD);
  240.     }
  241.     MPI_Finalize();
  242. }
Add Comment
Please, Sign In to add comment