calc.cpp

// Файл реалізації підпрограм користувача

#include "calc.h"

using namespace std;

unsigned int kolIter;			// кількість ітерацій методу
unsigned int kolOperation;		// кількість операцій (+,-,*,/) методу
unsigned _int64 calc_time;	// час обчислення (в мілісекундах)
//double EPS;

//=============== підпрограма зчитування елементів вхідної матриці з файлу ==============
double** inpMatrixInFile(int n, const char *str) {// розмірність матриці, покажчик на масив символів
	stringstream streamRead(str);		// конвертувати char * у стрім (потік) (щоб потім з нього читати числа)
	int kolRow = 0, kolCol = 0;			//	змінні підрахунку прочитаних рядків і колонок
	double **matrix = new double*[n];	// створення нової матриці
	bool err = false;		// прапорець помилки (помилка відсутня)
	for (int row = 0; row < n && !err; row++) {	// цикл по N рядках і поки немає помилки 
		matrix[row] = new double[n];			// створення нового рядка матриці
		kolCol = 0;			// обнулення кількості колонок
		for (int col = 0; col < n && !err; col++) {		// цикл по N колонках і поки немає помилки
			double el = 0;	// ініціалізація змінної (прочитане значення елемента матриці) 
			streamRead >> el;		// прочитати із потоку значення елемента 
			err = (streamRead) ? false : true;	// якщо вже кінець рядка, а ще не всі елементи прочитано (по кількості в рядку матриці), підняти прапорець  
			if (streamRead) {	// якщо покажчик у стрімі не вказує на NULL, то	
				matrix[row][col] = el;	// записати в поточний елемент матриці прочитаний із потоку елемент
				kolCol++;	// збільшити кількість прочитаних у рядку елементів на 1
			}
			else err = true;	// інакше підняти прапорець
		}
		err = !(kolCol == n);	// якщо кількість прочитаних у рядку елементів !=N, підняти прапорець
		kolRow++;				// збільшити кількість прочитаних рядків на 1
	}
	if (!err && kolRow == n)	// якщо не виникало помилки і кількість прочитаних рядків ==N, то
		return matrix;		// повернути адресу на динамічну матрицю із прочитаними елементами (прочитано вдало)
	else
		return NULL;		// повернути NULL (помилка читання даних)
}

