modernc.org/ccgo/v3@v3.16.14/lib/testdata/CompCert-3.6/test/c/bisect.c (about)

     1  /****
     2      Copyright (C) 1996 McGill University.
     3      Copyright (C) 1996 McCAT System Group.
     4      Copyright (C) 1996 ACAPS Benchmark Administrator
     5                         benadmin@acaps.cs.mcgill.ca
     6  
     7      This program is free software; you can redistribute it and/or modify
     8      it provided this copyright notice is maintained.
     9  
    10      This program is distributed in the hope that it will be useful,
    11      but WITHOUT ANY WARRANTY; without even the implied warranty of
    12      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
    13  ****/
    14  
    15  #include <stdio.h>
    16  #include <stdlib.h>
    17  #include <string.h>
    18  #include <math.h>
    19  
    20  #define DBL_EPSILON 0x1p-52
    21  
    22  void  *allocvector(size_t size) 
    23  {
    24    void *V;
    25  
    26    if ( (V = (void *) malloc((size_t) size)) == NULL ) {
    27      fprintf(stderr, "Error: couldn't allocate V in allocvector.\n");
    28      exit(2);
    29    }
    30    memset(V,0,size);
    31    return V;
    32  }
    33  
    34  void dallocvector(int n, double **V)
    35  {
    36    *V = (double *) allocvector((size_t) n*sizeof(double));
    37  }
    38  
    39  
    40  #define FUDGE  (double) 1.01
    41  
    42  
    43  
    44  int sturm(int n, double c[], double b[], double beta[], double x)
    45  
    46  /**************************************************************************
    47  
    48  Purpose:
    49  ------------
    50    Calculates the sturm sequence given by
    51  
    52      q_1(x)  =  c[1] - x
    53  
    54      q_i(x)  =  (c[i] - x) - b[i]*b[i] / q_{i-1}(x)
    55  
    56    and returns a(x) = the number of negative q_i. a(x) gives the number
    57    of eigenvalues smaller than x of the symmetric tridiagonal matrix
    58    with diagonal c[0],c[1],...,c[n-1] and off-diagonal elements
    59    b[1],b[2],...,b[n-1].
    60  
    61  
    62  Input parameters:
    63  ------------------------
    64    n :
    65           The order of the matrix.
    66  
    67    c[0]..c[n-1] : 
    68          An n x 1 array giving the diagonal elements of the tridiagonal matrix.
    69  
    70    b[1]..b[n-1] :   
    71          An n x 1 array giving the sub-diagonal elements. b[0] may be 
    72          arbitrary but is replaced by zero in the procedure.
    73  
    74    beta[1]..beta[n-1] :
    75           An n x 1 array giving the square of the  sub-diagonal elements, 
    76           i.e. beta[i] = b[i]*b[i]. beta[0] may be arbitrary but is replaced by
    77           zero in the procedure.
    78  
    79    x :
    80           Argument for the Sturm sequence.
    81  
    82  
    83  Returns:
    84  ------------------------
    85    integer a = Number of eigenvalues of the matrix smaller than x.
    86  
    87  
    88  Notes:
    89  ------------------------
    90  On SGI PowerChallenge this function should be compiled with option
    91  "-O3 -OPT:IEEE_arithmetic=3" in order to activate the optimization 
    92  mentioned in the code below.
    93  
    94  
    95  **********************************************************************/
    96  
    97  {
    98    int i;
    99    int a;
   100    double q;
   101    
   102    a = 0;
   103    q = 1.0;
   104  
   105    for(i=0; i<n; i++) {
   106  
   107  #ifndef TESTFIRST
   108  
   109      if (q != 0.0) {
   110  
   111  #ifndef RECIPROCAL
   112        q =  (c[i] - x) - beta[i]/q; 
   113  #else 
   114    /* A potentially NUMERICALLY DANGEROUS optimizations is used here.
   115     * The previous statement should read:
   116     *         q = (c[i] - x) - beta[i]/q 
   117     * But computing the reciprocal might help on some architectures
   118     * that have multiply-add and/or reciprocal instuctions.
   119     */
   120       iq = 1.0/q;  
   121        q =  (c[i] - x) - beta[i]*iq; 
   122  #endif
   123  
   124      }
   125      else {
   126        q = (c[i] - x) - fabs(b[i])/DBL_EPSILON;  
   127      }
   128  
   129      if (q < 0)
   130        a = a + 1;    
   131    }
   132  
   133  #else
   134      
   135      if (q < 0) {
   136        a = a + 1;    
   137        
   138  #ifndef RECIPROCAL
   139        q =  (c[i] - x) - beta[i]/q; 
   140  #else 
   141        /* A potentially NUMERICALLY DANGEROUS optimizations is used here.
   142         * The previous statement should read:
   143         *         q = (c[i] - x) - beta[i]/q 
   144         * But computing the reciprocal might help on some architectures
   145         * that have multiply-add and/or reciprocal instuctions.
   146         */
   147        iq = 1.0/q;  
   148        q =  (c[i] - x) - beta[i]*iq; 
   149  #endif
   150        
   151      }
   152      else if (q > 0.0) {
   153  #ifndef RECIPROCAL
   154        q =  (c[i] - x) - beta[i]/q; 
   155  #else 
   156        /* A potentially NUMERICALLY DANGEROUS optimizations is used here.
   157         * The previous statement should read:
   158         *         q = (c[i] - x) - beta[i]/q 
   159         * But computing the reciprocal might help on some architectures
   160         * that have multiply-add and/or reciprocal instuctions.
   161         */
   162        iq = 1.0/q;  
   163        q =  (c[i] - x) - beta[i]*iq; 
   164  #endif
   165      }
   166      else {
   167        q = (c[i] - x) - fabs(b[i])/DBL_EPSILON;  
   168      }
   169      
   170    }
   171    if (q < 0)
   172      a = a + 1;    
   173  #endif
   174  
   175    return a;
   176  }
   177    
   178  
   179  
   180  
   181  void dbisect(double c[], double b[], double beta[], 
   182  	     int n, int m1, int m2, double eps1, double *eps2, int *z,  
   183  	     double x[])
   184  
   185  
   186  /**************************************************************************
   187  
   188  Purpose:
   189  ------------
   190   
   191    Calculates eigenvalues lambda_{m1}, lambda_{m1+1},...,lambda_{m2} of
   192    a symmetric tridiagonal matrix with diagonal c[0],c[1],...,c[n-1]
   193    and off-diagonal elements b[1],b[2],...,b[n-1] by the method of
   194    bisection, using Sturm sequences.
   195  
   196  
   197    Input parameters:
   198  ------------------------
   199  
   200    c[0]..c[n-1] : 
   201          An n x 1 array giving the diagonal elements of the tridiagonal matrix.
   202  
   203    b[1]..b[n-1] :   
   204          An n x 1 array giving the sub-diagonal elements. b[0] may be 
   205          arbitrary but is replaced by zero in the procedure.
   206  
   207    beta[1]..beta[n-1] :
   208           An n x 1 array giving the square of the  sub-diagonal elements, 
   209           i.e. beta[i] = b[i]*b[i]. beta[0] may be arbitrary but is replaced by
   210           zero in the procedure.
   211  
   212    n :
   213           The order of the matrix.
   214  
   215    m1, m2 : 
   216           The eigenvalues lambda_{m1}, lambda_{m1+1},...,lambda_{m2} are 
   217           calculated (NB! lambda_1 is the smallest eigenvalue).
   218           m1 <= m2must hold otherwise no eigenvalues are computed.
   219           returned in x[m1-1],x[m1],...,x[m2-1]
   220   
   221    eps1 :
   222           a quantity that affects the precision to which eigenvalues are 
   223           computed. The bisection is continued as long as 
   224           x_high - x_low >  2*DBL_EPSILON(|x_low|  + |x_high|) + eps1       (1)
   225           When (1) is no longer satisfied, (x_high + x_low)/2  gives the
   226           current eigenvalue lambda_k. Here DBL_EPSILON (constant) is
   227           the machine accuracy, i.e. the smallest number such that
   228           (1 + DBL_EPSILON) > 1.
   229  
   230    Output parameters:
   231  ------------------------
   232  
   233    eps2 :
   234          The overall bound  for the error in any eigenvalue.
   235    z :
   236          The total number of iterations to find all eigenvalues.
   237    x : 
   238          The array  x[m1],x[m1+1],...,x[m2] contains the computed eigenvalues.
   239  
   240  **********************************************************************/
   241  {
   242    int i;
   243    double h,xmin,xmax;
   244    beta[0]  = b[0] = 0.0; 
   245  
   246  
   247    /* calculate Gerschgorin interval */
   248    xmin = c[n-1] - FUDGE*fabs(b[n-1]);
   249    xmax = c[n-1] + FUDGE*fabs(b[n-1]);
   250    for(i=n-2; i>=0; i--) { 
   251      h = FUDGE*(fabs(b[i]) + fabs(b[i+1]));
   252      if (c[i] + h > xmax)  xmax = c[i] + h;
   253      if (c[i] - h < xmin)  xmin = c[i] - h;
   254    }
   255  
   256    /*  printf("xmax = %lf  xmin = %lf\n",xmax,xmin); */
   257  
   258    /* estimate precision of calculated eigenvalues */  
   259    *eps2 = DBL_EPSILON * (xmin + xmax > 0 ? xmax : -xmin);
   260    if (eps1 <= 0)
   261      eps1 = *eps2;
   262    *eps2 = 0.5 * eps1 + 7 * *eps2;
   263    { int a,k;
   264      double x1,xu,x0;
   265      double *wu; 
   266  
   267      if( (wu = (double *) calloc(n+1,sizeof(double))) == NULL) {
   268        fputs("bisect: Couldn't allocate memory for wu",stderr);
   269        exit(1);
   270      }
   271  
   272      /* Start bisection process  */
   273      x0 = xmax;
   274      for(i=m2; i >= m1; i--) {
   275        x[i] = xmax;
   276        wu[i] = xmin;
   277      }
   278      *z = 0;
   279      /* loop for the k-th eigenvalue */
   280      for(k=m2; k>=m1; k--) {
   281        xu = xmin;
   282        for(i=k; i>=m1; i--) {
   283  	if(xu < wu[i]) {
   284  	  xu = wu[i];
   285  	  break;
   286  	}
   287        }
   288        if (x0 > x[k])
   289  	x0 = x[k];
   290  
   291        x1 = (xu + x0)/2;
   292        while ( x0-xu > 2*DBL_EPSILON*(fabs(xu)+fabs(x0))+eps1 ) {	
   293  	*z = *z + 1;
   294  	
   295  	/* Sturms Sequence  */
   296  
   297         	a = sturm(n,c,b,beta,x1); 
   298  
   299  	/* Bisection step */
   300  	if (a < k) {
   301  	  if (a < m1) 
   302  	    xu = wu[m1] = x1;
   303  	  else {
   304  	    xu = wu[a+1] = x1;
   305  	    if (x[a] > x1) x[a] = x1;
   306  	  }
   307  	}
   308  	else {
   309  	  x0 = x1;
   310  	}	
   311  	x1 = (xu + x0)/2;	
   312        }
   313        x[k] = (xu + x0)/2;
   314      }
   315      free(wu);
   316    }
   317  }     
   318  
   319  void test_matrix(int n, double *C, double *B)
   320  /* Symmetric tridiagonal matrix with diagonal
   321  
   322       c_i = i^4,  i = (1,2,...,n)
   323  
   324       and off-diagonal elements
   325  
   326       b_i = i-1,    i = (2,3,...n).
   327       It is possible to determine small eigenvalues of this matrix, with the
   328       same relative error as for the large ones. 
   329  */
   330  {
   331    int i;
   332      
   333    for(i=0; i<n; i++) {
   334      B[i] = (double) i;
   335      C[i] = (double ) (i+1)*(i+1);
   336      C[i] = C[i] * C[i];
   337    }
   338  }
   339  
   340  
   341  int main(int argc,char *argv[])
   342  {
   343    int rep,n,k,i,j;
   344    double eps,eps2;
   345    double *D,*E,*beta,*S;
   346  
   347    rep = 1;
   348    n = 500;
   349    eps = 2.2204460492503131E-16;
   350  
   351    dallocvector(n,&D);
   352    dallocvector(n,&E);
   353    dallocvector(n,&beta);
   354    dallocvector(n,&S);  
   355    
   356    for (j=0; j<rep; j++) {
   357      test_matrix(n,D,E);
   358      
   359      for (i=0; i<n; i++) {
   360        beta[i] = E[i] * E[i];
   361        S[i] = 0.0;
   362      }
   363      
   364      E[0] = beta[0] = 0;  
   365      dbisect(D,E,beta,n,1,n,eps,&eps2,&k,&S[-1]);    
   366      
   367    }
   368    
   369    for(i=1; i<20; i++)
   370      printf("%5d %.5e\n",i+1,S[i]); 
   371  
   372    printf("eps2 = %.5e,  k = %d\n",eps2,k);
   373  
   374    return 0;
   375  }
   376  
   377