SHARE
TWEET

Untitled

a guest Aug 10th, 2017 201 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2.  
  3. /*
  4. ************************************************************************
  5.                         Example 25 (congrad.c)
  6.  
  7.        Objective            : Conjugate Gradient Method to solve AX = b
  8.                               matrix system of linear equations.
  9.  
  10.        Input                : Real Symmetric Positive definite Matrix
  11.                               and the real vector - Input_A and Vector_B
  12.                               Read files (mdatcg.inp) for Matrix A
  13.                               and (vdatcg.inp) for Vector b
  14.  
  15.        Description          : Input matrix is stored in n by n format.
  16.                               Diagonal preconditioning matrix is used.
  17.                               Rowwise block striped partitioning matrix
  18.                               is used.Maximum iterations is given by
  19.                               MAX_ITERATIONS.Tolerance value is given
  20.                               by EPSILON
  21.                               Header file used  : cg_constants.h
  22.  
  23.        Output               : The solution of  Ax=b on process with
  24.                               Rank 0 and the number of iterations
  25.                               for convergence of the method.
  26.  
  27.        Necessary conditions : Number of Processes should be less than
  28.                               or equal to 8 Input Matrix Should be
  29.                               Square Matrix. Matrix size for Matrix A
  30.                               and vector size for vector b should be
  31.                               equally striped, that is Matrix size and
  32.                               Vector Size should be dividible by the
  33.                               number of processes used.
  34.  
  35. ***************************************************************
  36. */
  37.  
  38.  
  39.  
  40. #include<math.h>
  41. #include<stdio.h>
  42. #include<stdlib.h>
  43. #include<mpi.h>
  44.  
  45. #define EPSILON           1.0E-20
  46. #define MAX_ITERATIONS 10000
  47.  
  48. /******************************************************************************/
  49. void
  50. GetPreconditionMatrix(double **Bloc_Precond_Matrix, int NoofRows_Bloc,
  51.               int NoofCols)
  52. {
  53.  
  54.     /*... Preconditional Martix is identity matrix .......*/
  55.     int     Bloc_MatrixSize;
  56.     int     irow, icol, index;
  57.     double *Precond_Matrix;
  58.  
  59.     Bloc_MatrixSize = NoofRows_Bloc*NoofCols;
  60.  
  61.     Precond_Matrix = (double *) malloc(Bloc_MatrixSize * sizeof(double));
  62.  
  63.    index = 0;
  64.     for(irow=0; irow<NoofRows_Bloc; irow++){
  65.         for(icol=0; icol<NoofCols; icol++){
  66.             Precond_Matrix[index++] = 1.0;
  67.         }
  68.     }
  69.     *Bloc_Precond_Matrix = Precond_Matrix;
  70. }
  71. /******************************************************************************/
  72. double
  73. ComputeVectorDotProduct(double *Vector1, double *Vector2, int VectorSize)
  74. {
  75.     int     index;
  76.     double Product;
  77.  
  78.     Product = 0.0;
  79.     for(index=0; index<VectorSize; index++)
  80.         Product += Vector1[index]*Vector2[index];
  81.  
  82.     return(Product);
  83. }
  84. /******************************************************************************/
  85. void
  86. CalculateResidueVector(double *Bloc_Residue_Vector, double *Bloc_Matrix_A,
  87.                               double *Input_B, double *Vector_X, int NoofRows_Bloc,
  88.                               int VectorSize, int MyRank)
  89. {
  90.     /*... Computes residue = AX - b .......*/
  91.     int   irow, index, GlobalVectorIndex;
  92.     double value;
  93.  
  94.     GlobalVectorIndex = MyRank * NoofRows_Bloc;
  95.     for(irow=0; irow<NoofRows_Bloc; irow++){
  96.         index = irow * VectorSize;
  97.         value = ComputeVectorDotProduct(&Bloc_Matrix_A[index], Vector_X,
  98.                                                   VectorSize);
  99.         Bloc_Residue_Vector[irow] = value - Input_B[GlobalVectorIndex++];
  100.     }
  101. }
  102. /******************************************************************************/
  103. void
  104. SolvePrecondMatrix(double *Bloc_Precond_Matrix, double *HVector,
  105.                          double *Bloc_Residue_Vector, int Bloc_VectorSize)
  106. {
  107.     /*...HVector = Bloc_Precond_Matrix inverse * Bloc_Residue_Vector.......*/
  108.     int index;
  109.  
  110.     for(index=0; index<Bloc_VectorSize; index++){
  111.         HVector[index] = Bloc_Residue_Vector[index]/1.0;
  112.     }
  113. }
  114. /******************************************************************************/
  115. main(int argc, char *argv[])
  116. {
  117.     int        NumProcs, MyRank, Root=0;
  118.     int        NoofRows, NoofCols, VectorSize;
  119.     int        n_size, NoofRows_Bloc, Bloc_MatrixSize, Bloc_VectorSize;
  120.     int      Iteration = 0, irow, icol, index, CorrectResult;
  121.     double     **Matrix_A, *Input_A, *Input_B, *Vector_X, *Bloc_Vector_X;
  122.     double      *Bloc_Matrix_A, *Bloc_Precond_Matrix, *Buffer;
  123.     double    *Bloc_Residue_Vector ,*Bloc_HVector, *Bloc_Gradient_Vector;
  124.     double        *Direction_Vector, *Bloc_Direction_Vector;
  125.     double      Delta0, Delta1, Bloc_Delta0, Bloc_Delta1;
  126.     double        Tau, val, temp, Beta;
  127.  
  128.     double      *AMatDirect_local, *XVector_local;
  129.     double      *ResidueVector_local, *DirectionVector_local;
  130.     double     StartTime, EndTime;
  131.     MPI_Status status;
  132.     FILE      *fp;
  133.  
  134.     /*...Initialising MPI .......*/
  135.     MPI_Init(&argc,&argv);
  136.     MPI_Comm_size(MPI_COMM_WORLD,&NumProcs);
  137.     MPI_Comm_rank(MPI_COMM_WORLD,&MyRank);
  138.  
  139.  
  140.   /* .......Read the Input file ......*/
  141.   if(MyRank == Root) {
  142.  
  143.      if ((fp = fopen ("./data/mdatcg.inp", "r")) == NULL) {
  144.        printf("Can't open input matrix file");
  145.        exit(-1);
  146.      }
  147.      fscanf(fp, "%d %d", &NoofRows,&NoofCols);    
  148.      n_size=NoofRows;
  149.  
  150.      /* ...Allocate memory and read data .....*/
  151.       Matrix_A = (double **) malloc(n_size*sizeof(double *));
  152.  
  153.      for(irow = 0; irow < n_size; irow++){
  154.           Matrix_A[irow] = (double *) malloc(n_size * sizeof(double));
  155.         for(icol = 0; icol < n_size; icol++)
  156.             fscanf(fp, "%lf", &Matrix_A[irow][icol]);
  157.       }
  158.      fclose(fp);
  159.  
  160.      if ((fp = fopen ("./data/vdatcg.inp", "r")) == NULL){
  161.         printf("Can't open input vector file");
  162.         exit(-1);
  163.      }
  164.  
  165.      fscanf(fp, "%d", &VectorSize);    
  166.      n_size=VectorSize;
  167.      Input_B  = (double *)malloc(n_size*sizeof(double));
  168.      for (irow = 0; irow<n_size; irow++)
  169.         fscanf(fp, "%lf",&Input_B[irow]);
  170.      fclose(fp);
  171.  
  172.       /* ...Convert Matrix_A into 1-D array Input_A ......*/
  173.      Input_A  = (double *)malloc(n_size*n_size*sizeof(double));
  174.       index    = 0;
  175.       for(irow=0; irow<n_size; irow++)
  176.           for(icol=0; icol<n_size; icol++)
  177.               Input_A[index++] = Matrix_A[irow][icol];
  178.  
  179.   }
  180.  
  181.     MPI_Barrier(MPI_COMM_WORLD);
  182.     StartTime = MPI_Wtime();
  183.  
  184.     /*...Broadcast Matrix and Vector size and perform input validation tests...*/
  185.     MPI_Bcast(&NoofRows, 1, MPI_INT, Root, MPI_COMM_WORLD);
  186.     MPI_Bcast(&NoofCols, 1, MPI_INT, Root, MPI_COMM_WORLD);
  187.     MPI_Bcast(&VectorSize, 1, MPI_INT, Root, MPI_COMM_WORLD);
  188.  
  189.  
  190.    if(NoofRows != NoofCols){
  191.         MPI_Finalize();
  192.         if(MyRank == Root)
  193.             printf("Error : Coefficient Matrix Should be square matrix");
  194.         exit(-1);
  195.     }
  196.  
  197.    if(NoofRows != VectorSize){
  198.         MPI_Finalize();
  199.         if(MyRank == Root)
  200.             printf("Error : Matrix Size should be equal to VectorSize");
  201.         exit(-1);
  202.     }
  203.  
  204.     if(NoofRows % NumProcs != 0){
  205.         MPI_Finalize();
  206.         if(MyRank == Root)
  207.             printf("Error : Matrix cannot be evenly striped among processes");
  208.         exit(-1);
  209.     }
  210.  
  211.    /*...Allocate memory for Input_B and BroadCast Input_B.......*/
  212.    if(MyRank != Root)
  213.         Input_B = (double *) malloc(VectorSize*sizeof(double));
  214.     MPI_Bcast(Input_B, VectorSize, MPI_DOUBLE, Root, MPI_COMM_WORLD);
  215.  
  216.    /*...Allocate memory for Block Matrix A and Scatter Input_A .......*/
  217.    NoofRows_Bloc   = NoofRows / NumProcs;
  218.     Bloc_VectorSize = NoofRows_Bloc;
  219.     Bloc_MatrixSize = NoofRows_Bloc * NoofCols;
  220.     Bloc_Matrix_A = (double *) malloc (Bloc_MatrixSize*sizeof(double));
  221.     MPI_Scatter(Input_A, Bloc_MatrixSize, MPI_DOUBLE,
  222.                    Bloc_Matrix_A, Bloc_MatrixSize, MPI_DOUBLE, Root, MPI_COMM_WORLD);
  223.  
  224.  
  225.     /*... Allocates memory for solution vector and intialise it to zero.......*/
  226.     Vector_X = (double *) malloc(VectorSize*sizeof(double));
  227.     for(index=0; index<VectorSize; index++)
  228.         Vector_X[index] = 0.0;
  229.  
  230.     /*...Calculate RESIDUE = AX - b .......*/
  231.    Bloc_Residue_Vector = (double *) malloc(NoofRows_Bloc*sizeof(double));
  232.     CalculateResidueVector(Bloc_Residue_Vector,Bloc_Matrix_A, Input_B, Vector_X,
  233.                                   NoofRows_Bloc, VectorSize, MyRank);
  234.  
  235.     /*... Precondtion Matrix is identity matrix ......*/
  236.     GetPreconditionMatrix( &Bloc_Precond_Matrix, NoofRows_Bloc, NoofCols);
  237.  
  238.     /*...Bloc_HVector = Bloc_Precond_Matrix inverse * Bloc_Residue_Vector......*/
  239.     Bloc_HVector = (double *) malloc(Bloc_VectorSize*sizeof(double));
  240.     SolvePrecondMatrix(Bloc_Precond_Matrix, Bloc_HVector, Bloc_Residue_Vector,
  241.                              Bloc_VectorSize);
  242.  
  243.  
  244.    /*...Initailise Bloc Direction Vector = -(Bloc_HVector).......*/
  245.     Bloc_Direction_Vector = (double *) malloc(Bloc_VectorSize*sizeof(double));
  246.     for(index=0; index<Bloc_VectorSize; index++)
  247.         Bloc_Direction_Vector[index] = 0 - Bloc_HVector[index];
  248.  
  249.     /*...Calculate Delta0 and check for convergence .......*/
  250.     Bloc_Delta0 = ComputeVectorDotProduct(Bloc_Residue_Vector,
  251.                                         Bloc_HVector, Bloc_VectorSize);
  252.     MPI_Allreduce(&Bloc_Delta0, &Delta0, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  253.  
  254.     if(Delta0 < EPSILON){
  255.         MPI_Finalize();
  256.         exit(0);
  257.     }
  258.  
  259.     /*...Allocate memory for Direction Vector.......*/
  260.     Direction_Vector = (double *) malloc(VectorSize*sizeof(double));
  261.  
  262.     /*...Allocate temporary buffer to store Bloc_Matrix_A*Direction_Vector...*/
  263.     Buffer = (double *) malloc(Bloc_VectorSize*sizeof(double));
  264.  
  265.     /*...Allocate memory for Bloc_Vector_X .......*/
  266.     Bloc_Vector_X = (double *)malloc(Bloc_VectorSize*sizeof(double));
  267.  
  268.     Iteration = 0;
  269.     do{
  270.  
  271.         Iteration++;
  272.         /*
  273.         if(MyRank == Root)
  274.             printf("Iteration : %d\n",Iteration);
  275.         */
  276.  
  277.         /*...Gather Direction Vector on all processes.......*/
  278.         MPI_Allgather(Bloc_Direction_Vector, Bloc_VectorSize, MPI_DOUBLE,
  279.                     Direction_Vector, Bloc_VectorSize, MPI_DOUBLE, MPI_COMM_WORLD);
  280.  
  281.         /*...Compute Tau = Delta0 / (DirVector Transpose*Matrix_A*DirVector)...*/
  282.         for(index=0; index<NoofRows_Bloc; index++){
  283.             Buffer[index] = ComputeVectorDotProduct(&Bloc_Matrix_A[index*NoofCols],
  284.                                                       Direction_Vector, VectorSize);
  285.         }
  286.         temp = ComputeVectorDotProduct(Bloc_Direction_Vector, Buffer,
  287.                                                  Bloc_VectorSize);
  288.  
  289.         MPI_Allreduce(&temp, &val, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  290.         Tau = Delta0 / val;
  291.  
  292.         /*Compute new vector Xnew = Xold + Tau*Direction........................*/
  293.         /*Compute BlocResidueVec  = BlocResidueVect + Tau*Bloc_MatA*DirVector...*/
  294.         for(index = 0; index<Bloc_VectorSize; index++){
  295.             Bloc_Vector_X[index]       = Vector_X[MyRank*Bloc_VectorSize + index] +
  296.                                                   Tau*Bloc_Direction_Vector[index];
  297.             Bloc_Residue_Vector[index] = Bloc_Residue_Vector[index] +
  298.                                                   Tau*Buffer[index];
  299.         }
  300.         /*...Gather New Vector X at all  processes......*/
  301.         MPI_Allgather(Bloc_Vector_X, Bloc_VectorSize, MPI_DOUBLE, Vector_X,
  302.                           Bloc_VectorSize, MPI_DOUBLE, MPI_COMM_WORLD);
  303.  
  304.         SolvePrecondMatrix(Bloc_Precond_Matrix, Bloc_HVector, Bloc_Residue_Vector,
  305.                                  Bloc_VectorSize);
  306.         Bloc_Delta1 = ComputeVectorDotProduct(Bloc_Residue_Vector, Bloc_HVector,
  307.                                                           Bloc_VectorSize);
  308.        
  309.         MPI_Allreduce(&Bloc_Delta1, &Delta1, 1, MPI_DOUBLE, MPI_SUM,
  310.                           MPI_COMM_WORLD);
  311.  
  312.         if(Delta1 < EPSILON)
  313.             break;
  314.  
  315.         Beta   = Delta1 / Delta0;
  316.         Delta0 = Delta1;
  317.         for(index=0; index<Bloc_VectorSize; index++){
  318.             Bloc_Direction_Vector[index] = -Bloc_HVector[index] +
  319.                                                      Beta*Bloc_Direction_Vector[index];
  320.         }
  321.     }while(Delta0 > EPSILON && Iteration < MAX_ITERATIONS);
  322.  
  323.     MPI_Barrier(MPI_COMM_WORLD);
  324.     EndTime = MPI_Wtime();
  325.  
  326.   if (MyRank == 0) {
  327.  
  328.      printf ("\n");
  329.      printf(" ------------------------------------------- \n");
  330.      printf("Results of Jacobi Method on processor %d: \n", MyRank);
  331.      printf ("\n");
  332.  
  333.      printf("Matrix Input_A \n");
  334.      printf ("\n");
  335.      for (irow = 0; irow < n_size; irow++) {
  336.         for (icol = 0; icol < n_size; icol++)
  337.             printf("%.3lf  ", Matrix_A[irow][icol]);
  338.         printf("\n");
  339.      }
  340.      printf ("\n");
  341.      printf("Matrix Input_B \n");
  342.      printf("\n");
  343.      for (irow = 0; irow < n_size; irow++) {
  344.          printf("%.3lf\n", Input_B[irow]);
  345.      }
  346.      printf ("\n");
  347.      printf("Solution vector \n");
  348.       printf("Number of iterations = %d\n",Iteration);
  349.      printf ("\n");
  350.      for(irow = 0; irow < n_size; irow++)
  351.         printf("%.12lf\n",Vector_X[irow]);
  352.      printf(" --------------------------------------------------- \n");
  353.   }
  354.  
  355.  
  356.  
  357.     MPI_Finalize();
  358. }
RAW Paste Data
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
Top