Hello World!

Author: minar09

Diff for: MPI_python.py

@@ -0,0 +1,4 @@
+from mpi4py import MPI
+
+comm = MPI.COMM_WORLD
+print("%d of %d" % (comm.Get_rank(), comm.Get_size()))

Diff for: OMP_Fortran.f90

@@ -0,0 +1,115 @@
+  Program Main
+!====================================================================
+!  eigenvalues and eigenvectors of a real symmetric matrix
+!  Method: calls Jacobi
+!====================================================================
+implicit none
+integer, parameter :: n=3
+double precision :: a(n,n), x(n,n)
+double precision, parameter:: abserr=1.0e-09
+integer i, j
+
+! matrix A
+  data (a(1,i), i=1,3) /   1.0,  2.0,  3.0 /
+  data (a(2,i), i=1,3) /   2.0,  2.0, -2.0 /
+  data (a(3,i), i=1,3) /   3.0, -2.0,  4.0 /
+
+! print a header and the original matrix
+  write (*,200)
+  do i=1,n
+     write (*,201) (a(i,j),j=1,n)
+  end do
+
+  call Jacobi(a,x,abserr,n)
+
+! print solutions
+  write (*,202)
+  write (*,201) (a(i,i),i=1,n)
+  write (*,203)
+  do i = 1,n
+     write (*,201)  (x(i,j),j=1,n)
+  end do
+
+200 format (' Eigenvalues and eigenvectors (Jacobi method) ',/, &
+            ' Matrix A')
+201 format (6f12.6)
+202 format (/,' Eigenvalues')
+203 format (/,' Eigenvectors')
+end program main
+
+subroutine Jacobi(a,x,abserr,n)
+!===========================================================
+! Evaluate eigenvalues and eigenvectors
+! of a real symmetric matrix a(n,n): a*x = lambda*x 
+! method: Jacoby method for symmetric matrices 
+! Alex G. (December 2009)
+!-----------------------------------------------------------
+! input ...
+! a(n,n) - array of coefficients for matrix A
+! n      - number of equations
+! abserr - abs tolerance [sum of (off-diagonal elements)^2]
+! output ...
+! a(i,i) - eigenvalues
+! x(i,j) - eigenvectors
+! comments ...
+!===========================================================
+implicit none
+integer i, j, k, n
+double precision a(n,n),x(n,n)
+double precision abserr, b2, bar
+double precision beta, coeff, c, s, cs, sc
+
+! initialize x(i,j)=0, x(i,i)=1
+! *** the array operation x=0.0 is specific for Fortran 90/95
+x = 0.0
+do i=1,n
+  x(i,i) = 1.0
+end do
+
+! find the sum of all off-diagonal elements (squared)
+b2 = 0.0
+do i=1,n
+  do j=1,n
+    if (i.ne.j) b2 = b2 + a(i,j)**2
+  end do
+end do
+
+if (b2 <= abserr) return
+
+! average for off-diagonal elements /2
+bar = 0.5*b2/float(n*n)
+
+do while (b2.gt.abserr)
+  do i=1,n-1
+    do j=i+1,n
+      if (a(j,i)**2 <= bar) cycle  ! do not touch small elements
+      b2 = b2 - 2.0*a(j,i)**2
+      bar = 0.5*b2/float(n*n)
+! calculate coefficient c and s for Givens matrix
+      beta = (a(j,j)-a(i,i))/(2.0*a(j,i))
+      coeff = 0.5*beta/sqrt(1.0+beta**2)
+      s = sqrt(max(0.5+coeff,0.0))
+      c = sqrt(max(0.5-coeff,0.0))
+! recalculate rows i and j
+      do k=1,n
+        cs =  c*a(i,k)+s*a(j,k)
+        sc = -s*a(i,k)+c*a(j,k)
+        a(i,k) = cs
+        a(j,k) = sc
+      end do
+! new matrix a_{k+1} from a_{k}, and eigenvectors 
+      do k=1,n
+        cs =  c*a(k,i)+s*a(k,j)
+        sc = -s*a(k,i)+c*a(k,j)
+        a(k,i) = cs
+        a(k,j) = sc
+        cs =  c*x(k,i)+s*x(k,j)
+        sc = -s*x(k,i)+c*x(k,j)
+        x(k,i) = cs
+        x(k,j) = sc
+      end do
+    end do
+  end do
+end do
+return
+end

Diff for: OpenMP_C.c

