Тической физики Метод вращений (Якоби) для нахождения собственных значений и собственных векторов матриц

4. Программа на C++ для решения СЛАУ методом Якоби

#include

#include "nrutil.h" /* Здесь определяются некоторые утилиты типа выделения памяти */
/* Преобразование элементов при ротации */
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);a[k][l]=h+s*(g-h*tau)
/* максимальное число проходов */
#define MAXSWEEP 50
/* Программа jacobi вычисляет все собственные значения и собственные векторы
действительной симметричной матрицы a[1...n][1...n]. На выходе разрушаются
все наддиагональные элементы a. d[1...n] возвращает собственные значения
матрицы, v[1...n][1...n] -- в столбцах нормализованные собственные векторы,
nrot -- число ротаций Якоби, потребовавшихся для данной программы.
*/
void jacobi(float **a, int n, float *d, float **v, int *nrot) {
int j, iq, ip, i;
float tresh, theta, tau, t, sm, s, h, g, c, *b, *z;
/* выделить память для временно используемых векторов (декларации в nrutil.h) */
b=vector(1,n); z=vector(1,n);
/* инициализировать v как единичную матрицу */
for(ip=1;ip<=n;ip++) {
for(iq=1;iq<=n;iq++) v[ip][iq]=0.;
v[ip][ip]=1.;
}
/* инициализировать b диагональю a, z -- нулем */
for(ip=1;ip<=n;ip++) {b[ip]=a[ip][ip]; z[ip]=0.;}
/* вначале число ротаций нулевое */
*nrot=0;
/* делаем не более MAXSWEEP проходов */
for(i=1;i<=MAXSWEEP;i++) {
/* вычисляем сумму модулей внедиагональных элементов */
for(sm=0.,ip=1;ip<=n;ip++) for(iq=ip+1;iq<=n;iq++) sm += fabs(a[ip][iq]);
/* диагональная матрица -> нормальный выход */
if(sm==0.) {
free_vector(z,1,n); free_vector(b,1,n); return;
}
/* пороговое значение элемента для ротации */
tresh=(i<4 ? 0.2*sm/(n*n) : 0.);
/* проход осуществляется по строкам, в каждой строке по столбцам */
for(ip=1;ip<=n-1;ip++) for(iq=ip+1;iq<=n;iq++) {
/* отследить случай малого элемента после 4 проходов */
g=100.*fabs(a[ip][iq]);
if(i>4 && (float)fabs(d[ip]+g)==(float)fabs(d[ip])
&& (float)fabs(d[iq]+g)==(float)fabs(d[iq])) a[ip][iq]=0.;
/* и случай малого элемента на первых 3 проходах
(обработать только превысившие порог) */
else if(fabs(a[ip][iq])>tresh) {
h=d[ip]-d[iq];
/* вычислить значение t=s/c по формуле корня квадратного уравнения */
if((float)(fabs(h)+g)==(float)fabs(h)) t=a[ip][iq]/h;
else {
theta=0.5*h/a[ip][iq];
t=1./(fabs(theta)+sqrt(1.+theta*theta));
if(theta<0.) t = -t;
}
/* вычислить c, s, tau, и др. Изменить диагональ. Обнулить (ip,iq)-элемент */
c=1./sqrt(1+t*t); s=t*c; tau=s/(1.+c); h=t*a[ip][iq];
z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h;
a[ip][iq]=0.;
/* поворот при 1<=jfor(j=1;j<=ip-1;j++) {ROTATE(a,j,ip,j,iq);}
/* поворот при ipfor(j=ip+1;j<=iq-1;j++) {ROTATE(a,ip,j,iq,j);}
/* поворот при iqfor(j=iq+1;j<=n;j++) {ROTATE(a,ip,j,j,iq);}
/* добавка для матрицы собственных векторов */
for(j=1;j<=n;j++) {ROTATE(v,j,ip,j,iq);}
/* приращение счетчика ротаций */
++(*nrot);
}
}
/* добавить до диагонали и реинициализировать z */
for(ip=1;ip<=n;ip++) {
b[ip] += z[ip]; d[ip]=b[ip]; z[ip]=0.;
}
}
/* если мы здесь, то число ротаций превысило лимит. Функция nerror (выход с диагностикой ошибки) описана декларацией в nrutil.h. */
nerror("Too many iterations in the routine jacobi");
}