| |
| /**** |
| 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 "Triang.h" |
| #include "MM.h" |
| #include "print.h" |
| #include "QRfact.h" |
| |
| |
| double norm(Matrix A,int col,int rs,int re) |
| { |
| double t=0.0; |
| int i; |
| |
| for (i=rs;i<=re;i++) |
| t += A[i][col]*A[i][col]; |
| |
| return sqrt(t); |
| } |
| |
| void House(Matrix A,Vector v,int col,int sr,int er) |
| { |
| double a,b; |
| int m; |
| |
| a=norm(A,col,sr,er); |
| b=1/(A[sr][col]+sign(A[sr][col])*a); |
| |
| for (m=sr+1;m<=er;m++) |
| v[m]=A[m][col]*b; |
| |
| v[sr] = 1.0; |
| } |
| |
| double xty(Vector x,Vector y,int s,int e) |
| { |
| double t=0.0; |
| int i; |
| |
| for (i=s;i<=e;i++) |
| t += x[i]*y[i]; |
| |
| return t; |
| } |
| |
| Matrix Trianglelise(Matrix A,int i) |
| { |
| Matrix U,P,T; |
| |
| double a,b; |
| Vector v,w,p; |
| |
| int j,k,l,m,h; |
| |
| /* Initialize U to be the I */ |
| P = newMatrix(); |
| U = newIdMatrix(); |
| |
| v=(Vector)malloc(sizeof(double)*n); |
| w=(Vector)malloc(sizeof(double)*n); |
| p=(Vector)malloc(sizeof(double)*n); |
| |
| if (i<2) return A; |
| l=i; /* This is the number of non-zero off diagonal entries */ |
| for (j=0;j<n-2;j++) /* Iterate through the columns */ |
| { |
| |
| /* This is to save time, h = min(j+l+2,n-1) */ |
| h = (j+l+i<n-1)?j+l+i:n-1; |
| |
| /* Find the householder vektor */ |
| House(A,v,j,j+1,h); |
| |
| /* Now v[j+1:j+l] contains the householder vektor */ |
| |
| /* Well, we need to apply the thing pre and |
| post multiplikative */ |
| |
| /* b = 1/(v.t()*v) */ |
| b=1.0/xty(v,v,j+1,h); |
| |
| for (m=j;m<=h;m++) |
| p[m] = 2.0*xty(A[m],v,j+1,h)*b; |
| |
| a=xty(p,v,j+1,h)*b; |
| /* w = p - a*v */ |
| for (m=j+1;m<=h;m++) |
| w[m] = p[m] - a*v[m]; |
| |
| /* Update A[j+1:h][j+1:h], and utilize that we */ |
| /* know A to be symmetric before and after the update */ |
| for (m=j+1;m<=h;m++) |
| for (k=m;k<=h;k++) |
| { |
| A[m][k] -= v[m]*w[k] + w[m]*v[k]; |
| A[k][m] = A[m][k]; |
| } |
| |
| /* Calculating the off diagonal entry */ |
| |
| A[j+1][j] = A[j][j+1] = A[j][j+1] - p[j]; |
| |
| /* And zero the rest of the column/row */ |
| for (m=j+2;m<=h;m++) |
| A[j][m] = A[m][j] = 0.0; |
| |
| /* Calculate U*P */ |
| /* Only update colums j+1 through j+l */ |
| /* Remember that U is NOT symmetric */ |
| /* TIME l*n+n*2*l = 3*n*l */ |
| |
| for (m=0;m<n;m++) |
| w[m]=2.0*b*xty(U[m],v,j+1,h); |
| |
| for (m=0;m<n;m++) |
| for (k=j+1;k<=h;k++) |
| U[m][k] -= w[m]*v[k]; |
| |
| /* Update l */ |
| if (j+l+(i-1)<n-1) l += (i-1); else l = (n-1)-(j+1); |
| /* Check(A,U,i); */ |
| } |
| |
| free(v); |
| free(w); |
| free(p); |
| return U; |
| |
| } |
| |
| |