@@ -0,0 +1,17 @@
+#include<stdio.h>
+
+int main( int ac, char **av)
+
+{
+
+	#pragma omp parallel // specify the code between the curly brackets is part of an OpenMP 	parallel section.
+
+	{
+
+		printf("Hello World!!!\n");
+
+	}
+
+	return 0;
+
+}

Diff for: gauss_seidel.c

@@ -0,0 +1,113 @@
+/* L-20 MCS 572 Fri 7 Oct 2016 : gauss_seidel.c 
+ * This program runs the method of Gauss-Seidel on a test system A*x = b,
+ * where A is diagonally dominant and the exact solution consists
+ * of all ones.  The dimension can be supplied at the command line.
+ * This program is an adaption of the method of Jacobi of L-19,
+ * observe the improved convergence rate of Gauss-Seidel. */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+void test_system
+ ( int n, double **A, double *b );
+/*
+ * Given n on entry,
+ * On return is an n-by-n matrix A
+ * with n+1 on the diagonal and 1 elsewhere.
+ * The elements of the right hand side b
+ * all equal 2*n, so the exact solution x
+ * to A*x = b is a vector of ones. */
+
+void run_gauss_seidel_method
+ ( int n, double **A, double *b,
+   double epsilon, int maxit,
+   int *numit, double *x );
+/*
+ * Runs the  method of Gauss-Seidel for A*x = b.
+ *
+ * ON ENTRY :
+ *   n        the dimension of the system;
+ *   A        an n-by-n matrix A[i][i] /= 0;
+ *   b        an n-dimensional vector;
+ *   epsilon  accuracy requirement;
+ *   maxit    maximal number of iterations;
+ *   x        start vector for the iteration.
+ *
+ * ON RETURN :
+ *   numit    number of iterations used;
+ *   x        approximate solution to A*x = b. */
+
+int main ( int argc, char *argv[] )
+{
+   int n,i;
+   if(argc > 1)
+      n = atoi(argv[1]);
+   else
+   {
+      printf("give the dimension : ");
+      scanf("%d",&n);
+   }
+   {
+      double *b;
+      b = (double*) calloc(n,sizeof(double));
+      double **A;
+      A = (double**) calloc(n,sizeof(double*));
+      for(i=0; i<n; i++)
+         A[i] = (double*) calloc(n,sizeof(double));
+      test_system(n,A,b);
+      double *x;
+      x = (double*) calloc(n,sizeof(double));
+      /* we start at an array of all zeroes */
+      for(i=0; i<n; i++) x[i] = 0.0;
+      double eps = 1.0e-4;
+      int maxit = 2*n*n;
+      int cnt = 0;
+      run_gauss_seidel_method(n,A,b,eps,maxit,&cnt,x);
+      printf("computed %d iterations\n",cnt);
+      double sum = 0.0;
+      for(i=0; i<n; i++) /* compute the error */
+      {
+         double d = x[i] - 1.0;
+         sum += (d >= 0.0) ? d : -d;
+      }
+      printf("error : %.3e\n",sum);
+   }
+   return 0;
+}
+
+void test_system
+ ( int n, double **A, double *b )
+{
+   int i,j;
+   for(i=0; i<n; i++)
+   {
+      b[i] = 2.0*n;
+      for(j=0; j<n; j++) A[i][j] = 1.0;
+      A[i][i] = n + 1.0;
+   }
+}
+
+void run_gauss_seidel_method
+ ( int n, double **A, double *b,
+   double epsilon, int maxit,
+   int *numit, double *x )
+{
+   double *dx = (double*) calloc(n,sizeof(double));
+   int i,j,k;
+
+   for(k=0; k<maxit; k++)
+   {
+      double sum = 0.0;
+      for(i=0; i<n; i++)
+      {
+         dx[i] = b[i];
+         for(j=0; j<n; j++)
+            dx[i] -= A[i][j]*x[j]; 
+         dx[i] /= A[i][i]; x[i] += dx[i];
+         sum += ( (dx[i] >= 0.0) ? dx[i] : -dx[i]);
+      }
+      printf("%4d : %.3e\n",k,sum);
+      if(sum <= epsilon) break;
+   }
+   *numit = k+1; free(dx);
+}

Diff for: gaussseidel.py

