|
ru.algorithms
- RU.ALGORITHMS ----------------------------------------------------------------
From : Evgenij Masherov 2:5020/175.2 07 Feb 2002 11:10:23
To : Edward Avernin
Subject : Re: задача нахождения соб.векторов и соб.чисел
--------------------------------------------------------------------------------
Wed Feb 06 2002 22:31, Edward Avernin wrote to Evgenij Masherov:
EA> From: "Edward Avernin" <edward@atnet.ru>
EA> Здравствуйте Евгений!
EA> У меня задача - рассчитать любую матрицу(она также может быть и
EA> симметричной, но не всегда).
Вот программка, реализующая метод Якоби
(Матрица симметрична, хранится общим способом, в процессе счета слегка
портится...)
Оригинальность алгоритма - в реализации поиска максимального элемента за O(N)
шагов и использовании вместо обычных вращений - по Джентльмену.
void jacm(int n,double **a,double *d,double **v)
{
int i,j,p,q,*ind,*flag,fall;
double c,s,t,h,g,theta,maxa,t1,t2,tmp,*amax,*w,prec;
if((amax=(double*)calloc(n,sizeof(double)))==NULL)goto E0;
if((w=(double*)calloc(n,sizeof(double)))==NULL)goto E1;
if((ind=(int*)calloc(n,sizeof(int)))==NULL)goto E2;
if((flag=(int*)calloc(n,sizeof(int)))==NULL)goto E3;
prec=1e-12;
for(p=0;p<n;p++)
{
d[p]=a[p][p];
w[p]=1e0;
for(q=0;q<n;q++) v[p][q]=0;
v[p][p]=1e0;
}
for(p=0;p<n-1;p++)
{ind[p]=0;amax[p]=0;flag[p]=1;
for(q=p+1;q<n;q++) {if (fabs(a[p][q])>=amax[p])
{amax[p]=fabs(a[p][q]);ind[p]=q;}}
}
amax[n-1]=-1;ind[n-1]=0;fall=1;
while(fall!=0)
{maxa=0;p=-1;q=0;
for(i=0;i<n;i++)
{if ((amax[i]>maxa)&&(flag[i]!=0)){p=i;maxa=amax[i];q=ind[p];}}
if(p<0) fall=0;
else if (amax[p]<((fabs(d[p]+1e0))*prec)) flag[p]=0;
else {
h=d[q]-d[p];
if (fabs(a[p][q])<(fabs(h)*prec)) t=a[p][q]/h;
else
{theta=0.5*h/a[p][q];t=1e0/(fabs(theta)+sqrt((double)(1e0+theta*theta)));
if (theta<0) t=-t;}
c=1e0/sqrt((double)(1e0+t*t));
s=t*c;
h=t*a[p][q];
d[p]-=h;d[q]+=h;
a[p][q]=0;
amax[p]=0;amax[q]=0;
for(j=0;j<p;j++)
{g=a[j][p];h=a[j][q];
a[j][p]=c*g-s*h;
a[j][q]=c*h+s*g;}
for(j=p+1;j<q;j++)
{g=a[p][j];h=a[j][q];
a[p][j]=c*g-s*h;
a[j][q]=c*h+s*g;}
for(j=q+1;j<n;j++)
{g=a[p][j];h=a[q][j];
a[p][j]=c*g-s*h;
a[q][j]=c*h+s*g;}
if (fabs(c)>fabs(s))
{t1=t*w[q]/w[p];t2=t*w[p]/w[q];
w[p]=w[p]*c;w[q]=w[q]*c;
for(j=0;j<n;j++)
{g=v[j][p];h=v[j][q];
v[j][p]=g-t1*h;
v[j][q]=t2*g+h;}}
else {t1=w[p]/w[q]/t;t2=w[q]/w[p]/t;
tmp=w[q]*s;w[q]=w[p]*s;w[p]=tmp;
for(j=0;j<n;j++)
{g=v[j][p];h=v[j][q];
v[j][p]=t1*g-h;
v[j][q]=g+t2*h;}}
for(j=0;j<n-1;j++)
{
flag[j]=1;
if((j==p)||(j==q)||(ind[j]==p)||(ind[j]==q))
{
amax[j]=0;ind[j]=0;
for(i=j+1;i<n;i++)
{
tmp=fabs(a[j][i]);
if(amax[j]<tmp)
{
amax[j]=tmp;
ind[j]=i;
}
}
}
}
}
}
for(p=0;p<n;p++)
for(q=0;q<n;q++)
{v[p][q]*=w[q];
if(p>q) a[q][p]=a[p][q];}
for(p=0;p<n-1;p++)
for(q=p+1;q<n;q++)
if(d[p]<d[q])
{tmp=d[p];d[p]=d[q];d[q]=tmp;
for(j=0;j<n;j++)
{tmp=v[j][p];v[j][p]=v[j][q];v[j][q]=tmp;}
}
free(flag);
E3:
free(ind);
E2:
free(w);
E1:
free(amax);
E0:;
return;
}
EA> Мне нужно вывести все СЗ, и обязательно все СВ, интервал не задан. И
EA> кажется
EA> у меня даже описание
EA> итерационного метода не полное(многое не понятно). Hет ли у Вас
EA> алгоритмов и примеров с решением этой задачи,
EA> сейчас подойдет даже итерационный метод(кстати он случайно не для
EA> нахождения max СЗ ?).
Метод прямых итераций сходится к максимальному СЗ. Hебольшой модификацией
исходной матрицы он приспосабливается к нахождению минимального СЗ.
Метод обратных итераций дает СВ для заданного СЗ (которое находится QR или
бисекциями).
Евгений Машеров АКА СанитарЖеня
--- ifmail v.2.15
* Origin: FidoNet Online - http://www.fido-online.com (2:5020/175.2)
Вернуться к списку тем, сортированных по: возрастание даты уменьшение даты тема автор
|
