From f9ebf19ba3ca4c3ee67cc88bbea407d4dd734249 Mon Sep 17 00:00:00 2001 From: xleroy Date: Wed, 17 Feb 2010 13:25:03 +0000 Subject: 3 more benchmarks git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@1252 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e --- test/c/bisect.c | 376 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 376 insertions(+) create mode 100644 test/c/bisect.c (limited to 'test/c/bisect.c') diff --git a/test/c/bisect.c b/test/c/bisect.c new file mode 100644 index 00000000..759f9978 --- /dev/null +++ b/test/c/bisect.c @@ -0,0 +1,376 @@ +/**** + 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 +#include +#include +#include +#include + +void *allocvector(size_t size) +{ + void *V; + + if ( (V = (void *) malloc((size_t) size)) == NULL ) { + fprintf(stderr, "Error: couldn't allocate V in allocvector.\n"); + exit(2); + } + memset(V,0,size); + return V; +} + +void dallocvector(int n, double **V) +{ + *V = (double *) allocvector((size_t) n*sizeof(double)); +} + + +#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; + + a = 0; + q = 1.0; + + for(i=0; i 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); + } +} + +void test_matrix(int n, double *C, double *B) +/* Symmetric tridiagonal matrix with diagonal + + c_i = i^4, i = (1,2,...,n) + + and off-diagonal elements + + b_i = i-1, i = (2,3,...n). + It is possible to determine small eigenvalues of this matrix, with the + same relative error as for the large ones. +*/ +{ + int i; + + for(i=0; i