#include <stdio.h>

#include <stdlib.h>

#include <windows.h>

#include <math.h>

#include <time.h>

#define _CRT_SECURE_NO_WARNINGS_

#define B(i,j) MM[i*(n+2)+j]

#define A(i,j) M[i*(n+2)+j]

float *M, *MM, eps, epst, t;

int n, iter, it;

int iter_Up();

void var_M();

void Ax_b();

int main()

{

	int i, j, code, fl = 0, c;

	int fl1 = 0, iopt;

	float topt, dt = 0.05;

	char s[2];

	SetConsoleCP(1251); //чтобы воспринимала русские буквы

	SetConsoleOutputCP(1251);

	printf(" *** Программа находит решение системы методом верхней релаксации *** ");

	printf("\nВведите размерность системы: ");

	scanf("%d", &n);

	double *M = (double*)malloc(n*(n+2));

	double *MM = (double*)malloc(n*(n + 2));

	printf("Введите 0/1 - задать систему случайным образом или вручную ");

	scanf("%d", &c);

	switch (c) {

	case 0:

		//инициализация датчика п.с. чисел текущим временем

		srand(time(NULL));

		for (i = 0; i < n; i++)for (j = 0; j < n + 1; j++)B(i, j) = 0.5 - rand() / (RAND_MAX + 1.0);

		break;

	default:

		for (i = 0; i < n; i++) {

			for (j = 0; j < n; j++) {

				printf("Введите B(%d,%d): ", i, j);

				scanf("%f", &B(i, j));

			}

			printf("Введите b(%d): ", i);

			scanf("%f", &B(i, n));

		}

	}

	printf("Введенная система:\n");

	for (i = 0; i < n; i++) {

		for (j = 0; j < n; j++) {

			printf("%f ", B(i, j));

		}

		printf("%f\n", B(i, n));

	}

	var_M(); //умножаем систему на At

	printf("Преобразованная система:\n");

	for (i = 0; i < n; i++) {

		for (j = 0; j < n; j++) {

			printf("%f ", A(i, j));

		}

		printf("%f\n", A(i, n));

	}

	while (fl == 0) {

		printf("Введите число итераций,точность и параметр t: ");

		scanf("%d %f %f", &iter, &eps, &t);

		//printf("Введенные точность, число итераций и t = %f %d %.2f",eps,iter,t);

		var_M();

		code = iter_Up();

		switch (code)

		{

		case 0: //все хорошо

			printf("\nКод = %d, Число использованных итераций = %d, Точность = %.8f", code, it, epst);

			printf("\nРешение:\n");

			for (i = 0; i < n; i++) {

				printf("%f ", A(i, n + 1));

			}

			Ax_b();

			printf("\nНевязка:\n");

			for (i = 0; i < n; i++) {

				printf("%f ", B(i, n + 1));

			}

			printf("\nНайти отимальное t -> Введите 0, иначе 1 - ввод новых параметров, 2 - выход: "); scanf("%d", &fl1);

			if (fl1 == 1) break;

			if (fl1 == 2) { fl = 1; break; }

			if (fl1 == 0) { //нужно найти оптимальное t

				iopt = iter; t = 0;

				for (i = 0; i < 40; i++) {

					t = t + dt; // dt = 0.05

					var_M();

					code = iter_Up();

					if ((code == 0) && (it < iopt)) { iopt = it; topt = t; }

				}

				printf("\nПолучены следующие оптимальные значения t и iter: %.3f; %d", topt, iopt);

				printf("\nПовторить вычисления для новых параметров итерационного процесса - введите 0, иначе 1: ");

				scanf("%d", &fl);

			}

			break;

		case 1:

			printf("\nВ матрице системы есть диагональный компонент = 0 !");

			fl = 1;

			break;

		case 2:

			printf("\nМетод разошелся !");

			// fl=1;

			break;

		}

	}

	printf("\nДля завершения программы нажмите любую клавишу и ENTER ");

	scanf("%s", s); printf("\n");

	return 0;

}

//п/п реализует итерационный метод верхней релаксации

int iter_Up() {

	int i, j, k, fl = 0;

	float nt, buf, norm, norm_old = 0, t1 = 1 - t;

	//нормировка системы

	for (i = 0; i < n; i++) {

		A(i, n + 1) = A(i, n);

		if (fabs(A(i, i) < 1.E-30))return 1; //есть A(i,i)=0

		A(i, n) = t * A(i, n) / A(i, i);

		for (j = 0; j < n; j++)if (i != j)A(i, j) = t * A(i, j) / A(i, i);

	}

	it = 0; //число выполненных итераций

	for (k = 1; k <= iter; k++) { //цикл по итерациям

		it = it + 1; norm = 0;

		for (i = 0; i < n; i++) { //цикл по уравнениям

			buf = t1 * A(i, n + 1) + A(i, n);

			for (j = 0; j < n; j++) if (i != j)buf = buf - A(i, j)*A(j, n + 1);

			nt = fabs(buf - A(i, n + 1)); A(i, n + 1) = buf;

			if (nt > norm) norm = nt;

		}

		epst = norm;

		if (norm <= eps)return 0; //все хорошо

		if (norm > norm_old)fl = fl + 1; else fl = 0;

		if (fl > 3)return 2; //процесс расходится

		norm_old = norm;

	}

	return 0;

}

//п/п симметрезует систему

void var_M() {

	int i, j, k;

	for (i = 0; i < n; i++) {

		for (j = 0; j <= n; j++) {

			A(i, j) = 0.;

			for (k = 0; k < n; k++)
				A(i, j) = A(i, j) + B(k, i)*B(k, j);

		}

	}

}

//п/п вычисляет невязку

void Ax_b() {

	int i, j;

	for (i = 0; i < n; i++) {

		B(i, n + 1) = B(i, n);

		for (j = 0; j < n; j++) B(i, n + 1) = B(i, n + 1) - B(i, j)*A(j, n + 1);

	}

}