[llvm-commits] CVS: llvm/test/Programs/MultiSource/McCat-04-bisect/Makefile README allocvector.c allocvector.h bisect_test.c bisect_test.in dbisect.c dbisect.h
Chris Lattner
lattner at cs.uiuc.edu
Mon May 12 14:00:00 PDT 2003
Changes in directory llvm/test/Programs/MultiSource/McCat-04-bisect:
Makefile added (r1.1)
README added (r1.1)
allocvector.c added (r1.1)
allocvector.h added (r1.1)
bisect_test.c added (r1.1)
bisect_test.in added (r1.1)
dbisect.c added (r1.1)
dbisect.h added (r1.1)
---
Log message:
Initial Checkin
---
Diffs of the changes:
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/Makefile
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/Makefile:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/Makefile Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,7 ----
+ LEVEL = ../../../..
+ PROG = bisect
+ LDFLAGS = -lm
+ #RUN_OPTIONS +=
+ INPUT_FILENAME = bisect_test.in
+ include ../Makefile.multisrc
+
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/README
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/README:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/README Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,34 ----
+ Program for the calculation of the eigenvalues
+ of a symetric tridiagonal matrix by the method of bisection:
+ --------------------------------------------------------------
+ 1) To compile using gcc type:
+
+ gcc -o bisect_test bisect_test.c dbisect.c allocvector.c
+
+ --------------------------------------------------------------
+ 2) To run the benchmark type:
+
+ time bisect_test < bisect_test.in1 > bisect_test.out
+ --------------------------------------------------------------
+ 3) To check whether the results are identical to the results
+ computed by the original program, run
+
+ diff bisect_test.out.orig bisect_test.out
+
+
+ Two possibly useful modifications to the source code can be
+ incorporated by setting the symbols TESTFIRST and RECIPROCAL. Thus if
+ a version of the program where the test (q<0) is moved t the top of
+ the loop, you should compile the program using the command:
+
+ gcc -DTESTFIRST -o bisect_test bisect_test.c dbisect.c allocvector.c
+
+ Similarly, if a version where division is relaced by reiprocal
+ calculation is required, issue the command:
+
+ gcc -DRECIPROCAL -o bisect_test bisect_test.c dbisect.c allocvector.c
+
+ The two changes can also be applied simultaniously:
+
+ gcc -DRECIPROCAL -DTESTFIRST -o bisect_test bisect_test.c dbisect.c
+ allocvector.c
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/allocvector.c
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/allocvector.c:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/allocvector.c Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,60 ----
+
+ /****
+ Copyright (C) 1996 McGill University.
+ Copyright (C) 1996 McCAT System Group.
+ Copyright (C) 1996 ACAPS Benchmark Administrator
+ benadmin at 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 <string.h>
+ #include <unistd.h>
+
+
+ void *allocvector(size_t);
+
+ void *allocvector(size_t size)
+ {
+ void *V;
+
+ /* if (size > sysconf(_SC_PAGESIZE)/2 ) {
+ if ( (V = (void *) valloc(size) ) == NULL ) {
+ fprintf(stderr, "Error: couldn't allocate V in allocvector.\n");
+ exit(-1);
+ }
+ }
+ else { */
+ if ( (V = (void *) malloc((size_t) size)) == NULL ) {
+ fprintf(stderr, "Error: couldn't allocate V in allocvector.\n");
+ exit(-1);
+ }
+ /* } */
+ memset(V,0,size);
+ return V;
+ }
+
+
+ void dallocvector(int n, double **V)
+ {
+ *V = (double *) allocvector((size_t) n*sizeof(double));
+ }
+
+
+ void sallocvector(int n, float **V)
+ {
+ *V = (float *) allocvector((size_t) n*sizeof(float));
+ }
+
+
+ void iallocvector(int n, int **V)
+ {
+ *V = (int *) allocvector((size_t) n*sizeof(int));
+ }
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/allocvector.h
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/allocvector.h:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/allocvector.h Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,9 ----
+ #ifndef __ALLOCVECTOR_H_DEF
+ #define __ALLOCVECTOR_H_DEF
+
+ extern void dallocvector(int n, double **V);
+ extern void sallocvector(int n, float **V);
+ extern void iallocvector(int n, int **V);
+
+
+ #endif
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/bisect_test.c
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/bisect_test.c:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/bisect_test.c Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,83 ----
+
+ /****
+ Copyright (C) 1996 McGill University.
+ Copyright (C) 1996 McCAT System Group.
+ Copyright (C) 1996 ACAPS Benchmark Administrator
+ benadmin at 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"
+ #include "allocvector.h"
+
+
+ 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<n; i++) {
+ B[i] = (double) i;
+ C[i] = (double ) (i+1)*(i+1);
+ C[i] = C[i] * C[i];
+ }
+ }
+
+
+ main(int argc,char *argv[])
+ {
+ int rep,n,k,i,j;
+ double eps,eps2;
+ double *D,*E,*beta,*S;
+
+ scanf("%d",&rep);
+ scanf("%d",&n);
+ scanf("%f",&eps);
+
+ dallocvector(n,&D);
+ dallocvector(n,&E);
+ dallocvector(n,&beta);
+ dallocvector(n,&S);
+
+ for (j=0; j<rep; j++) {
+ test_matrix(n,D,E);
+
+ for (i=0; i<n; i++) {
+ beta[i] = E[i] * E[i];
+ S[i] = 0.0;
+ }
+
+ E[0] = beta[0] = 0;
+ dbisect(D,E,beta,n,1,n,eps,&eps2,&k,&S[-1]);
+
+ }
+
+ for(i=1; i<n; i++)
+ printf("%5d %.15e\n",i+1,S[i]);
+
+ printf("eps2 = %e, k = %d\n",eps2,k);
+
+ return 0;
+ }
+
+
+
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/bisect_test.in
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/bisect_test.in:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/bisect_test.in Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,3 ----
+ 1
+ 500
+ 2.2204460492503131E-16
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/dbisect.c
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/dbisect.c:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/dbisect.c Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,299 ----
+
+ /****
+ Copyright (C) 1996 McGill University.
+ Copyright (C) 1996 McCAT System Group.
+ Copyright (C) 1996 ACAPS Benchmark Administrator
+ benadmin at 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);
+ }
+ }
Index: llvm/test/Programs/MultiSource/McCat-04-bisect/dbisect.h
diff -c /dev/null llvm/test/Programs/MultiSource/McCat-04-bisect/dbisect.h:1.1
*** /dev/null Mon May 12 13:59:16 2003
--- llvm/test/Programs/MultiSource/McCat-04-bisect/dbisect.h Mon May 12 13:59:06 2003
***************
*** 0 ****
--- 1,12 ----
+ #ifndef __DBISECT_H_DEF
+
+ #define __DBISECT_H_DEF
+
+ #define TRACE_ITERATION 0x1
+ #define TRACE_INTERVALS 0x2
+ #define TRACE_INPUT 0x4
+
+ void dbisect(double c[], double b[], double beta[], int n, int m1, int m2,
+ double eps1, double *eps2, int *z, double x[]);
+
+ #endif
More information about the llvm-commits
mailing list