• API
• FAQ
• Tools
• Archive
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. ************************************************************************
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;
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.
Top