Added elliptic equation

.gitignore

+compiled
+*.o
+*.pyc
+*.out

README.md

@@ -6,13 +6,14 @@
 
 Где
 * **-w** - ожидание резльтатов(может быть долго, если много потоков, и лучше его не выставлять)
-* **-t matrix/thermal/oscillation** - задание которое запускать
+* **-t matrix/thermal/oscillation/elliptic** - задание которое запускать
 * **-np число_потоков**
 * **-o вариант** - вариант по методичке(нумерация от 1). Для матрицы не актуально
 * Варинаты аргументов для разных задач:
     * **-n размер_матрицы** - для матриц
-    * **-N число_частей_по_оси_X -t_therm временной_отрезок** - для уравннения теплопроводности
-    * **-t_osc временной_отрезок -dt шаг_по_времени -dx шаг_по_X** - для уравннения колебания струны
+    * **-N число_частей_по_оси_X -t_therm временной_отрезок** - для уравнения теплопроводности
+    * **-t_osc временной_отрезок -dt_osc шаг_по_времени -dx_osc шаг_по_X** - для уравнения колебания струны
+    * **-dx_osc шаг_по_X -dy_osc шаг_по_Y** - для эллиптического уравнения
 
 ### Справка скрипта ###
 

elliptic.cc

+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+
+#include "mpi.h"
+
+const int MASTER_RANK = 0;
+const int MASTER_SEND_DATA = 1;
+const int MASTER_RECEIVE_DATA = 2;
+
+double x_y_min(int option)
+{
+	switch(option)
+	{
+		case 1: return 0.0;
+		case 2: return 0.0;
+		case 3: return -1.0;
+		case 4: return 0.0;
+		case 5: return 0.0;
+		case 6: return 0.0;
+		default: return .0;
+	}
+}
+
+double x_y_max(int option)
+{
+	switch(option)
+	{
+		case 1: return 1.0;
+		case 2: return 1.0;
+		case 3: return 1.0;
+		case 4: return 2.0;
+		case 5: return 1.0;
+		case 6: return 1.0;
+		default: return .0;
+	}
+}
+
+
+double u_x_min(double x, int option)
+{
+	switch(option)
+	{
+		case 1: return .0;
+		case 2: return x + x_y_min(option);
+		case 3: return .0;
+		case 4: return x;
+		case 5: return .0;
+		case 6: return x;
+		default: return .0;
+	}
+}
+
+double u_x_max(double x, int option)
+{
+	switch(option)
+	{
+		case 1: return .0;
+		case 2: return x + x_y_max(option);
+		case 3: return .0;
+		case 4: return x;
+		case 5: return x;
+		case 6: return x;
+		default: return .0;
+	}
+}
+
+double u_min_y(double y, int option)
+{
+	switch(option)
+	{
+		case 1: return .0;
+		case 2: return y + x_y_min(option);
+		case 3: return .0;
+		case 4: return x_y_min(option);
+		case 5: return .0;
+		case 6: return x_y_min(option);
+		default: return .0;
+	}
+}
+
+double u_max_y(double y, int option)
+{
+	switch(option)
+	{
+		case 1: return sin(M_PI * y);
+		case 2: return y + x_y_max(option);
+		case 3: return .0;
+		case 4: return x_y_max(option);
+		case 5: return y;
+		case 6: return x_y_max(option);
+		default: return .0;
+	}
+}
+
+double matrix(int i, int j, int n, int m, double dx, double dy)
+{
+    if((i / n) == (j / m)) // diagonal block
+    {
+        if(i % n == j % m) // diagonal element
+            return 2.0 * dx * dx + 2.0 * dy * dy;
+        else if((i % n == j % m - 1) || (i % n - 1 == j % m)) // near diagonal element
+            return -(dy * dy);
+        else
+            return 0.0;
+    }
+    else if((i / n == j / m - 1) || (i / n - 1 == j / m)) // near diagonal block
+    {
+        if(i % n == j % m) // diagonal element
+            return -(dx * dx);
+        else
+            return 0.0;
+    }
+
+	return 0.0;
+}
+
+double b_value(int i, int n, int m, double dx, double dy, int option)
+{
+	double ret = .0;
+
+	if(i % n == 0) // first elemnt in part
+		ret += dy * dy * u_min_y(x_y_min(option) + (i / n + 1) * dy, option);
+	else if(i % n == n - 1) // last elemnt in part
+		ret += dy * dy * u_max_y(x_y_min(option) + (i / n + 1) * dy, option);
+
+	if(i / n == 0) // first part
+		ret += dx * dx * u_x_min(x_y_min(option) + (i % n + 1) * dx, option);
+	if(i / n == m - 1) // last part
+		ret += dx * dx * u_x_max(x_y_min(option) + (i % n + 1) * dx, option);
+
+	return ret;
+}
+
+void fillData(double *A, double *b, double *u, int n, int m, double dx, double dy, int option)
+{
+	for(int i = 0; i < n * m; ++i)
+	{
+		for(int j = 0; j < n * m; ++j)
+			A[i * (n * m) + j] = matrix(i, j, n, m, dx, dy);
+
+		b[i] = b_value(i, n, m, dx, dy, option);
+		u[i] = .0;
+	}
+}
+
+int main(int argc, char *argv[])
+{
+	int rank;
+	int size;
+	double dy, dx;
+	unsigned int n, m;
+
+	int option;
+
+	double *A;
+	double *b;
+    double *u;
+
+	if(argc > 3) {
+		option = atoi(argv[1]);
+		dx = atof(argv[2]);
+		dy = atof(argv[3]);
+	}
+	else {
+		exit(1);
+	}
+
+	double x_y_diff = x_y_max(option) - x_y_min(option);
+
+    n = (int) floor(x_y_diff / dx);
+	m = (int) floor(x_y_diff / dy);
+
+	//eclude bound values from equations
+	int N = n - 1;
+	int M = m - 1;
+	int ASize = N * M;
+
+	MPI_Init(&argc, &argv);
+	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+	MPI_Comm_size(MPI_COMM_WORLD, &size);
+
+	int thCount = size - 1;
+
+	//Master thread allocates data and send to worker threads
+	if (rank == MASTER_RANK) {
+		double start;
+		double end;
+
+		int stop = 0;
+
+		printf("dx: %.16lf\n", dx);
+		printf("dy: %.16lf\n", dy);
+		printf("X part num: %u\n", n);
+		printf("Y part num: %u\n", m);
+		printf("Thread count: %d\n", thCount);
+
+		A = new double[ASize * ASize];
+		b = new double[ASize];
+		u = new double[ASize];
+
+		start = MPI_Wtime();
+
+		fillData(A, b, u, N, M, dx, dy, option);
+		//Solve system by Gauss method
+		for(int step = 0; step < ASize; ++step)
+		{
+			double max = 0.0;
+			int maxRowIndex = step;
+			// find row with max firt element in submatrix
+			for(int i = step; i < ASize; ++i)
+			{
+				double absValue = fabs(A[i * ASize + step]);
+
+				if(absValue > max)
+				{
+					maxRowIndex = i;
+					max = absValue;
+				}
+			}
+			double tmp;
+			// swap row with max element to the current row
+			if(maxRowIndex != step)
+			{
+				for(int j = step; j < ASize; ++j)
+				{
+					tmp = A[step * ASize + j];
+					A[step * ASize + j] = A[maxRowIndex * ASize + j];
+					A[maxRowIndex * ASize + j] = tmp;
+				}
+
+				tmp = b[step];
+				b[step] = b[maxRowIndex];
+				b[maxRowIndex] = tmp;
+			}
+
+			double p = 1.0 / A[step * ASize + step];
+
+			for (int j = step; j < ASize; ++j)
+				A[step * ASize + j] *= p;
+
+			b[step] *= p;
+
+			int workRows = ASize - step;
+			unsigned int rowsPerTh = workRows / thCount;
+			unsigned int extraRowCount = workRows % thCount;
+			unsigned int workThreadCount = thCount;
+
+			if(rowsPerTh == 0)// rows < threads
+			{
+				rowsPerTh = 1;
+				workThreadCount = extraRowCount;
+				extraRowCount = 0;
+			}
+
+			int rows = 0;
+			int offset = 0;
+
+			for (int th = 1; th <= workThreadCount; ++th)
+			{
+				rows = (th <= extraRowCount) ? rowsPerTh + 1 : rowsPerTh;
+
+				MPI_Send(&stop, 1, MPI_INT, th, MASTER_SEND_DATA, MPI_COMM_WORLD);
+				MPI_Send(&rows, 1, MPI_INT, th, MASTER_SEND_DATA, MPI_COMM_WORLD);
+				MPI_Send(A + offset * ASize, rows * ASize, MPI_DOUBLE, th, MASTER_SEND_DATA, MPI_COMM_WORLD);
+				MPI_Send(b + offset, rows, MPI_DOUBLE, th, MASTER_SEND_DATA, MPI_COMM_WORLD);
+
+				offset = offset + rows;
+			}
+
+			offset = 0;
+
+			for (int th = 1; th <= workThreadCount; ++th)
+			{
+				rows = (th <= extraRowCount) ? rowsPerTh + 1 : rowsPerTh;
+
+				MPI_Recv(A + offset * ASize, rows * ASize, MPI_DOUBLE, th, MASTER_RECEIVE_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+				MPI_Recv(b + offset, rows, MPI_DOUBLE, th, MASTER_RECEIVE_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+				offset = offset + rows;
+			}
+		}
+
+		stop = 1;
+
+		for (int th = 1; th <= thCount; ++th)
+			MPI_Send(&stop, 1, MPI_INT, th, MASTER_SEND_DATA, MPI_COMM_WORLD);
+
+		//Gauss reverse steps
+		for (int i = ASize - 1; i >= 0; --i)
+	    {
+			double p = b[i];
+
+			for (int j = i + 1; j < ASize; ++j)
+				p -= u[j] * A[i * ASize + j];
+
+			u[i] = p / A[i * ASize + i];
+	    }
+
+		end = MPI_Wtime();
+
+		printf("Work time: %lf\n", end - start);
+
+		delete A;
+		delete b;
+		delete u;
+	}
+	else
+	{
+		for(int step = 0; step < ASize; ++step)
+		{
+			int stop;
+
+			MPI_Recv(&stop, 1, MPI_INT, MASTER_RANK, MASTER_SEND_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+			if(stop > 0)
+				break;
+
+			int rows = 0;
+
+			MPI_Recv(&rows, 1, MPI_INT, MASTER_RANK, MASTER_SEND_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+			A = new double[rows * ASize];
+			b = new double[rows];
+
+			MPI_Recv(A, rows * ASize, MPI_DOUBLE, MASTER_RANK, MASTER_SEND_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+			MPI_Recv(b, rows, MPI_DOUBLE, MASTER_RANK, MASTER_SEND_DATA, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+
+			//substruct first row of part from others
+			for (int i = 1; i < rows; ++i)
+		    {
+		        double p = A[i * ASize];
+
+		        for (int j = 0; j < ASize; ++j)
+					A[i * ASize + j] -= p * A[j];
+
+		        b[i] -= p * b[0];
+		    }
+
+			MPI_Send(A, rows * ASize, MPI_INT, MASTER_RANK, MASTER_RECEIVE_DATA, MPI_COMM_WORLD);
+			MPI_Send(b, rows, MPI_DOUBLE, MASTER_RANK, MASTER_RECEIVE_DATA, MPI_COMM_WORLD);
+
+			delete A;
+			delete b;
+		}
+	}
+
+	MPI_Finalize();
+}

run.py

@@ -11,7 +11,7 @@ def _parse_arguments(argv):
 
     parser = ArgumentParser(prog="run.py",
                             description='\nCompiles and run MPI tasks')
-    parser.add_argument("-t", "--task",dest="task", choices=['matrix','thermal','oscillation'], default='matrix',
+    parser.add_argument("-t", "--task",dest="task", choices=['matrix','thermal','oscillation', 'elliptic'], default='matrix',
                             help='Task to run', required=True)
 
     parser.add_argument("-w", "--wait",action="store_true", dest="wait",
@@ -34,9 +34,15 @@ def _parse_arguments(argv):
 
     group.add_argument("-t_osc", dest="T", type=float, help='Specify the T (time length)')
 
-    group.add_argument("-dt", dest="dt", type=float, help='Specify the dt')
+    group.add_argument("-dt_osc", dest="dt", type=float, help='Specify the dt')
 
-    group.add_argument("-dx", dest="dx", type=float, help='Specify the dx')
+    group.add_argument("-dx_osc", dest="dx", type=float, help='Specify the dx')
+
+    group = parser.add_argument_group('Elliptic equation options')
+
+    group.add_argument("-dx_ell", dest="dx", type=float, help='Specify the dx')
+
+    group.add_argument("-dy_ell", dest="dy", type=float, help='Specify the dy')
 
     return parser.parse_known_args(argv)
 
@@ -66,6 +72,8 @@ if __name__ == "__main__":
         run_args = '%d %d %lf' % (args.option, args.N, args.T)
     elif args.task == 'oscillation':
         run_args = '%d %lf %lf %lf' % (args.option, args.T, args.dx, args.dt)
+    elif args.task == 'elliptic':
+        run_args = '%d %lf %lf' % (args.option, args.dx, args.dt)
 
     if not args.wait:
         os.system('mpirun -np %d %s %s' % (args.threads + 1, run_file, run_args))