diff --git a/src/amuse/community/bhtree/src/BHTC.C b/src/amuse/community/bhtree/src/BHTC.C index 5e23003b92..35d51dbe18 100644 --- a/src/amuse/community/bhtree/src/BHTC.C +++ b/src/amuse/community/bhtree/src/BHTC.C @@ -5,7 +5,7 @@ // // Change in the calling format for apply_vf, from function name to // add explicit address operator. This was necessary to pass g++ 2.8.1 -// +// // Version 1.1 Jun Makino : take parameters from command line argumensts // Coding for 1.1 started on 1998/12/31 // Version 1.0 Jun Makino : First publisized version, parameters @@ -50,7 +50,7 @@ int main(int argc, char ** argv) extern char *poptarg; int c; - char* param_string = "i:o:d:D:T:e:t:n:w:cx:v:s:S:h"; + char const * param_string = "i:o:d:D:T:e:t:n:w:cx:v:s:S:h"; foname[0] = '?'; foname[1] = '\0'; real eps = 0.025; @@ -103,7 +103,7 @@ int main(int argc, char ** argv) break; case 'S': vel_scale = atof(poptarg); break; - case 'h': + case 'h': cerr << "list of options\n"; cerr << "-i name of snapshot input file (no default)\n"; cerr << "-o name of snapshot output file (default: no output)\n"; @@ -141,7 +141,7 @@ int main(int argc, char ** argv) cout << "n= " << pb.n << endl; pb.use_self_gravity = 1; fsnapin.close(); - if (snapout) fsnapout.open(foname,ios::out); + if (snapout) fsnapout.open(foname,ios::out); PR(eps); PR(theta); PR(ncrit); cout << "eps= " << eps << " theta=" << theta << " ncrit=" < l*0.5) return 1; } return 0; } - - + + int bhnode::sanity_check() { int i; @@ -343,7 +343,7 @@ int bhnode::sanity_check() return iret; } -#ifdef SPH +#ifdef SPH void bhnode::set_hmax_for_sph() { int i; @@ -453,7 +453,7 @@ int check_and_set_nbl(bhnode * p1,bhnode * p2) if (p1->isleaf || (p2->isleaf == 0) && p2->l > p1->l){ // // p1 is already leaf, or both are node but p2 is larger. - // go down p2 by + // go down p2 by for(int i = 0; i<8;i++) iret |= check_and_set_nbl(p1, p2->child[i]); }else{ @@ -473,7 +473,7 @@ int check_and_set_nbl(bhnode * p1,bhnode * p2) // // here, comparison of key guarantee that a pair // is evaluated only once - + iret |=check_and_set_nbl(*((bpi+i)->get_rp()), *((bpj+j)->get_rp())); } @@ -508,19 +508,19 @@ void real_system::setup_tree() bhnode * btmp = bn+1; int heap_remainder = bnsize-1; BHlong key = 0; - + bn->create_tree_recursive(btmp,heap_remainder,key, default_key_length, ncrit_for_tree); - //cello: last argument was 12 for some reason - - + //cello: last argument was 12 for some reason + + set_cm_quantities_for_default_tree(); // PR(bnsize); PRL(heap_remainder); // PRL(bn->sanity_check()); } - -#ifdef SPH + +#ifdef SPH int sph_system::set_nnb_using_tree() { setup_tree(); @@ -533,7 +533,7 @@ int sph_system::set_nnb_using_tree() } #endif -void accumulate_force_from_point(vec dx, real r2, real eps2, +void accumulate_force_from_point(vec dx, real r2, real eps2, vec & acc, real & phi, real jmass) @@ -562,7 +562,7 @@ void print_tree_counters() } void calculate_force_from_interaction_list(const vec & pos, - real eps2, + real eps2, vec & acc, real & phi, vec * poslist, @@ -596,11 +596,11 @@ void bhnode::accumulate_force_from_tree(vec & ipos, real eps2, real theta2, int i; if (isleaf){ bhparticle * bp = bpfirst; - + for(i = 0; i < nparticle; i++){ vec dx = (bp+i)->get_rp()->get_pos()-ipos; real r2 = dx*dx; - + // Added by Steve (8/07). // Turned off by AVE (12/12) if (r2 < r_nn_2) { @@ -611,12 +611,12 @@ void bhnode::accumulate_force_from_tree(vec & ipos, real eps2, real theta2, nn_ptr = (bp+i)->get_rp(); } } - - + + accumulate_force_from_point(dx, r2, eps2, acc, phi, (bp+i)->get_rp()->get_mass()); total_interactions += 1; - + } }else{ for(i=0;i<8;i++){ @@ -636,7 +636,7 @@ void bhnode::add_to_interaction_list(bhnode & dest_node, real theta2, int list_max, int & first_leaf) { -#if 0 +#if 0 real r2 = separation_squared(this, &dest_node); if (r2*theta2 > l*l){ #else @@ -650,7 +650,7 @@ void bhnode::add_to_interaction_list(bhnode & dest_node, real theta2, cerr << "List length exceeded\n"; exit(1); } - + }else{ int i; if (isleaf || (this == (&dest_node))){ @@ -724,11 +724,11 @@ void calculate_force_from_interaction_list_using_grape4(vec * pos_list, real * m #ifdef GPU extern "C" { - + int g6_open_(int *clusterid); int g6_close_(int *clusterid); - + int g6_set_j_particle_(int *cluster_id, int *address, int *index, @@ -736,26 +736,26 @@ extern "C" { double *mass, double k18[3], double j6[3], double a2[3], double v[3], double x[3]); - + int g6_set_ti_(int *id, double *ti); - - + + void g6calc_firsthalf_(int *cluster_id, int *nj, int *ni, - int index[], + int index[], double xi[][3], double vi[][3], double aold[][3], double j6old[][3], - double phiold[3], + double phiold[3], double *eps2, double h2[]); - - + + int g6calc_lasthalf_(int *cluster_id, int *nj, int *ni, - int index[], + int index[], double xi[][3], double vi[][3], double *eps2, double h2[], - double acc[][3], double jerk[][3], double pot[]); - + double acc[][3], double jerk[][3], double pot[]); + int g6_reset_(int* cluster_id); int g6_npipes_(); } @@ -774,18 +774,18 @@ void calculate_force_from_interaction_list_using_grape6(vec * pos_list, real * m } else { //g6_reset_(&clusterid); } - + double ti = 0.0; double zero = 0.0; g6_set_ti_(&clusterid,&ti); int max_ni = g6_npipes_(); - + nisum += ni; tree_walks += 1; total_interactions += ((real)ni)*list_length; //accel_by_harp3_separate_noopen_(&ni,pos_list+first_leaf, &list_length,pos_list, mass_list, // acc_list, phi_list, &eps2); -// +// call_count += ni; double k18[3]; double j6[3]; @@ -832,32 +832,32 @@ void calculate_force_from_interaction_list_using_grape6(vec * pos_list, real * m jerkout[i][k]=0.0; } zero1[i] = 0; - index[i] =i + first_leaf; - + index[i] =i + first_leaf; + } - + for(int pi = 0; pi < ni; pi+=max_ni) { int actualni = pi + max_ni > ni ? ni - pi : max_ni; - + g6calc_firsthalf_( - &clusterid, - &list_length, - &actualni, - index + pi, + &clusterid, + &list_length, + &actualni, + index + pi, positions + pi, - zero3, zero3, - zero3, - zero1, + zero3, + zero3, + zero1, &eps2, zero1 - ); - + ); + g6calc_lasthalf_( - &clusterid, - &list_length, - &actualni, + &clusterid, + &list_length, + &actualni, index + pi, positions + pi, zero3, @@ -895,7 +895,7 @@ void bhnode::evaluate_gravity_using_tree_and_list(bhnode & source_node, static real mass_list[list_max]; static vec pos_list[list_max]; real epsinv = 1.0/sqrt(eps2); -#ifdef HARP3 +#ifdef HARP3 static vec * acc_list = NULL; static real * phi_list = NULL; if (acc_list == NULL){ @@ -903,7 +903,7 @@ void bhnode::evaluate_gravity_using_tree_and_list(bhnode & source_node, phi_list = new real[ncrit + 100]; } #endif -#ifdef GPU +#ifdef GPU static vec * acc_list = NULL; static real * phi_list = NULL; if (acc_list == NULL){ @@ -939,7 +939,7 @@ void bhnode::evaluate_gravity_using_tree_and_list(bhnode & source_node, exit(1); } bhparticle * bp = bpfirst; -#ifndef HARP3 +#ifndef HARP3 #ifndef GPU for(int i = 0; i < nparticle; i++){ real_particle * p = (bp+i)->get_rp(); @@ -967,7 +967,7 @@ void bhnode::evaluate_gravity_using_tree_and_list(bhnode & source_node, p->set_acc_gravity(acc_list[i]); p->set_phi_gravity(phi_list[i] + p->get_mass()*epsinv); } -#endif +#endif } } @@ -995,7 +995,7 @@ void real_particle::calculate_gravity_using_tree(real eps2, real theta2) bn->accumulate_force_from_tree(pos, eps2, theta2, acc_gravity, phi_gravity, index); nisum += 1; - + // AMUSE STOPPING CONDITIONS SUPPORT (AVE May 2010) //if(nn_index >= 0 && (COLLISION_DETECTION_BITMAP & enabled_conditions) ) { // compare indices to prevent two detections for one collision: (NdV Feb 2012) if(nn_index > my_index && (COLLISION_DETECTION_BITMAP & enabled_conditions) ) { @@ -1024,7 +1024,7 @@ vec real_system::calculate_gravity_at_point(vec pos, real eps2, real theta2) r_nn_2 = 0.0; // turn-off collision detection bn->accumulate_force_from_tree(pos, eps2, theta2, acc_gravity, phi_gravity, -1); - + return acc_gravity; } @@ -1038,7 +1038,7 @@ real real_system::calculate_potential_at_point(vec pos, real eps2, real theta2) r_nn_2 = 0.0; // turn-off collision detection bn->accumulate_force_from_tree(pos, eps2, theta2, acc_gravity, phi_gravity, -1); - + return phi_gravity; } @@ -1075,7 +1075,7 @@ int main() if((psph+i)->get_nnb() != (psphcopy+i)->get_nnb()){ cerr << "Neighbour count differs for "; PRL(i); err = 1; - + } if (err == 0){ for(int j = 0; (j< (psph+i)->get_nnb()) && (err == 0); j++){ @@ -1145,7 +1145,7 @@ void main() if((psph+i)->get_nnb() != (psphcopy+i)->get_nnb()){ cerr << "Neighbour count differs for "; PRL(i); err = 1; - + } if (err == 0){ for(int j = 0; (j< (psph+i)->get_nnb()) && (err == 0); j++){ @@ -1162,8 +1162,8 @@ void main() PRL(error); pb.use_self_gravity = 1; pb.eps2_for_gravity = 0.01; -#define COMPARISON_WITH_DIRECT -#ifdef COMPARISON_WITH_DIRECT +#define COMPARISON_WITH_DIRECT +#ifdef COMPARISON_WITH_DIRECT pb.calculate_uncorrected_gravity_direct(); copy_sph_particles(&pb, &pbcopy); psphcopy = pbcopy.get_particle_pointer(); @@ -1176,7 +1176,7 @@ void main() #endif - + cerr << "Tree force \n"; for(int j = 0; j<10; j++){ PRL(j); @@ -1187,7 +1187,7 @@ void main() } pb.apply_vf(real_particle::clear_acc_phi_gravity); bn->evaluate_gravity_using_tree_and_list(*bn,0.4,pb.eps2_for_gravity,1); -#ifdef COMPARISON_WITH_DIRECT +#ifdef COMPARISON_WITH_DIRECT real perrmax = 0; real ferrmax = 0; for(int i = 0; iget_phi_gravity(); vec acc = (psph+i)->get_acc_gravity(); - PR(i); PR(phi); PRL(acc); + PR(i); PR(phi); PRL(acc); } - -#endif - + +#endif + } #endif - + #ifdef TESTXX // @@ -1248,14 +1248,14 @@ void main() real_particle * psph = pb.get_particle_pointer(); real h0 = pow(1.0/n,0.33333); for(int i = 0; i<10; i++){ - + pb.apply_vf(real_particle::set_h, h0); pb.apply_vf(real_particle::clear_nnb); bn->set_hmax_for_sph(); // bn->dump(); PRL(check_and_set_nbl(bn, bn)); pb.apply_vf(real_particle::sort_nblist); - + } } #endif diff --git a/src/amuse/community/bhtree/src/pgetopt.C b/src/amuse/community/bhtree/src/pgetopt.C index e7a7843a1b..7b56b2e881 100644 --- a/src/amuse/community/bhtree/src/pgetopt.C +++ b/src/amuse/community/bhtree/src/pgetopt.C @@ -6,7 +6,7 @@ * version 2: Dec 1992 Piet Hut -- adopted to the new C++-based starlab * version 2.0A: Dec 1998 Jun Makino -- pskipopt() function added. *............................................................................. - * non-local function: + * non-local function: * pgetopt * pskipopt *............................................................................. @@ -26,7 +26,7 @@ * mumble -c -a -b 10 * mumble -b 10 -ca * mumble -acb 10 - * but the following versions are illegal, and will give error messages: + * but the following versions are illegal, and will give error messages: * mumble -a -b10 -c * mumble -ab10c * mumble -a -c -b @@ -53,15 +53,15 @@ char *poptarg; /*----------------------------------------------------------------------------- * pgetopt -- each call to pgetopt() returns the next option encountered - * on the command line, starting from the beginning. If the - * option occurs in the optstr string as well, the option - * itself is returned (as the int value of the character; + * on the command line, starting from the beginning. If the + * option occurs in the optstr string as well, the option + * itself is returned (as the int value of the character; * options are limited to one char). Otherwise, the character '?' - * is returned. If the end of the string is reached, the value + * is returned. If the end of the string is reached, the value * -1 is returned. If an option is followed by the character ':' - * in optstr , then a command line argument is expected, + * in optstr , then a command line argument is expected, * separated by spaces. If such an argument if found, the pointer - * poptarg is pointed at that string, so that the calling + * poptarg is pointed at that string, so that the calling * function can access the argument. If such an argument is not * found, an error message is given. * @@ -74,24 +74,26 @@ char *poptarg; */ static int argv_counter = 1; /* skip argv[0], the command name */ static int argv_offset = 0; -int pgetopt(int argc, char ** argv, char * optstr) +int pgetopt(int argc, char ** argv, char const * optstr) { int optstr_counter; char option_char; - + if (argv_counter >= argc) return(-1); /* signal that we've run out of options */ - if (argv_offset == 0) + if (argv_offset == 0) { if (argv[argv_counter][argv_offset] != '-') { cerr << "pgetopt: command line argument does not begin with -\n"; exit(1); } - else + else { argv_offset++; + } + } -#ifdef IN_STARLAB +#ifdef IN_STARLAB // We have a legal switch. First check to see if all we want to // know is the STARLAB version number. @@ -110,7 +112,7 @@ int pgetopt(int argc, char ** argv, char * optstr) { /* * prepare for the next call to pgetopt(): - */ + */ if (argv[argv_counter][++argv_offset] == '\0') { argv_counter++; @@ -138,7 +140,7 @@ int pgetopt(int argc, char ** argv, char * optstr) poptarg = argv[argv_counter]; /* * prepare for the next call to pgetopt(): - */ + */ argv_counter++; argv_offset = 0; } @@ -146,14 +148,14 @@ int pgetopt(int argc, char ** argv, char * optstr) { /* * prepare for the next call to pgetopt(): - */ + */ if (argv[argv_counter][++argv_offset] == '\0') { argv_counter++; argv_offset = 0; } } - + return(option_char); } diff --git a/src/amuse/community/bhtree/src/second.c b/src/amuse/community/bhtree/src/second.c index e23c713b0c..5b707569dc 100644 --- a/src/amuse/community/bhtree/src/second.c +++ b/src/amuse/community/bhtree/src/second.c @@ -10,7 +10,7 @@ * void second(&double): returns the cpu seconds * void cpumin(&double): returns the cpu minutes */ - + #ifdef _WIN32 #include @@ -21,16 +21,14 @@ void timer_init() time(&start_time); } -second(dtime) -double * dtime; +second(double * dtime) { long int timenow; time(&timenow); *dtime = timenow - start_time; } -cpumin(dmin) -double *dmin ; +cpumin(double * dmin) { double sec; second(&sec); @@ -60,7 +58,7 @@ static double tstart; static struct timeval timearg; static struct timezone zonearg; -#define RETURN_CPU +#define RETURN_CPU void timer_init() { @@ -89,38 +87,12 @@ void tminit_() timer_init(); } -cpumin(t) -double * t; -{ - xcpumin_(t); -} - -double cpusec() -{ - double t; - second(&t); - return t; -} -second(t) -double * t; -{ - second_(t); -} -xcpumin_(t) - double*t; -{ - double sec; - second_(&sec); - *t=sec/60.0; -} - -second_(t) -double *t; +void second_(double * t) { -#ifdef RETURN_CPU +#ifdef RETURN_CPU #ifdef SOLARIS struct tms buffer; - + if (times(&buffer) == -1) { printf("times() call failed\n"); exit(1); @@ -132,7 +104,7 @@ double *t; fprintf(stderr,"getrusage failed\n"); } *t = usage.ru_utime.tv_sec + usage.ru_utime.tv_usec*1e-6; -#endif +#endif #else #ifndef xxxx if(gettimeofday(&timearg,&zonearg)){ @@ -148,10 +120,34 @@ double *t; } *t = it0*0.01 - tstart; #endif -#endif +#endif } #endif +void xcpumin_(double * t) +{ + double sec; + second_(&sec); + *t=sec/60.0; +} + +void cpumin(double * t) +{ + xcpumin_(t); +} + +void second(double * t) +{ + second_(t); +} + +double cpusec() +{ + double t; + second(&t); + return t; +} + #ifdef TEST main() { @@ -172,10 +168,8 @@ typedef union double_and_int{ unsigned int iwork; COMPTYPE fwork; } DOUBLE_INT; -void c_assignbody_(pos,cpos,bsub) -register COMPTYPE pos[]; -register COMPTYPE cpos[]; -int * bsub; + +void c_assignbody_(COMPTYPE pos[], COMPTYPE cpos[], int * bsub) { register int k,l,k1,k2; register DOUBLE_INT tmpx,tmpy,tmpz ; @@ -209,7 +203,7 @@ int * bsub; tmpx.iwork,tmpy.iwork,tmpz.iwork, tmpx.iwork>>31,tmpy.iwork>>30,tmpz.iwork>>29, *bsub,k); -#endif +#endif #endif }