1 #include "mrilib.h"
2
3 #define TEST_TIMER
4
main(int argc,char ** argv)5 int main( int argc , char **argv )
6 {
7 double a[9] , e[3] ; int ii ;
8 double aa[9] , ee[3] ;
9 double am[9] ;
10 double ss , smax=0.0 , dd , dmax=0.0 ;
11 double c1,c2,c3 ;
12
13 srand48((long)time(NULL)) ; symeig_forbid_23(1) ;
14
15 #ifdef TEST_TIMER
16 am[0] = 2.0*drand48() ;
17 am[4] = 2.0*drand48() ;
18 am[8] = 2.0*drand48() ;
19 am[1] = am[3] = drand48() ;
20 am[2] = am[6] = drand48() ;
21 am[5] = am[7] = drand48() ;
22 c1 = COX_cpu_time() ;
23 for( ii=0 ; ii < 1000000 ; ii++ ){
24 memcpy( a , am , sizeof(double)*9 ) ;
25 symeig_3( a , e , 1 ) ;
26 }
27 c2 = COX_cpu_time() ;
28 for( ii=0 ; ii < 1000000 ; ii++ ){
29 memcpy( aa , am , sizeof(double)*9 ) ;
30 symeig_double( 3 , aa , ee ) ;
31 }
32 c3 = COX_cpu_time() ;
33
34 fprintf(stderr,"CPU: New=%g Old=%g Ratio=%g\n",c2-c1 , c3-c2 , (c3-c2)/(c2-c1) ) ;
35
36 fprintf(stderr,"New #1: %10.5f [ %10.5f %10.5f %10.5f ]\n",
37 e[0], a[0],a[1],a[2] ) ;
38
39 fprintf(stderr,"Old #1: %10.5f [ %10.5f %10.5f %10.5f ]\n",
40 ee[0], aa[0],aa[1],aa[2] ) ;
41
42 fprintf(stderr,"New #2: %10.5f [ %10.5f %10.5f %10.5f ]\n",
43 e[1], a[3],a[4],a[5] ) ;
44
45 fprintf(stderr,"Old #2: %10.5f [ %10.5f %10.5f %10.5f ]\n",
46 ee[1], aa[3],aa[4],aa[5] ) ;
47
48 fprintf(stderr,"New #3: %10.5f [ %10.5f %10.5f %10.5f ]\n",
49 e[2], a[6],a[7],a[8] ) ;
50
51 fprintf(stderr,"Old #3: %10.5f [ %10.5f %10.5f %10.5f ]\n",
52 ee[2], aa[6],aa[7],aa[8] ) ;
53
54 #else
55
56 #if 0
57 for( ii=0 ; ii < 300000 ; ii++ ){
58 am[0] = drand48() ;
59 am[4] = drand48() ;
60 am[8] = drand48() ;
61 am[1] = am[3] = drand48() ;
62 am[2] = am[6] = drand48() ;
63 am[5] = am[7] = drand48() ;
64 memcpy( a , am , sizeof(double)*9 ) ;
65 symeig_3( a , e , 1 ) ;
66 memcpy( aa , am , sizeof(double)*9 ) ;
67 symeig_double( 3 , aa , ee ) ;
68 ss = fabs(ee[0]-e[0]) + fabs(ee[1]-e[1]) + fabs(ee[2]-e[2]) ;
69 if( ss > 1.e-6 ) fprintf(stderr,"*") ;
70 if( ss > smax ) smax = ss ;
71
72 dd = fabs( fabs(a[0]*aa[0] + a[1]*aa[1] + a[2]*aa[2])
73 +fabs(a[3]*aa[3] + a[4]*aa[4] + a[5]*aa[5])
74 +fabs(a[6]*aa[6] + a[7]*aa[7] + a[8]*aa[8]) - 3.0 ) ;
75
76 if( dd > 1.e-6 ) fprintf(stderr,"+") ;
77 if( dd > dmax ) dmax = dd ;
78 }
79 #else
80 for( ii=0 ; ii < 900000 ; ii++ ){
81 am[0] = drand48() ;
82 am[3] = drand48() ;
83 am[1] = am[2] = ( lrand48()%131 != 0 ) ? drand48() : 0.0 ;
84 memcpy( a , am , sizeof(double)*4 ) ;
85 symeig_2( a , e , 1 ) ;
86 memcpy( aa , am , sizeof(double)*4 ) ;
87 symeig_double( 2 , aa , ee ) ;
88 ss = fabs(ee[0]-e[0]) + fabs(ee[1]-e[1]) ;
89 if( ss > 1.e-6 ) fprintf(stderr,"* OLD=%g,%g NEW=%g,%g ",e[0],e[1],ee[0],ee[1]) ;
90 if( ss > smax ) smax = ss ;
91
92 dd = fabs( fabs(a[0]*aa[0] + a[1]*aa[1])
93 +fabs(a[2]*aa[2] + a[3]*aa[3]) - 2.0 ) ;
94
95 if( dd > 1.e-6 ) fprintf(stderr,"+") ;
96 if( dd > dmax ) dmax = dd ;
97 }
98 #endif
99 fprintf(stderr," smax=%g dmax=%g\n",smax,dmax) ;
100 #endif
101
102 exit(0) ;
103 }
104