pp-48.c

// 48. Целочисленное линейное программирование. Заданы целочисленная матрица A,
//     целочисленные векторы a,b. Нужно найти вектор x из нулей и единиц такой,
//     что Ax<b и велична ax максимальна.

#include <assert.h>
#include <math.h>
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

static int *Alloc_vector(int size) {
	int *v;

	if ((v = calloc(size, sizeof(*v))) == NULL) {
		fprintf(stderr, "Couldn't allocate %d bytes\n", size * sizeof(*v));
		exit(1);
	}

	return v;
}

static void Free_vector(int *v) {
	assert(v != NULL);
	free(v);
}

static void Read_vector(FILE *stream, int *v, int size) {
	int i;

	for (i = 0; i < size; i++) {
		fscanf(stream, "%d", &v[i]);
	}
}

static void Write_vector(FILE *stream, int *v, int size) {
	int i;

	for (i = 0; i < size; i++) {
		if (i > 0) {
			fprintf(stream, " ");
		}
		fprintf(stream, "%d", v[i]);
	}
	fprintf(stream, "\n");
}

static void Write_matrix(FILE *stream, int *v, int nrows, int ncols) {
	int i, n;

	for (i = 0, n = nrows * ncols; i < n; i += ncols) {
		Write_vector(stream, &v[i], ncols);
	}
	fprintf(stream, "\n");
}

static int Multiply_vectors(int *v1, int *v2, int size) {
	int i, p;

	for (i = 0, p = 0; i < size; i++) {
		p += v1[i] * v2[i];
	}

	return p;
}

static double Get_vector_length(int *v, int size) {
	int i, acc;

	for (i = 0, acc = 0; i < size; i++) {
		acc += v[i] * v[i];
	}

	return sqrt(acc);
}

static int Compare_vectors(int *v1, int *v2, int size) {
	double v1_len, v2_len;

	v1_len = Get_vector_length(v1, size);
	v2_len = Get_vector_length(v2, size);

	if (v1_len < v2_len) {
		return -1;
	}
	if (v1_len > v2_len) {
		return 1;
	}

	return 0;
}

static int Compare_vectors_element_wise(int *v1, int *v2, int size) {
	int i;

	for (i = 0; i < size; i++) {
		if (v1[i] < v2[i]) {
			return -1;
		}
		if (v1[i] > v2[i]) {
			return 1;
		}
	}

	return 0;
}

static void Generate_vector(int *v, int size, int min, int max) {
	int i;

	for (i = 0; i < size; i++) {
		v[i] = min + rand() % (max - min);
	}
}

static int Next_x(int *x, int size, int base) {
	int i, j;

	for (i = 0; i < size; i++) {
		if (x[i] < base - 1) {
			x[i]++;
			for (j = i - 1; j >= 0; j--) {
				x[j] = 0;
			}
			return 1;
		}
	}

	return 0;
}

