harp3fastsim-ok.c Source File

00001 /* 
00002  *
00003  * HARP3FASTSIM.C : harp3 simulater which performs the operations all using
00004  * IEEE-754 arithmetics
00005  *
00006  */
00007 
00008 #include <stdio.h>
00009 #include <sys/types.h>
00010 #ifdef INTERNAL_OUT
00011 #   undef INTERNAL_OUT
00012 #   define  INTERNAL_OUT 1
00013 #else
00014 #   define  INTERNAL_OUT 0
00015 #endif
00016 
00017 unsigned int h3wait_();
00018 
00019 
00020 #define harpreal double
00021 
00022 #define NMAX 1000
00023 #define IPMAX 3
00024 #define JPMAX 26
00025 #define BUFMAX  50
00026 #define NDIM 3
00027 
00028 #define NNBMAX 1024
00029 
00030 int npipes = IPMAX;
00031 double xjmem[NMAX][NDIM];
00032 harpreal vjmem[NMAX][NDIM];
00033 double xjpred[NMAX][NDIM];
00034 harpreal vjpred[NMAX][NDIM];
00035 harpreal ajmem[NMAX][NDIM];
00036 harpreal jjmem[NMAX][NDIM];
00037 harpreal mjmem[NMAX];
00038 double tjmem[NMAX];
00039 
00040 double xjbuf[BUFMAX][JPMAX][NDIM];
00041 harpreal vjbuf[BUFMAX][JPMAX][NDIM];
00042 harpreal ajbuf[BUFMAX][JPMAX][NDIM];
00043 harpreal jjbuf[BUFMAX][JPMAX][NDIM];
00044 harpreal mjbuf[BUFMAX][JPMAX];
00045 double tjbuf[BUFMAX][JPMAX];
00046 int indexbuf[BUFMAX][JPMAX];
00047 int njbuf[BUFMAX];
00048 double timem;
00049 
00050 double ximem[IPMAX][NDIM];
00051 harpreal vimem[IPMAX][NDIM];
00052 double hmem[IPMAX];
00053 double epsmem[IPMAX];
00054 
00055 int nnb[IPMAX];
00056 int nb[IPMAX][NNBMAX];
00057 static int debug_level = 0;
00058 void h3open_()
00059 {
00060     
00061 }
00062 
00063 void h3close_()
00064 {
00065 }
00066 
00067 int h3npipe_()
00068 {
00069     return npipes;
00070 }
00071 
00072 void h3setnboards_(nb)
00073     int * nb;
00074 {
00075     
00076 }
00077 
00078 int h3getnboards_()
00079 {
00080     return 1;
00081     /* simulator always return 1 for now... */
00082 }
00083 
00084 unsigned int h3wait_()
00085 {
00086     return 0;
00087 }
00088 
00089 void h3setti_(ti)
00090     double * ti;
00091 {
00092     timem= *ti;
00093 }
00094 
00095 int h3jpmax_()
00096 {
00097 
00098     return JPMAX;
00099 }
00100 
00101 
00102 #ifdef G77
00103 #define h3jpdma_indirect_ h3jpdma_indirect__
00104 #define h3mjpdma_indirect_ h3mjpdma_indirect__
00105 #define h3mjpdma_start_ h3mjpdma_start__
00106 #define h3mjpdma_flush_ h3mjpdma_flush__
00107 #endif
00108 
00109 void h3setmode_()
00110 {}
00111 
00112 
00113 void h3jpdma_indirect_(int *,int*,double[][3],double[][3],double[][3],
00114                   double[][3],double*,double* ,int*);
00115 
00116 void h3mjpdma_indirect_(int *,int*,double[][3],double[][3],double[][3],
00117                   double[][3],double*,double* ,int*,int*);
00118 void h3mjpdma_start_(int *);
00119 void h3mjpdma_flush_();
00120 
00121 void h3jpdma_indirect_(nj,hostindex,xj,vj,aj,jj,mj,tj,mode)
00122     int * nj;
00123     int hostindex[];
00124     double xj[][3];
00125     double vj[][3];
00126     double aj[][3];
00127     double jj[][3];
00128     double mj[];
00129     double tj[];
00130     int *mode;
00131 {
00132     int zero = 0;
00133     h3mjpdma_indirect_(nj,hostindex,xj,vj,aj,jj,mj,tj,mode,&zero);
00134     h3mjpdma_start_(&zero);
00135 }
00136 
00137 
00138 static int buff_used = 0;
00139 void h3mjpdma_indirect_(nj,hostindex,xj,vj,aj,jj,mj,tj,mode,buff_id)
00140     int * nj;
00141     int hostindex[];
00142     double xj[][3];
00143     double vj[][3];
00144     double aj[][3];
00145     double jj[][3];
00146     double mj[];
00147     double tj[];
00148     int *mode;
00149     int * buff_id;
00150 {
00151     register int i;
00152     int bufno;
00153     if (debug_level) fprintf(stderr,"(h3mjpdma) called %d %d\n",*nj, *mode);
00154     if(*nj > JPMAX){
00155         fprintf(stderr,"(h3jpdma) Too large nj (%d > %d\n",nj, JPMAX);
00156         exit(-1);
00157     }
00158     if(*buff_id >= 0){
00159         bufno = *buff_id;
00160     }else{
00161         bufno = buff_used;
00162         buff_used ++;
00163     }
00164     njbuf[bufno] = *nj;
00165     for(i=0;i<*nj;i++){
00166         int k;
00167         int hostid;
00168         hostid = hostindex[i]-1;
00169         for(k=0;k<NDIM;k++){
00170             xjbuf[bufno][i][k] = xj[hostid][k];
00171             vjbuf[bufno][i][k] = vj[hostid][k];
00172             ajbuf[bufno][i][k] = aj[hostid][k];
00173             jjbuf[bufno][i][k] = jj[hostid][k];
00174         }
00175         tjbuf[bufno][i] = tj[hostid];
00176         mjbuf[bufno][i] = mj[hostid];
00177         indexbuf[bufno][i] = hostid;
00178     }
00179 }
00180 
00181 void h3mjpdma_start_(buff_id)
00182     int * buff_id;
00183 {
00184     register int i,k, j, buff;
00185     if (debug_level> 0) fprintf(stderr,"(h3mjpdma_start) called %d \n",*buff_id);
00186     buff = *buff_id;
00187     for(i=0;i<njbuf[buff];i++){
00188         int k;
00189         int j;
00190         j= indexbuf[buff][i];
00191         if (debug_level > 2){
00192             printf("jpdma_start %d %d %d %f\n",
00193                    i, j, buff, xjbuf[buff][i][0]);
00194         }
00195         for(k=0;k<NDIM;k++){
00196             xjmem[j][k] = xjbuf[buff][i][k];
00197             vjmem[j][k] = vjbuf[buff][i][k];
00198             ajmem[j][k] = ajbuf[buff][i][k];
00199             jjmem[j][k] = jjbuf[buff][i][k];
00200         }
00201         tjmem[j]= tjbuf[buff][i];
00202         mjmem[j]= mjbuf[buff][i];
00203     }
00204 
00205 }
00206 
00207 void h3jpdma_flush_()
00208 {
00209     int i;
00210     for(i=0;i<buff_used; i++){
00211         h3mjpdma_start_(&i);
00212     }
00213     buff_used = 0;
00214 }
00215 
00216 void h3jpdma_reset_()
00217 {
00218     buff_used = 0;
00219 }
00220 
00221 static int njset;
00222 h3calc_firsthalf_(nj,ni,xi,vi,eps2,h2)
00223     int * nj;
00224     int * ni;
00225     double xi[][3];
00226     double vi[][3];
00227     double eps2[];
00228     double h2[];
00229 {
00230     int i,j,k;
00231     njset = *nj;
00232     if(*ni > IPMAX){
00233         fprintf(stderr,"(h3ipdma) Too large ni (%d > %d\n",ni, IPMAX);
00234         exit(-1);
00235     }
00236     for (i=0;i < *ni; i++){
00237         for(k=0;k<3;k++)ximem[i][k] = xi[i][k];
00238         for(k=0;k<3;k++)vimem[i][k] = vi[i][k];
00239         hmem[i] = 2*h2[i];
00240         epsmem[i] = 2*eps2[i];
00241         /* factor of two just to mimic HARP chips...*/
00242     }
00243     if (debug_level > 0){
00244         for (i=0; i< *ni; i+= 2){
00245             if ((hmem[i] != hmem[i+1]) ||(hmem[i] != hmem[i+1])){
00246                 fprintf(stderr,"h3calc_firsthalf, eps/h2 error for %d\n", i);
00247             }
00248         }
00249     }
00250     for (j=0;j< *nj; j++){
00251         double s, s1, s2;
00252         s = timem - tjmem[j];
00253         s1 = 1.5*s;
00254         s2 = 2.0*s;
00255         for (k=0;k<3;k++){
00256             xjpred[j][k]= ((jjmem[j][k]*s + ajmem[j][k])*s + vjmem[j][k])*s
00257             + xjmem[j][k];
00258             vjpred[j][k]= (jjmem[j][k]*s1 + ajmem[j][k])*s2 + vjmem[j][k];
00259             if (debug_level > 2){
00260                 printf("predict %d %d %22.14g %22.14g\n",
00261                        j,k,xjpred[j][k], vjpred[j][k]);
00262             }
00263         }
00264     }
00265 }
00266 
00267 h3calc_lasthalf_(ni,acc,jerk,pot)
00268     int * ni;
00269     double acc[][3];
00270     double jerk[][3];
00271     double pot[];
00272 {
00273     register int i, j, k;
00274     if(*ni > IPMAX){
00275         fprintf(stderr,"(h3ipdma) Too large ni (%d > %d\n",ni, IPMAX);
00276         exit(-1);
00277     }
00278     if (debug_level) fprintf(stderr,"lasthalf ni=%d\n",*ni);
00279     for (i=0; i< *ni; i++){
00280         double sqrt();
00281         register double dx[3], dv[3], dr2, drdv, dr2i, dr3i;
00282         for(k=0;k<3 ; k++){
00283             acc[i][k] = 0.0;
00284             jerk[i][k] = 0.0;
00285         }
00286         nnb[i] = 0;
00287         pot[i] = 0.0;
00288         for(j=0;j<njset;j++){
00289             dr2 = epsmem[i];
00290             drdv = 0.0;
00291             for(k=0;k<3 ; k++){
00292                 dx[k] = xjpred[j][k] - ximem[i][k];
00293                 dv[k] = vjpred[j][k] - vimem[i][k];
00294                 dr2 +=  dx[k]*dx[k];
00295                 drdv += dx[k]*dv[k];
00296                 if (debug_level > 2){
00297                     fprintf(stderr,"rijk %d %d %d %g %21.14g %21.14g\n",i, j, k, dx[k],
00298                        xjpred[j][k],ximem[i][k]);
00299                 }
00300             }
00301             if (dr2 < hmem[i]){
00302                 nb[i][nnb[i]] = j+1;
00303                 nnb[i]++;
00304             }
00305             if (debug_level > 1){
00306                 fprintf(stderr, "i, j, nnb[i], dr2, h = %d %d %d %g %g\n",
00307                        i, j, nnb[i], dr2, hmem[i]);
00308             }
00309             if (dr2 > 0.0){
00310                 dr2i = 1.0/dr2;
00311                 dr3i = mjmem[j]*dr2i*sqrt(dr2i);
00312                 drdv = 3.*drdv*dr2i;
00313                 for(k=0;k<3 ; k++){
00314                     acc[i][k] +=  dx[k]*dr3i;
00315                     jerk[i][k] += (dv[k] - dx[k]*drdv)*dr3i;
00316                     if (debug_level > 1){
00317                         fprintf(stderr,"force i,j,k %d %d %d %20.10g \n",i, j, k, acc[i][k]);
00318                     }
00319                 }
00320             }else{
00321                 if (debug_level > 1){
00322                     fprintf(stderr,"force skip i,j %d %d  \n",i, j);
00323                 }
00324             }
00325         }
00326     }   
00327 }
00328 
00329 h3calc_(nj,ni,xi,vi,eps2,h2, acc, jerk, pot)
00330     int * nj;
00331     int * ni;
00332     double xi[][3];
00333     double vi[][3];
00334     double eps2[];
00335     double h2[];
00336     double acc[][3];
00337     double jerk[][3];
00338     double pot[];
00339 
00340 {
00341     h3calc_firsthalf_(nj,ni,xi,vi,eps2,h2);
00342     h3calc_lasthalf_(ni,acc,jerk,pot);
00343 }
00344 
00345 void h3nbread_(nboards)
00346     int * nboards;
00347 {
00348 
00349 }
00350 #ifdef G77
00351 # define set_debug_level_ set_debug_level__
00352 #endif
00353 void set_debug_level_(int * level)
00354 {
00355     debug_level = *level;
00356 }
00357 
00358 void h3setdebuglevel_(int * level)
00359 {
00360     set_debug_level_(level);
00361 }
00362 int compare(i,j)
00363     int *i;
00364     int *j;
00365 {
00366     return *i - *j;
00367 }
00368 
00369 int h3nblist_(board, chip, nblist)
00370     int * board;
00371     int * chip;
00372     int * nblist;
00373 {
00374 
00375     int i;
00376     for(i=0; i< nnb[*chip]; i++){
00377         nblist[i] = nb[*chip][i];
00378     }
00379     return nnb[*chip];
00380 }
00381 
00382