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