| |
| /**** |
| 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 "print.h" |
| |
| void ApplyGivens(Matrix A,double s,double c, int i, int j,int start,int end) |
| { /* A = G.t()AG, G = Givens(i,j,theta); */ |
| double t1,t2; |
| int k; |
| |
| for (k=start;k<=end;k++) |
| { /* A = G.t()*A */ |
| t1=A[i][k]; t2=A[j][k]; |
| A[i][k] = c*t1-s*t2; A[j][k] = s*t1+c*t2; |
| } |
| |
| for (k=start;k<=end;k++) |
| { /* A = A*G */ |
| t1=A[k][i]; t2=A[k][j]; |
| A[k][i] = c*t1-s*t2; A[k][j] = s*t1+c*t2; |
| } |
| } |
| |
| Matrix Jacobi(Matrix A,int bw) |
| { |
| double a,b,s,c; |
| int i,j,k,l,m; |
| Matrix U; |
| |
| U = newIdMatrix(); |
| |
| for (i=bw;i>=2;i--) |
| { /* Remove one bandwith at a time */ |
| for (j=0;j<n-i;j++) |
| { /* Iterate through the diagonal */ |
| |
| /* Calculate the Givens to zero A[j][j+i] */ |
| Givens(A[j][j+i-1],A[j][j+i],&s,&c); |
| |
| /* Now apply that givens */ |
| |
| /* A = G.t()*A*G */ |
| ApplyGivens(A,s,c,j+i-1,j+i,j,(j+2*i<n)?j+2*i:n-1); |
| /* U = U*G */ |
| ApplyRGivens(U,s,c,j+i-1,j+i); |
| |
| /* Now the bandwith is i+1, and we want it to return to i */ |
| /* so we chase the unwanted non-zero element down the diagonal */ |
| /* The unwanted non-zero element is A[j+i-1][j+2i+1] */ |
| |
| /* Check(A,U,bw); */ |
| |
| m = j+i; |
| while (m<n-i) |
| { |
| /* We determine the givens that will zero the unwanted |
| non-zero element */ |
| |
| Givens(A[m-1][m+i-1],A[m-1][m+i],&s,&c); |
| |
| /* And apply that givens */ |
| |
| /* A = G.t()*A*G */ |
| ApplyGivens(A,s,c,m+i-1,m+i,m-1,(m+2*i<n)?m+2*i:n-1); |
| /* U = U*G */ |
| ApplyRGivens(U,s,c,m+i-1,m+i); |
| |
| /* Check(A,U,bw); */ |
| |
| m += i; |
| } |
| |
| /* Now the ith band is zero in the first j positions */ |
| } |
| } |
| return U; |
| } |