| |
| /**** |
| Copyright (C) 1996 McGill University. |
| Copyright (C) 1996 McCAT System Group. |
| Copyright (C) 1996 ACAPS Benchmark Administrator |
| benadmin@acaps.cs.mcgill.ca |
| |
| This program is free software; you can redistribute it and/or modify |
| it provided this copyright notice is maintained. |
| |
| This program is distributed in the hope that it will be useful, |
| but WITHOUT ANY WARRANTY; without even the implied warranty of |
| MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. |
| ****/ |
| |
| /************************************************************************/ |
| /* author : Mikkel Damsgaard */ |
| /* Kirsebaerhaven 85b */ |
| /* */ |
| /* DK-8520 Lystrup */ |
| /* email mikdam@daimi.aau.dk */ |
| /* */ |
| /* files : */ |
| /* Divsol.c QRfact.h Divsol.h Jacobi.c */ |
| /* Jacobi.h Triang.c print.c MM.c */ |
| /* Triang.h print.h MM.h QRfact.c */ |
| /* main.c main.h */ |
| /* */ |
| /* It calculates the eigenvalues for 4 different matrixes. It does not */ |
| /* take any input; those 4 matrixes are calculated by MakeMatrix */ |
| /* function. Output is given as 4 files: val2, val3, val4, val5, that */ |
| /* contains the eigenvalues for each of the matrixes. */ |
| /* */ |
| /************************************************************************/ |
| #include "QRfact.h" |
| #include "MM.h" |
| |
| /* Calculate givens transformation */ |
| void Givens(double a,double b,double *s,double *c) |
| { |
| double t; |
| if (b==0.0) {*c=1;*s=0;} |
| else if (fabs(b)>fabs(a)) |
| { |
| t = -a/b; |
| *s = 1/sqrt(1+t*t); |
| *c = (*s)*t; |
| } |
| else |
| { |
| t = -b/a; |
| *c = 1/sqrt(1+t*t); |
| *s = (*c)*t; |
| } |
| /* printf("%f %f %f %f\n",a,b,(*c)*a-(*s)*b,(*c)*b+(*s)*a); */ |
| } |
| |
| int sign(double a) |
| { |
| if (a<0) return -1; |
| return 1; |
| } |
| |
| |
| void ApplyRGivens(Matrix U,double s, double c,int i,int j) |
| { /* U = U*G, G=Givens(i,j,theta) */ |
| int k; |
| double t1,t2; |
| |
| for (k=0;k<n;k++) |
| { |
| t1=U[k][i]; t2=U[k][j]; |
| U[k][i] = c*t1-s*t2; U[k][j] = s*t1+c*t2; |
| } |
| } |
| |
| |
| Matrix QRiterate(Matrix A, Matrix U) |
| { |
| double c,s,t,a,b,app,aqq,apq,a1p,a1q,a4p,a4q; |
| double d,mu,x,z,t1,t2; |
| |
| int i,j,k,notdone=1,p,q,l,m; |
| |
| while (notdone) |
| { |
| /* First examine if any offdiagonal is small enough to elimante */ |
| for (i=0;i<n-1;i++) |
| if (fabs(A[i+1][i])<(fabs(A[i][i])+fabs(A[i+1][i+1]))*epsilon) |
| A[i+1][i]=A[i][i+1]=0.0; |
| |
| /* Find the boundaries for the QR step */ |
| q=n-1; |
| while ((q>0) && (A[q-1][q]==0.0)) q--; |
| if (q==0) |
| { |
| notdone=0; |
| } |
| else |
| { |
| p=q; |
| while ((p>0) && (A[p-1][p]!=0.0)) p--; |
| } |
| |
| if (!notdone) break; |
| /* printf("%e %i\n",A[q-1][q],q); */ |
| /* Now we do the QR step on the submatrice A[p..q][p..q] */ |
| |
| /* First calculate the shift */ |
| d = (A[q-1][q-1]-A[q][q])/2; t=A[q][q-1]; t=t*t; |
| mu = A[q][q]-(t/(d+sign(d)*sqrt(d*d+t))); |
| x = A[p][p]-mu; |
| z = A[p+1][p]; |
| |
| /* Now QR faktorise */ |
| for (i=p;i<q;i++) |
| { |
| /* Find the givens rotation */ |
| Givens(x,z,&s,&c); |
| |
| /* l=max(p,i-1); */ |
| |
| l = (i-1>p)?i-1:p; |
| |
| /* m=min(q,i+2); */ |
| |
| m = (q<i+2)?q:i+2; |
| |
| for (k=l;k<=m;k++) |
| { |
| t1=A[i][k]; t2=A[i+1][k]; |
| A[i][k] = c*t1-s*t2; A[i+1][k] = s*t1+c*t2; |
| } |
| |
| for (k=l;k<=m;k++) |
| { |
| t1=A[k][i]; t2=A[k][i+1]; |
| A[k][i] = c*t1-s*t2; A[k][i+1] = s*t1+c*t2; |
| } |
| |
| |
| /* Store the givens in U */ |
| /* U = G.T()*U */ |
| ApplyRGivens(U,s,c,i,i+1); |
| |
| /* And find the next pair of troublemakers */ |
| if (i<q-1) |
| { |
| x=A[i+1][i]; |
| z=A[i+2][i]; |
| } |
| } |
| } |
| |
| } |