// ============= підпрограма реалізації метода Шульца ============================================
double** Shulca(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	double EPS = 1e-5;		// точність
	kolOperation = 0;			// кількість операцій метода	
	kolIter = 0;				// кількість ітерацій метода
	double **I2 = oneMatrix(n);	// ініциалізація одиничної матриці однакового розміру з вхідною матрицею
	A_scalar_Num(I2, n, 2);		// добуток одиничної матриці на 2
	/* ----------------------------------------------------------------------------------------
	Обчислення оберненої матриці за методом Ньютона-Шульца
	1. Записати перше наближення: а) транспонувати вхідну матрицю;	б) нормувати її за колонками та рядками
	2. Повторювати ці дії до досягнення заданої точності
	-------------------------------------------------------------------------------------------- 	*/
	// (1) - транспонувати вхідну матрицю і нормувати її за колонками та рядками - перше наближення	 або
	// (2) - транспонувати вхідну матрицю, нормувати її нормою Фробеніуса з добутку вхідної матриці і транспонованої - перше наближення
	double **inv = firstApproxim1(matrix, n); // або =firstApproxim2(matrix, n);		// вибір першого наближення U^(-1)
	while (fabs(detMatrix(A_scalar_B(matrix, inv, n), n) - 1) >= EPS) {	// пока визначник матриці (A * U^(-1) - I) по модулю >= похибки
		double **prev = clone(inv, n);		//A[k-1]	 - запам'ятати попереднє наближення
		inv = A_scalar_B(matrix, prev, n);	//A*(A[k-1]^(-1)) - inv= добуток А і попереднього наближення
		A_scalar_Num(inv, n, -1);			//-A-(A[k-1]^(-1))	- inv помножити на (-1)
		sumMatrix(inv, I2, n);				//2I - A*(A[k-1]^(-1))	- сума (inv + 2*I)
		inv = A_scalar_B(prev, inv, n);		//(A[k-1]^(-1))*(2I - A*(A[k-1]^(-1)))	- добуток попереднього наближення і (inv + 2I) - наше чергове inv і є нашим U^(-1)
		kolIter++;			// збільшити к-сть ітерaцій
		clear(prev, n);		// звільнити пам'ять від попереднього наближення
	}
	return inv;		// повернути обернену матрицю
}
//=============== підпрограма реалізації метода Жардана-Гауса  ================
double** JordanaGausa(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	double **B = clone(matrix, n);	// копія вхідної матриці
	double **I = oneMatrix(n);		// ініциалізація одиничної матриці однакового розміру з вхідною матрицею
	double coeff = 1.0;
	kolIter = 0;			// кількість ітерацій метода
	kolOperation = 0;		// кількість операцій метода		
	for (int step = 0; step < n; step++) {	// цикл по всіх колонках вхідної матриці 
		if (!B[step][step]) {		// якщо ключовий елемент == 0
			bool flag = false;		// прапор обміну рядків
			for (int row = step; row < n && !flag; row++) {		// цикл по всіх рядках, що нижче даного або до першого обміну 
				if (B[row][step]) {	// якщо цей елемент !=0
					swap(B, n, step, row);	// обмін місцями рядків вхідної матриці 
					swap(I, n, step, row);	// обмін місцями рядків оберненої матриці 
					flag = true;			// підняти прапорець
				}
			}
			kolIter += (n - step);	// збільшити к-сть ітерaцій
		}
		coeff = B[step][step];		// отримати ключовий елемент (елемент з головної діагоналі ключової колонки вхідної матриці) 
		if (coeff != 1.0) {
			for (int col = step; col < n; col++)
				B[step][col] /= coeff;		// розділити всі елементи рядка вхідної матриці, в якому розміщений ключовий елемент на ключовий елемент (отримаємо 1 на діагоналі)
			for (int col = 0; col < n; col++)
				I[step][col] /= coeff;		// розділити всі елементи рядка одиничної матриці на той самий ключовий елемент
			kolOperation += (n - step) * 4;	// збільшити к-сть оперaцій
			kolIter += (n - step);	// збільшити к-сть ітерaцій
		}
		for (int row = 0; row < n; row++) {	// цикл по всіх рядках вхідної і одиничної матриць 
			if (row != step) {
				coeff = -B[row][step];		// коефіцієнт - протилежне значення елемента в і-му рядку ключової колонки 
				for (int col1 = step; col1 < n; col1++)
					B[row][col1] += B[step][col1] * coeff;		// сума решти елементів ключового рядка вхідної матриці, помноженого на коефіціент, та елементів і-го рядка (отримаємо 0 у всіх рядках ключової колонки, крім головної діагоналі)
				for (int col1 = 0; col1 < n; col1++)
					I[row][col1] += I[step][col1] * coeff;		// сума всіх елементів ключового рядка одиничної матриці, помноженого на коефіціент, та елементів і-го рядка одиничної матриці
			}
		}
		kolOperation += (n - step) * 2 + n * 2;	// збільшити к-сть оперaцій
		kolIter += n * (n-1);	// збільшити к-сть ітерaцій
	}
	clear(B, n);
	return I;	// повернути адресу на обернену матрицю
}
//=============== підпрограма генерації одиничної в матриці ==============
double** oneMatrix(int n) { // розмірність матриці
	double **matrix = new double*[n];		// створення нової матриці
	for (int row = 0; row < n; row++) {
		matrix[row] = new double[n];		// створення нового рядка матриці
		for (int col = 0; col < n; col++)
			matrix[row][col] = (row == col) ? 1 : 0;	// запис "1" у головну діагональ
	}
	return matrix;		// повернути адресу одиничну матрицю
}
//=============== підпрограма транспонування матриці ================
double** transMatrix(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	double **A0 = clone(matrix, n);			// копія матриці
	for (int row = 0; row < n - 1; row++) 	// саме транспонування 
	for (int col = row + 1; col < n; col++)
		swap(A0[col][row], A0[row][col]);	// поміняти місцями рядки і колонки
	kolIter += (n * n);		// збільшити к-сть ітерaцій
	return A0;		// повернути адресу на транспоновану матрицю 
}
//=============== підпрограма отримання першого наближення (через норми матриці) ================
double** firstApproxim1(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	double N1 = 0, Ninf = 0;			// норми матриці за колонками та рядками
	double **A0 = clone(matrix, n);     // копія матриці
	for (int row = 0; row < n; row++) {
		double colSum = 0, rowSum = 0;	// суми модулів в колонці та рядку
		for (int col = 0; col < n; col++) {
			rowSum += fabs(A0[row][col]);	// обчислити суми модулів в рядку
			colSum += fabs(A0[col][row]);	// обчислити суми модулів в колонці
		}
		Ninf = max(rowSum, Ninf);			// максимальне з сум модулів рядків
		N1 = max(colSum, N1);				// максимальне з сум модулів колонок		
	}
	kolIter += (n * n);				// збільшити к-сть ітерaцій
	kolOperation += (2 * n * n);	// збільшити к-сть оперaцій
	A0 = transMatrix(A0, n);	// транспонувати матрицю
	A_scalar_Num(A0, n, (1 / (N1 * Ninf)));		// нормувати матрицю
	return A0;		// повернути адресу на нормовану матрицю 
}
//=============== підпрограма обчислення норми Фробеніуса матриці ================
double normFrobenius(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	double res = 0.0;		// норма Фробеніуса
	for (int row = 0; row < n; row++)
	for (int col = 0; col < n; col++)
		res += pow(matrix[row][col], 2);		// сума квадратів елементів матриці
	kolIter += (n * n);				// збільшити к-сть ітерaцій
	kolOperation += (2 * n * n);	// збільшити к-сть оперaцій
	return sqrt(res);		// повернути значення норми Фробеніуса
}
//=============== підпрограма отримання першого наближення (через норму Фробеніуса) ================
double** firstApproxim2(double **matrix, int n) { // https://ru.wikipedia.org/wiki/%D0%9E%D0%B1%D1%80%D0%B0%D1%82%D0%BD%D0%B0%D1%8F_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D0%B0
	double **A0 = clone(matrix, n);		// копія матриці
	A0 = transMatrix(A0, n);			// транспонувати матрицю	
	double N_Froben = normFrobenius(A_scalar_B(matrix, A0, n), n);		// норма Фробеніуса добутку матриць (||A*A(t)||)
	A_scalar_Num(A0, n, (1 / N_Froben));		// нормування A(t) за нормою Фробеніуса
	return A0;		// повернути адресу на нормовану матрицю 
}
//=============== підпрограма очистки виділеної пам'яті ================
void clear(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	for (int i = 0; i < n; i++) delete[] matrix[i];	 // очистка пам'яті від рядків
	delete[] matrix;	// очистка пам'яті від змінної матриці
}
//=============== підпрограма створення копії матриці ================
double** clone(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	double **newArr = new double*[n];		// створення нової матриці
	for (int row = 0; row < n; row++) {
		newArr[row] = new double[n];		// створення нового рядка матриці
		for (int col = 0; col < n; col++)
			newArr[row][col] = matrix[row][col];	// скопіювати всі елементи в нову матрицю
	}
	return newArr;		// повернути адресу на копію матриці
}
//=============== підпрограма добутку двох матриць ================
double** A_scalar_B(double **A, double **B, int n) { // адреса матриці А, адреса матриці В, розмірність матриць
	double **result = new double*[n];		// створення нової матриці
	for (int row = 0; row < n; row++) {
		result[row] = new double[n];		// створення нового рядка матриці
		for (int col = 0; col < n; col++) result[row][col] = 0;	// заповнити нулями всі її елементи
	}
	for (int row = 0; row < n; row++) 		// скалярний добуток двох матриць
	for (int col = 0; col < n; col++)
	for (int j = 0; j < n; j++)
		result[row][col] += A[row][j] * B[j][col];	// перемножити дві матриці, результат помістити в result  
	kolOperation += (2 * n * n * n);	// збільшити к-сть оперaцій
	kolIter += (n * n * n);				// збільшити к-сть ітерaцій
	return result;		// повернути адресу на результат множення
}
//=============== підпрограма добутку матриці на число ================
void A_scalar_Num(double **matrix, int n, double a) { // адреса вхідної матриці,  розмірність матриці, множник
	for (int row = 0; row < n; row++)
	for (int col = 0; col < n; col++)
		matrix[row][col] *= a;
	kolOperation += n * n;	// збільшити к-сть оперaцій
	kolIter += (n * n);		// збільшити к-сть ітерaцій
}
//=============== підпрограма обчислення суми двох квадратних матриць ================
void sumMatrix(double **A, double **B, int n) { // адреса матриці А, адреса матриці В, розмірність матриць
	for (int row = 0; row < n; row++)
	for (int col = 0; col < n; col++)
		A[row][col] += B[row][col];		// поелементне додавання, результат буде в першій матриці
	kolOperation += (n * n);	// збільшити к-сть оперaцій
	kolIter += (n * n);			// збільшити к-сть ітерaцій
}
//=============== підпрограма обчислення визначника матриці ================
double detMatrix(double **matrix, int n) {	// адреса вхідної матриці,  розмірність матриці
	double **B = GausStraight(matrix, n);	// привести матрицю до трикутного вигляду методом Гауса (прямим ходом)
	double Det = 1;			// визначник матриці
	for (int i = 0; i < n; i++)
		Det *= B[i][i];		// визначник є добутком елементів головної діагоналі
	kolOperation += n;		// збільшити к-сть оперaцій
	kolIter += n;			// збільшити к-сть ітерaцій
	clear(B, n);			// очистить память від матриці В
	return Det;				// повернути значення визначника
}
//=============== підпрограма приведення матриці до трикутного вигляду методом Гауса (прямий хід) ================
double** GausStraight(double **matrix, int n) { // адреса вхідної матриці,  розмірність матриці
	double **B = clone(matrix, n);		// копія матриці
	for (int step = 0; step < n - 1; step++)
	for (int row = step + 1; row < n; row++) {
		if (!B[step][step]) {	// якщо елемент на головній діагоналі == 0, то шукати серед рядків, що нижче, той, у якого елемент в цій колонці <>0 
			bool flag = false;	// прапорець пошуку
			for (int row1 = step; row1 < n && !flag; row1++) {	// цикл по всіх рядках, що нижче даного
				if (B[row1][step]) {	// якщо знайшли ненульовий елемент, 
					swap(B, n, step, row1);	// то поміняти рядки місцями
					flag = true;		// підняти прапорець
				}
			}
			if (!flag) return B;	// якщо в кожному з рядків, що нижче даного, немає ненульового елементу в цій колонці, повернути матрицю (хоть і не приведену до трикутної - визначник у неї буде ==0)
			kolIter += (n - step);		// збільшити к-сть ітерaцій
		}
		double coeff = -B[row][step] / B[step][step];	// коефіцієнт 
		for (int col = step; col < n; col++)
			B[row][col] += B[step][col] * coeff;		// сума елементів
		kolOperation += (2 * (n - step) + 1);			// збільшити к-сть оперaцій
		kolIter += (n - step);							// збільшити к-сть ітерaцій
	}
	kolIter += (n * n);		// збільшити к-сть ітерaцій
	return B;	// повернути адресу на матрицю
}
//=============== підпрограма обміну місцями двох рядків матриці  ================
void swap(double **matrix, int n, int row_in, int row_to) { // адреса вхідної матриці,  розмірність, № рядка, № рядка
	for (int col = 0; col < n; col++)
		swap(matrix[row_in][col], matrix[row_to][col]);
	kolIter += n;		// збільшити к-сть ітерaцій
}

//=============== підпрограма округлення елементів матриці до вказаного розряду ================
double** RoundMatrix(double **inv, int n, int rozr) { // адреса матриці,  розмірність матриці
	for (int i = 0; i < n; i++) {	// цикл по всіх елементах матриці
		for (int j = 0; j < n; j++)
			inv[i][j] = ((int)(inv[i][j] * rozr + (inv[i][j] >= 0 ? 0.5 : -0.5)) / (float)rozr);
	}
	return inv;
}