@@ -0,0 +1,40 @@
+def gaussSeidel(A, b, x, N, tol):
+    maxIterations = 1000000
+    xprev = [0.0 for i in range(N)]
+    for i in range(maxIterations):
+        for j in range(N):
+            xprev[j] = x[j]
+        for j in range(N):
+            summ = 0.0
+            for k in range(N):
+                if (k != j):
+                    summ = summ + A[j][k] * x[k]
+            x[j] = (b[j] - summ) / A[j][j]
+        diff1norm = 0.0
+        oldnorm = 0.0
+        for j in range(N):
+            diff1norm = diff1norm + abs(x[j] - xprev[j])
+            oldnorm = oldnorm + abs(xprev[j])  
+        if oldnorm == 0.0:
+            oldnorm = 1.0
+        norm = diff1norm / oldnorm
+        if (norm < tol) and i != 0:
+            print("Sequence converges to [", end="")
+            for j in range(N - 1):
+                print(x[j], ",", end="")
+            print(x[N - 1], "]. Took", i + 1, "iterations.")
+            return
+    print("Doesn't converge.")
+
+matrix2 = [[3.0, 1.0], [2.0, 6.0]]
+vector2 = [5.0, 9.0]
+guess = [0.0, 0.0]
+
+matrix3 = [[9.0, -3.0], [-2.0, 8.0]]
+vector3 = [6.0, -4.0]
+
+matrix4 = [[1.0, -3.0], [-2.0, 8.0]]
+
+gaussSeidel(matrix2, vector2, guess, 2, 0.00000000000001)
+gaussSeidel(matrix3, vector3, guess, 2, 0.00000000000001)
+gaussSeidel(matrix4, vector3, guess, 2, 0.00000000000001)

Diff for: gs_seq_mpi.c

@@ -0,0 +1,93 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <math.h>
+#include "mpi.h"
+
+#define MAX_ITER 100
+#define MAX 100 //maximum value of the matrix element
+#define TOL 0.000001
+
+// Generate a random float number with the maximum value of max
+float rand_float(int max){
+  return ((float)rand()/(float)(RAND_MAX)) * max;
+}
+
+
+// Allocate 2D matrix
+void allocate_init_2Dmatrix(float ***mat, int n, int m){
+  int i, j;
+  *mat = (float **) malloc(n * sizeof(float *));
+  for(i = 0; i < n; i++) {
+    (*mat)[i] = (float *)malloc(m * sizeof(float));
+    for (j = 0; j < m; j++)
+      (*mat)[i][j] = rand_float(MAX);
+  }
+} 
+
+// solver
+void solver(float ***mat, int n, int m){
+  float diff = 0, temp;
+  int done = 0, cnt_iter = 0, i, j;
+
+  while (!done && (cnt_iter < MAX_ITER)){
+    diff = 0;
+    for (i = 1; i < n - 1; i++)
+      for (j = 1; j < m - 1; j++){
+	temp = (*mat)[i][j];
+	(*mat)[i][j] = 0.2 * ((*mat)[i][j] + (*mat)[i][j - 1] + (*mat)[i - 1][j] + (*mat)[i][j + 1] + (*mat)[i + 1][j]);
+	diff += abs((*mat)[i][j] - temp);
+      }
+    if (diff/n/n < TOL)
+      done = 1; 
+    cnt_iter ++;
+  }
+
+  if (done)
+    printf("Solver converged after %d iterations\n", cnt_iter);
+  else
+    printf("Solver not converged after %d iterations\n", cnt_iter);
+  
+}
+
+int main(int argc, char *argv[]) {
+  int n, communication;
+  float **a;
+  
+  MPI_Init(&argc, &argv);
+
+  if (argc < 3) {
+    printf("Call this program with two parameters: matrix_size communication \n");
+    printf("\t matrix_size: Add 2 to a power of 2 (e.g. : 18, 1026)\n");
+    printf("\t communication:\n");
+    printf("\t\t 0: initial and final using point-to-point communication\n");
+    printf("\t\t 1: initial and final using collective communication\n");
+    
+    exit(1);
+  }
+
+  n = atoi(argv[1]);
+  communication =  atoi(argv[2]);
+  printf("Matrix size = %d communication = %d\n", n, communication);
+
+  allocate_init_2Dmatrix(&a, n, n);
+
+  double tsop = MPI_Wtime();
+  solver(&a, n, n);
+  double tfop = MPI_Wtime();
+
+  FILE *f;
+		if (access("seq.csv", F_OK) == -1) {
+ 			f = fopen("seq.csv", "a");
+			fprintf(f, "Operations-time\n");
+		}
+		else {
+			f = fopen("seq.csv", "a");
+		}
+
+		fprintf(f, "%f;\n", tfop - tsop);
+		fclose(f);
+
+  MPI_Finalize();
+  return 0;
+}