Шпоры по МПиПА / Численные методы / метод вращений Якоби / Исходник / SLAU

#include "stdio.h"
#include "iostream.h"
#include "math.h"

const double PI = 3.1415926536;

// проверка матрицы на диагональность
bool isSimmetrial( double **coefficients, int numberOfEquation ) {
	bool result = true;
	int i, j;

	for ( i = 0; i < numberOfEquation; i++ ) {
		for ( j = i + 1; j < numberOfEquation; j ++ ) {
			if ( coefficients[i][j] != coefficients[j][i] ) {
				result = false;
				break;
			}
		}
		if ( !result ) {
			break;
		}
	}

	return result;
}

// приведение матрицы к диагональному виду "вращением" - на диагонали собственные значения, результат - за сколько шагов с заданнй точностью
int wrachenie( double **coefficients, int numberOfEquation, double **solution, double precision ) {
	int result = 1;
	int i, j, k;
	int maxI, maxJ;
	double max, fi;

	// матрица с помощью которой осуществляем поворот
	double** matricaPoworota;
	matricaPoworota = new double*[numberOfEquation];
	for ( i = 0; i < numberOfEquation; i ++ ) {
		matricaPoworota[i] = new double[numberOfEquation];
	}

	// вспомгательная матрица для вращения
	double** temp;
	temp = new double*[numberOfEquation];
	for ( i = 0; i < numberOfEquation; i ++ ) {
		temp[i] = new double[numberOfEquation];
	}

	// вычисление погрешности - корня из суммы квадратов внедиагональных элементов
	double fault = 0.0;
	for ( i = 0; i < numberOfEquation; i ++ ) {
		for ( j = i + 1; j < numberOfEquation; j ++ ) {
			fault = fault + coefficients[i][j]*coefficients[i][j];
		}
	}
	fault = sqrt( 2*fault );

	// выполняем "вращение" пока погрешность больше заданной точности
	while ( fault > precision ) {

		// ищем максимальный внедиагональный элемент
		max = 0.0;
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = i + 1; j < numberOfEquation; j ++ ) {
				if ( coefficients[i][j] > 0 && coefficients[i][j] > max ) {
					max = coefficients[i][j];
					maxI = i;
					maxJ = j;
				}
				else if ( coefficients[i][j] < 0 && - coefficients[i][j] > max ) {
					max = - coefficients[i][j];
					maxI = i;
					maxJ = j;
				}
			}
		}

		// заполняем матрицу поворота
		for ( i = 0; i < numberOfEquation; i ++ ){
			for ( j = 0; j < numberOfEquation; j ++ ){
				matricaPoworota[i][j] = 0;
			}
			matricaPoworota[i][i] = 1;
		}

		if ( coefficients[maxI][maxI] == coefficients[maxJ][maxJ] ) {
			matricaPoworota[maxI][maxI] = matricaPoworota[maxJ][maxJ] = matricaPoworota[maxJ][maxI] = sqrt( 2.0 ) / 2.0;
			matricaPoworota[maxI][maxJ] = - sqrt( 2.0 ) / 2.0;	
		}
		else {
			fi = 0.5 * atan( ( 2.0 * coefficients[maxI][maxJ] ) / ( coefficients[maxI][maxI] - coefficients[maxJ][maxJ] ) );
			matricaPoworota[maxI][maxI] = matricaPoworota[maxJ][maxJ] = cos( fi );
			matricaPoworota[maxI][maxJ] = - sin( fi );
			matricaPoworota[maxJ][maxI] = sin( fi );
		}
		// закончили заполнять матрицу поворота


		// обнуление вспомогательной матрицы
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				temp[i][j] = 0.0;
			}
		}

		// первая фаза поворота (одно из двух умножений)
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				for ( k = 0; k < numberOfEquation; k ++ ) { 
					temp[i][j] = temp[i][j] + matricaPoworota[k][i] * coefficients[k][j];
				}
			}
		}

		// обнуление поворачиваемой матрицы
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				coefficients[i][j] = 0.0;
			}
		}

		// вторая фаза поворота (пределение новой версии повернутой матрицы)
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				for ( k = 0; k < numberOfEquation; k ++ ) { 
					coefficients[i][j] = coefficients[i][j] + temp[i][k] * matricaPoworota[k][j];
				}
			}
		}
		
		// вычисляем погрешность
		fault = 0.0;
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = i + 1; j < numberOfEquation; j ++ ) {
				fault = fault + coefficients[i][j]*coefficients[i][j];
			}
		}
		fault = sqrt( 2*fault );

		// с помощью матрицы поворота обновляем матрицу содержащую в столбцах собственные векторы исходной матрицы....
		// обнуление вспомогательной матрицы
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				temp[i][j] = 0.0;
			}
		}

		// определяем новые значения матрицы
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				for ( k = 0; k < numberOfEquation; k ++ ) { 
					temp[i][j] = temp[i][j] + solution[i][k] * matricaPoworota[k][j];
				}
			}
		}

		// сохраняем вспомгательную матрицу на нужное место 
		for ( i = 0; i < numberOfEquation; i ++ ) {
			for ( j = 0; j < numberOfEquation; j ++ ) {
				solution[i][j] = temp[i][j];
			}
		}
		result++;
	}

	return result;
}


void main() {
	int i, j;
	int size;
	double **coefficients, **solution, precision;

    cout << "Metod wracheniya.\nWweddite razmernost' matrici: ";
    cin >> size;

	coefficients = new double*[size];
	solution = new double*[size];
	for ( i = 0; i < size; i++ ) {
		coefficients[i] = new double[size];
		solution[i] = new double[size];
	}

	for ( i = 0; i < size; i ++ ){
		for ( j = 0; j < size; j ++ ){
			solution[i][j] = 0;
		}
		solution[i][i] = 1;
	}

	for ( i = 0; i < size; i ++ ){
		cout << "Enter " << i + 1 << " row: ";
		for ( j = 0; j < size; j ++ ){
			cin >> coefficients[i][j];
		}
	}

	cout << "Wwedite tochnost' rascheta: ";
	cin >> precision;
	
	if ( !isSimmetrial( coefficients, size ) ) {
		cout << "Matrica ne simmetrichna";
	}
	else {
		int steps = wrachenie( coefficients, size, solution, precision );
		cout << "Reshenie:\n";
		for ( i = 0; i < size; i++ ) {
			cout << "Sobstwennii wektor nomer " << i + 1 << ":\n";
			for ( j = 0; j < size; j ++ ){
				cout << solution[j][i] << "\n";
			}
		}
		cout << "Sobstwennie znacheniya:\n";
		for ( i = 0; i < size; i++ ) {
			cout << coefficients[i][i] << "\n";
		}

		cout << "Obchee chislo shagow: " << steps;
	}

    cout << "\nPress \"Enter\" to continue..." << endl; 
    getchar();	
}