lab2_schedule

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

#define N 1024

const double t = 1e-5;
const double e = 1e-9;

#define schedule_type static
#define schedule_arg 1

typedef struct matrix {
	double* data;
	int w, h;
} matrix;

void setInitData(matrix *A, matrix *B, matrix *x) {
	A->data = (double*)malloc(N * N * sizeof(double));
	A->w = N;
	A->h = N;

	B->data = (double*)malloc(N * sizeof(double));
	B->w = 1;
	B->h = N;

	x->data = (double*)malloc(N * sizeof(double));
	x->w = 1;
	x->h = N;

	for (int i = 0; i < A->w; i++) {
		for (int j = 0; j < A->h; j++) {
			A->data[i * N + j] = 1.0;
			if (i == j) A->data[i * N + j] = 2.0;
		}
	}

	for (int i = 0; i < B->h; i++) {
		B->data[i] = N + 1;
		x->data[i] = 0;
	}
}

void mult(matrix *a, matrix *b, matrix *res) {
	double temp = 0;

#pragma omp parallel for reduction(+:temp) schedule(schedule_type, schedule_arg)
	for (int i = 0; i < a->h; i++) {
		temp = 0;

		for (int k = 0; k < b->h; k++)
			temp += a->data[(i * a->w) + k] * b->data[k];

		res->data[i] = temp;
	}
}

void sub(matrix *a, matrix *b, matrix *res) {
#pragma omp parallel for schedule(schedule_type, schedule_arg)
	for (int i = 0; i < a->h; i++)
		res->data[i] = a->data[i] - b->data[i];
}

void sum(matrix *a, matrix *b, matrix *res) {
#pragma omp parallel for schedule(schedule_type, schedule_arg)
	for (int i = 0; i < a->h; i++)
		res->data[i] = a->data[i] + b->data[i];
}

void scalarMult(matrix *a, double k, matrix *res) {
#pragma omp parallel for schedule(schedule_type, schedule_arg)
	for (int i = 0; i < a->h; i++)
		res->data[i] = a->data[i] * k;
}

double squaredNorma(matrix *a) {
	double res = 0;

#pragma omp parallel for reduction(+:res) schedule(schedule_type, schedule_arg)
	for (int i = 0; i < a->h; i++)
		res += a->data[i] * a->data[i];

	return res;
}


int main(int argc, char** argv) {
	if (argc > 1) {
		omp_set_num_threads(atoi(argv[1]));
	}
	else {
		printf("Error: specify number of threads!\n");
		return 0;
	}

	double timeStart, Ax_b_norma = N, b_norma;
	int counter = 0;
	matrix A, B, x, temp;

	setInitData(&A, &B, &x);
	temp.data = (double*)malloc(N * sizeof(double));
	temp.h = N;
	temp.w = 1;
	b_norma = sqrt(squaredNorma(&B));

	timeStart = omp_get_wtime();

	while (sqrt(Ax_b_norma) / b_norma >= e) {
		mult(&A, &x, &temp);
		sub(&temp, &B, &temp);
		Ax_b_norma = squaredNorma(&temp);
		scalarMult(&temp, t, &temp);
		sub(&x, &temp, &x);
		counter++;
	}

	printf("x[0]: %lf\nCounter: %d\nTime: %lf s\n", x.data[0], counter, omp_get_wtime() - timeStart);

	free(A.data);
	free(B.data);
	free(x.data);
	free(temp.data);

	return 0;
}