modernc.org/ccgo/v3@v3.16.14/lib/testdata/CompCert-3.6/test/c/spectral.c (about) 1 /* -*- mode: c -*- 2 * 3 * The Great Computer Language Shootout 4 * http://shootout.alioth.debian.org/ 5 * 6 * Contributed by Sebastien Loisel 7 */ 8 9 #include <stdio.h> 10 #include <stdlib.h> 11 #include <math.h> 12 13 static inline double eval_A(int i, int j) { return 1.0/((i+j)*(i+j+1)/2+i+1); } 14 15 void eval_A_times_u(int N, const double u[], double Au[]) 16 { 17 int i,j; 18 for(i=0;i<N;i++) 19 { 20 Au[i]=0; 21 for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j]; 22 } 23 } 24 25 void eval_At_times_u(int N, const double u[], double Au[]) 26 { 27 int i,j; 28 for(i=0;i<N;i++) 29 { 30 Au[i]=0; 31 for(j=0;j<N;j++) Au[i]+=eval_A(j,i)*u[j]; 32 } 33 } 34 35 void eval_AtA_times_u(int N, const double u[], double AtAu[]) 36 { 37 double *v = malloc(N * sizeof(double)); 38 eval_A_times_u(N,u,v); 39 eval_At_times_u(N,v,AtAu); 40 free(v); 41 } 42 43 int main(int argc, char *argv[]) 44 { 45 int i; 46 int N = ((argc == 2) ? atoi(argv[1]) : 1000); 47 double * u, * v, vBv, vv; 48 u = malloc(N * sizeof(double)); 49 v = malloc(N * sizeof(double)); 50 for(i=0;i<N;i++) u[i]=1; 51 for(i=0;i<10;i++) 52 { 53 eval_AtA_times_u(N,u,v); 54 eval_AtA_times_u(N,v,u); 55 } 56 vBv=vv=0; 57 for(i=0;i<N;i++) { vBv+=u[i]*v[i]; vv+=v[i]*v[i]; } 58 printf("%0.9f\n",sqrt(vBv/vv)); 59 return 0; 60 } 61 62 /******** 63 build & benchmark results 64 65 BUILD COMMANDS FOR: spectralnorm.gcc 66 67 Fri Sep 15 20:48:08 PDT 2006 68 69 /usr/bin/gcc -pipe -Wall -O3 -fomit-frame-pointer -funroll-loops -march=pentium4 -lm spectralnorm.c -o spectralnorm.gcc_run 70 71 ================================================================= 72 COMMAND LINE (%A is single numeric argument): 73 74 spectralnorm.gcc_run %A 75 76 N=2500 77 78 PROGRAM OUTPUT 79 ============== 80 1.274224153 81 ********/