Метод Гивенса для симметричных матриц

Метод Гивенса основан на преобразовании подобия, аналогич-ном применяемому в методе Якоби. Однако в этом случае алго-ритм построен таким образом, что вновь образованные нулевые элементы при всех последующих преобразованиях сохраняются. Поэтому метод Гивенса требует выполнения конечного числа преобразований и по сравнению с методом Якоби связан с мень-шими затратами машинного времени. Его единственный недоста-ток состоит в том, что симметричная матрица приводится не к диагональному, а к трехдиагональному виду. Ниже будет пока-зано, что такая форма матрицы может быть весьма полезной и оправдывает усилия, затраченные на ее получение.

В случае матрицы размерности п х п метод Гивенса требует п -- 2 основных шагов, на каждом из которых выполняется ряд преобразований, число которых зависит от числа нулей, кото-рое хотят получить в данном столбце или строке. На k -м шаге обращают в нули элементы, стоящие вне трех диагоналей k -й строки и k -го столбца, сохраняя в то же время нулевые элементы, полученные на предыдущих шагах. Таким образом, перед нача-лом k -го шага преобразованная матрица является трехдиа-гональной, если ограничиться рассмотрением ее первых k -- 1 строк и столбцов. По мере преобразований симметричная матри-ца размерности 5х5 приобретает следующие формы:

	*	*	*	*	*
	*	*	*	*	*
A0=	*	*	*	*	*	исходная матрица,
	*	*	*	*	*
	*	*	*	*	*

	*	*
	*	*	*	*	*
A1=		*	*	*	*	после первого основного шага,
		*	*	*	*	состоящего из трех преобразований,
		*	*	*	*

*

*

*

*

*

A2=

*

*

*

*

после второго основного шага,

*

*

*

состоящего из двух преобразований,

*

*

*

	*	*
	*	*	*			после третьего основного шага,
A3=		*	*	*		состоящего из одного преобразования.
			*	*	*	Теперь матрица име-ет трехдиагональный вид.
				*	*

На каждом основном шаге изменяются лишь те элементы мат-рицы а_ij, которые расположены в ее правой нижней (заштрихо-ванной) части. Таким образом на k-м шаге преобразуется только матрица порядка (п -- k + 1), занимающая правый нижний угол исходной матрицы. Ясно, что на каждой следующей стадии вы-полняется меньшее число преобразований, чем на предыдущей. Всего для приведения матрицы к трехдиагональному виду тре-буется выполнить (n ² -- Зп + 2)/2 преобразований.

Наш опыт применения метода Гивенса показывает, что можно при выполнении одного шага преобразований обратить в нуль сразу все элементы целой строки и столбца, стоящие вне трех диагоналей матрицы. Метод, позволяющий выполнить такое преобразование, предложил Хаусхолдер.

Код на Си метода Гивенса:

#include "stdafx.h"

#include <stdlib.h>

#include <math.h>

int _tmain(int argc, _TCHAR* argv[])

{

int N, sum=0;

float c, s, h, g;

bool besc =false,fl=true,fl2=true;

/*Выделяем память под динамические массивы*/

printf("N=");

scanf("%d",&N);

float * b=(float*)malloc(sizeof(float)*N);

float * x=(float*)malloc(sizeof(float)*N);

float ** A;

A=(float**)malloc(sizeof(float*)*N);

for(int i=0; i<N; i++)

{

A[i]=(float*)malloc(sizeof(float)*N);

printf("A[%d]:",i);

for(int k=0; k<N; k++)

{

scanf("%f",&A[i][k]);

}

}

/*----------------------------------------*/

for(int i=0; i<N; i++)

x[i]=0;

printf("b[]:");

for(int i=0; i<N; i++)

{

scanf("%f",&b[i]);

}

for(int i=0; i<N-1; i++)

{

for(int j=i+1; j<N; j++)

{

if(A[i][i]==0 && A[j][i]==0)

int y;

else

{

c=A[i][i]/sqrt(A[i][i]*A[i][i] + A[j][i]*A[j][i]);

s=A[j][i]/sqrt(A[i][i]*A[i][i] + A[j][i]*A[j][i]);

for(int m=i; m<N; m++)

{

h=A[i][m]; g=A[j][m];

A[i][m]=h*c + g*s;

A[j][m]=g*c - h*s;

}

h=b[i]; g=b[j];

b[i]=c*h+s*g;

b[j]=c*g-s*h;

for(int r=0; r<N; r++)

{

for(int t=0; t<N-1; t++)

{

printf("(%f)x%d + ",A[r][t],t);

}

printf("(%f)x%d=%f\n",A[r][N-1],N-1,b[r]);

}

puts("--------------------------------------");

}

}

}

for(int i=N-1; i>=0; i--)

{

sum=0;

for(int j=N-1; j>=i; j--)

sum+=x[j]*A[i][j];

if((A[i][i]==0 && sum==0 && b[i]!=0) || (A[i][i]==0 && sum!=0 && b[i]==0))

{

/*Освобождаем память*/

free(b);

for(int i=0; i<N; i++)

{

free(A[i]);

}

free(A);

free(x);

/*------------------*/

puts("No decisions.");

return 0;

}

else

{

if(A[i][i]==0 || (A[i][i]==0 && sum==0 && b[i]==0))

{

x[1]=0;

besc=true;

}

else

{

x[i]=(b[i]-sum)/A[i][i];

}

}

}

for(int r=0; r<N; r++)

{

printf("x%d=%f\n",r,x[r]);

}

for(int r=0; r<N; r++)

{

if(b[r]!=0)

fl=false;

}

for(int r=0; r<N; r++)

{

for(int t=0; t<N; t++)

{

if(A[r][t]!=0)

fl2=false;

}

}

if(fl)

besc=false;

if(fl && fl2)

besc=true;

if(besc)

puts("...and infinite set of decisions");

/*Освобождаем память*/

free(b);

for(int i=0; i<N; i++)

{

free(A[i]);

}

free(A);

free(x);

/*------------------*/

return 0;

}

Исходные матрицы:

А= В=

Введите размерность квадратной матрицы N: 5

Введите коэффициенты 1 строки матрицы A:

9 -9 -5 -4 -7

Введите коэффициенты 2 строки матрицы A:

9 -5 1 9 -6

Введите коэффициенты 3 строки матрицы A:

-5 7 5 5 6

Введите коэффициенты 4 строки матрицы A:

2 8 -1 3 -2

Введите коэффициенты 5 строки матрицы A:

-4 -9 -8 -6 5

Введите матрицу вектор B:

5 -4 -6 -9 -7

Решение системы:

X[1] = -0,53601

X[2] = -0,57755

X[3] = 1,96475

X[4] = -1,32067

X[5] = -1,30960

Проверка 1 уравнения:

Левая часть уравнения = 5

Правая часть уравнения = 5

Проверка 2 уравнения:

Левая часть уравнения = -4

Правая часть уравнения = -4

Проверка 3 уравнения:

Левая часть уравнения = -6

Правая часть уравнения = -6

Проверка 4 уравнения:

Левая часть уравнения = -9

Правая часть уравнения = -9

Проверка 5 уравнения:

Левая часть уравнения = -7

Правая часть уравнения = -7

Невязки:

X[1] = 8.88178e-016

X[2] = 2.66454e-015

X[3] = 8.88178e-016

X[4] = 0

X[5] = 1.77636e-015

Для продолжения нажмите любую клавишу...