int main(int argc, char **argv) {
	int root = 0;                   /* ранг главного процесса */
	MPI_Comm comm = MPI_COMM_WORLD; /* коммуникатор */
	int rank;                       /* ранг данного процесса */ 
	int nprocs;                     /* общее количество процессов */
	int M;                          /* количество строк в матрице A */
	int N;                          /* число столбцов в оной же (и длина остальных векторов) */
	int *A, *a, *b, *x;             /* исходные данные... */
	int *r;                         /* кусок матрицы A, над которым парится данный процесс */
	int nrows, count;               /* количество строк и элементов из A, отведенных этому процессу */
	int *Ax, *rx;                   /* произведение A на x и r на x */
	int ax, ax_max = 0;             /* произведение векторов a и x и его максимальное значение */
	int *x_best;                    /* вектор x, соответствующий ax_max */
	int *counts, *displs;           /* для MPI_Scatterv() и MPI_Gatherv() */
	int i, j;                       /* ну это понятно чё */
	double t_start;

	MPI_Init(&argc, &argv);

	MPI_Comm_rank(comm, &rank);
	MPI_Comm_size(comm, &nprocs);

	t_start = MPI_Wtime();

	if (rank == root) {
		/* Размер матрицы и т.п. берется из командной строки. */
		if (argc < 3) {
			fprintf(stderr, 
				"Usage: %s <M> <N> <file>\n"
				" <M> - the number of rows in A\n"
				" <N> - the number of columns in A (and the length of a, b and x)\n"
				" <file> - an input file to read data from\n"
			, argv[0]);
			MPI_Abort(comm, 1);
		}

		M = atoi(argv[1]);
		N = atoi(argv[2]);
	}

	/* Разсылаем всем процессам значения M и N т.к. они пока только у мастера. */
	MPI_Bcast(&M, 1, MPI_INT, root, comm);
	MPI_Bcast(&N, 1, MPI_INT, root, comm);

	nrows = M / nprocs;
	if (rank == root) {
		nrows += M % nprocs;
	}
	count = nrows * M;

	/* Выделение памяти под матрицу и вектора. */
	a = Alloc_vector(N);
	b = Alloc_vector(N);
	x = Alloc_vector(N);
	r = Alloc_vector(N * nrows);
	rx = Alloc_vector(N);

	if (rank == root) {
		A = Alloc_vector(N * M);
		Ax = Alloc_vector(N);
		x_best = Alloc_vector(N);
	} else {
		A = Ax = x_best = 0;
	}

	if (rank == root) {
		/* Если задано имя файла, данные считываются из него, если нет - генерируются случайно. */
		if (argc >= 4) {
			FILE *fp;

			if ((fp = fopen(argv[3], "r")) == NULL) {
				fprintf(stderr, "Couldn't read from file '%s'\n", argv[3]);
				MPI_Abort(comm, 2);
			}

			Read_vector(fp, a, N);
			Read_vector(fp, b, N);
			Read_vector(fp, A, N * M);
		} else {
			srand((unsigned)time(NULL) * rank);

			Generate_vector(a, N, -10, 10);
			Generate_vector(b, N, -100, 100);
			Generate_vector(A, N * M, -10, 10);
		}

		/*printf("a:\n");
		Write_vector(stdout, a, N);

		printf("b:\n");
		Write_vector(stdout, b, N);

		printf("A:\n");
		Write_matrix(stdout, A, M, N); */
	}

	counts = Alloc_vector(nprocs);
	displs = Alloc_vector(nprocs);

	MPI_Gather(&count, 1, MPI_INT, counts, 1, MPI_INT, root, comm);

	for (i = 0, j = 0; i < nprocs; i++) {
		displs[i] = j;
		j += counts[i];
	}

	/* Разсылаем матрицу A по кускам всем процессам (собственно, для этого
	 * и забивали массивы counts и displs чуть выше. */
	MPI_Scatterv(A, counts, displs, MPI_INT, r, count, MPI_INT, root, comm);

	/* Оба вектора a и b отсылаются целиком. */
	MPI_Bcast(a, N, MPI_INT, root, comm);
	MPI_Bcast(b, N, MPI_INT, root, comm);

	do {
		/*if (rank == root) {
			printf("x\n");
			Write_vector(stdout, x, N);
		}*/

		/* Передаем всем вектор x. */
		MPI_Bcast(x, N, MPI_INT, root, comm);

		/* Перемножаем свой кусок матрицы с вектором x. */
		for (i = 0; i < nrows; i++) {
			for (j = 0; j < M; j++) {
				rx[i] = Multiply_vectors(&r[i * M], x, M);
			}
		}

		MPI_Gather(&nrows, 1, MPI_INT, counts, 1, MPI_INT, root, comm);

		for (i = 0, j = 0; i < nprocs; i++) {
			displs[i] = j;
			j += counts[i];
		}

		/* Отослать главарю свой rx, чтобы тот мог собрать вектор Ax. */
		MPI_Gatherv(rx, nrows, MPI_INT, Ax, counts, displs, MPI_INT, root, comm);

		if (rank == root) {
			/*printf("Ax:\n");
			Write_vector(stdout, Ax, N);*/

			/* Короче тут сравнивается Ax и b. Compare_vectors() возвращает -1 если первый вектор
			 * меньше второго (по длине). Если выяснится, что вектора надо сравнивать поэлементно,
			 * замени на Compare_vectors_element_wise() */
			if (Compare_vectors(Ax, b, N) < 0) {
				ax = Multiply_vectors(a, x, N);
				if (ax > ax_max) {
					ax_max = ax;
					/* О! Новый максимум - надо сохранить x. */
					memcpy(x_best, x, N * sizeof(*x));
				}
			}
		}
	} while (Next_x(x, N, 2)); /* цикл выполняется до тех пор, пока мы не перебрали все значения x. */

	/* Вывести время выполнения каждого процесса. */
	printf("[%d] Elapsed time: %lf\n", rank, MPI_Wtime() - t_start);

	if (rank == root) {
		printf("\nThe maximum of ax was:\n  %d\nfor the following x:\n  ", ax_max);
		Write_vector(stdout, x_best, N);
	}

	if (rank == root) {
		Free_vector(A);
		Free_vector(Ax);
		Free_vector(x_best);
	}

	Free_vector(a);
	Free_vector(b);
	Free_vector(x);
	Free_vector(r);
	Free_vector(rx);

	Free_vector(counts);
	Free_vector(displs);

	return MPI_Finalize();
}