| |
| /**** |
| 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. |
| ****/ |
| |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <math.h> |
| #include <float.h> |
| #include "dbisect.h" |
| |
| #define FUDGE (double) 1.01 |
| |
| |
| |
| int sturm(int n, double c[], double b[], double beta[], double x) |
| |
| /************************************************************************** |
| |
| Purpose: |
| ------------ |
| Calculates the sturm sequence given by |
| |
| q_1(x) = c[1] - x |
| |
| q_i(x) = (c[i] - x) - b[i]*b[i] / q_{i-1}(x) |
| |
| and returns a(x) = the number of negative q_i. a(x) gives the number |
| of eigenvalues smaller than x of the symmetric tridiagonal matrix |
| with diagonal c[0],c[1],...,c[n-1] and off-diagonal elements |
| b[1],b[2],...,b[n-1]. |
| |
| |
| Input parameters: |
| ------------------------ |
| n : |
| The order of the matrix. |
| |
| c[0]..c[n-1] : |
| An n x 1 array giving the diagonal elements of the tridiagonal matrix. |
| |
| b[1]..b[n-1] : |
| An n x 1 array giving the sub-diagonal elements. b[0] may be |
| arbitrary but is replaced by zero in the procedure. |
| |
| beta[1]..beta[n-1] : |
| An n x 1 array giving the square of the sub-diagonal elements, |
| i.e. beta[i] = b[i]*b[i]. beta[0] may be arbitrary but is replaced by |
| zero in the procedure. |
| |
| x : |
| Argument for the Sturm sequence. |
| |
| |
| Returns: |
| ------------------------ |
| integer a = Number of eigenvalues of the matrix smaller than x. |
| |
| |
| Notes: |
| ------------------------ |
| On SGI PowerChallenge this function should be compiled with option |
| "-O3 -OPT:IEEE_arithmetic=3" in order to activate the optimization |
| mentioned in the code below. |
| |
| |
| **********************************************************************/ |
| |
| { |
| int i; |
| int a; |
| double q,iq; |
| |
| a = 0; |
| q = 1.0; |
| |
| for(i=0; i<n; i++) { |
| |
| #ifndef TESTFIRST |
| |
| if (q != 0.0) { |
| |
| #ifndef RECIPROCAL |
| q = (c[i] - x) - beta[i]/q; |
| #else |
| /* A potentially NUMERICALLY DANGEROUS optimizations is used here. |
| * The previous statement should read: |
| * q = (c[i] - x) - beta[i]/q |
| * But computing the reciprocal might help on some architectures |
| * that have multiply-add and/or reciprocal instuctions. |
| */ |
| iq = 1.0/q; |
| q = (c[i] - x) - beta[i]*iq; |
| #endif |
| |
| } |
| else { |
| q = (c[i] - x) - fabs(b[i])/DBL_EPSILON; |
| } |
| |
| if (q < 0) |
| a = a + 1; |
| } |
| |
| #else |
| |
| if (q < 0) { |
| a = a + 1; |
| |
| #ifndef RECIPROCAL |
| q = (c[i] - x) - beta[i]/q; |
| #else |
| /* A potentially NUMERICALLY DANGEROUS optimizations is used here. |
| * The previous statement should read: |
| * q = (c[i] - x) - beta[i]/q |
| * But computing the reciprocal might help on some architectures |
| * that have multiply-add and/or reciprocal instuctions. |
| */ |
| iq = 1.0/q; |
| q = (c[i] - x) - beta[i]*iq; |
| #endif |
| |
| } |
| else if (q > 0.0) { |
| #ifndef RECIPROCAL |
| q = (c[i] - x) - beta[i]/q; |
| #else |
| /* A potentially NUMERICALLY DANGEROUS optimizations is used here. |
| * The previous statement should read: |
| * q = (c[i] - x) - beta[i]/q |
| * But computing the reciprocal might help on some architectures |
| * that have multiply-add and/or reciprocal instuctions. |
| */ |
| iq = 1.0/q; |
| q = (c[i] - x) - beta[i]*iq; |
| #endif |
| } |
| else { |
| q = (c[i] - x) - fabs(b[i])/DBL_EPSILON; |
| } |
| |
| } |
| if (q < 0) |
| a = a + 1; |
| #endif |
| |
| return a; |
| } |
| |
| |
| |
| |
| void dbisect(double c[], double b[], double beta[], |
| int n, int m1, int m2, double eps1, double *eps2, int *z, |
| double x[]) |
| |
| |
| /************************************************************************** |
| |
| Purpose: |
| ------------ |
| |
| Calculates eigenvalues lambda_{m1}, lambda_{m1+1},...,lambda_{m2} of |
| a symmetric tridiagonal matrix with diagonal c[0],c[1],...,c[n-1] |
| and off-diagonal elements b[1],b[2],...,b[n-1] by the method of |
| bisection, using Sturm sequences. |
| |
| |
| Input parameters: |
| ------------------------ |
| |
| c[0]..c[n-1] : |
| An n x 1 array giving the diagonal elements of the tridiagonal matrix. |
| |
| b[1]..b[n-1] : |
| An n x 1 array giving the sub-diagonal elements. b[0] may be |
| arbitrary but is replaced by zero in the procedure. |
| |
| beta[1]..beta[n-1] : |
| An n x 1 array giving the square of the sub-diagonal elements, |
| i.e. beta[i] = b[i]*b[i]. beta[0] may be arbitrary but is replaced by |
| zero in the procedure. |
| |
| n : |
| The order of the matrix. |
| |
| m1, m2 : |
| The eigenvalues lambda_{m1}, lambda_{m1+1},...,lambda_{m2} are |
| calculated (NB! lambda_1 is the smallest eigenvalue). |
| m1 <= m2must hold otherwise no eigenvalues are computed. |
| returned in x[m1-1],x[m1],...,x[m2-1] |
| |
| eps1 : |
| a quantity that affects the precision to which eigenvalues are |
| computed. The bisection is continued as long as |
| x_high - x_low > 2*DBL_EPSILON(|x_low| + |x_high|) + eps1 (1) |
| When (1) is no longer satisfied, (x_high + x_low)/2 gives the |
| current eigenvalue lambda_k. Here DBL_EPSILON (constant) is |
| the machine accuracy, i.e. the smallest number such that |
| (1 + DBL_EPSILON) > 1. |
| |
| Output parameters: |
| ------------------------ |
| |
| eps2 : |
| The overall bound for the error in any eigenvalue. |
| z : |
| The total number of iterations to find all eigenvalues. |
| x : |
| The array x[m1],x[m1+1],...,x[m2] contains the computed eigenvalues. |
| |
| **********************************************************************/ |
| { |
| int i; |
| double h,xmin,xmax; |
| beta[0] = b[0] = 0.0; |
| |
| |
| /* calculate Gerschgorin interval */ |
| xmin = c[n-1] - FUDGE*fabs(b[n-1]); |
| xmax = c[n-1] + FUDGE*fabs(b[n-1]); |
| for(i=n-2; i>=0; i--) { |
| h = FUDGE*(fabs(b[i]) + fabs(b[i+1])); |
| if (c[i] + h > xmax) xmax = c[i] + h; |
| if (c[i] - h < xmin) xmin = c[i] - h; |
| } |
| |
| /* printf("xmax = %lf xmin = %lf\n",xmax,xmin); */ |
| |
| /* estimate precision of calculated eigenvalues */ |
| *eps2 = DBL_EPSILON * (xmin + xmax > 0 ? xmax : -xmin); |
| if (eps1 <= 0) |
| eps1 = *eps2; |
| *eps2 = 0.5 * eps1 + 7 * *eps2; |
| { int a,k; |
| double x1,xu,x0; |
| double *wu; |
| |
| if( (wu = (double *) calloc(n+1,sizeof(double))) == NULL) { |
| fputs("bisect: Couldn't allocate memory for wu",stderr); |
| exit(1); |
| } |
| |
| /* Start bisection process */ |
| x0 = xmax; |
| for(i=m2; i >= m1; i--) { |
| x[i] = xmax; |
| wu[i] = xmin; |
| } |
| *z = 0; |
| /* loop for the k-th eigenvalue */ |
| for(k=m2; k>=m1; k--) { |
| xu = xmin; |
| for(i=k; i>=m1; i--) { |
| if(xu < wu[i]) { |
| xu = wu[i]; |
| break; |
| } |
| } |
| if (x0 > x[k]) |
| x0 = x[k]; |
| |
| x1 = (xu + x0)/2; |
| while ( x0-xu > 2*DBL_EPSILON*(fabs(xu)+fabs(x0))+eps1 ) { |
| *z = *z + 1; |
| |
| /* Sturms Sequence */ |
| |
| a = sturm(n,c,b,beta,x1); |
| |
| /* Bisection step */ |
| if (a < k) { |
| if (a < m1) |
| xu = wu[m1] = x1; |
| else { |
| xu = wu[a+1] = x1; |
| if (x[a] > x1) x[a] = x1; |
| } |
| } |
| else { |
| x0 = x1; |
| } |
| x1 = (xu + x0)/2; |
| } |
| x[k] = (xu + x0)/2; |
| } |
| free(wu); |
| } |
| } |
