#include "config.h" #ifdef USE_MPI #include #endif #include #include #include #include #include #include /* turn on a harness main() (bottom of the file) */ #define _JV_STANDALONE // used to 'simulate' 2D arrays using linear arrays #define sqrmat(x,y,a) ((y)*(a) + (x)) #define cpuClock() ((double)clock()/CLOCKS_PER_SEC) using namespace std; #include "emil.h" #include "linearAssignment.h" float LinearAssignment::getPairCost(int partId, int wellId){ float dist2, dx, dy, dz; Particle *p; Well *w; p = part[partId]; w = well[wellId]; /* just do the distance calc locally..easier than navigating the now-quite-bloated Hsc class structure, and also saves casting between doubles and floats. */ dx = (p->R[0] - w->R[0]); dy = (p->R[1] - w->R[1]); dz = (p->R[2] - w->R[2]); //it is possible to do this more efficiently if( periodic ) { while( dx >= halfBox[0] ){ dx -= box[0]; } while( dx < -1*halfBox[0] ){ dx += box[0]; } while( dy >= halfBox[1] ){ dy -= box[1]; } while( dy < -1*halfBox[1] ){ dy += box[1]; } while( dz >= halfBox[2] ){ dz -= box[2]; } while( dz < -1*halfBox[2] ){ dz += box[2]; } } dist2 = dx*dx + dy*dy + dz*dz; return( dist2 ); } // Adapted from C++ code by Rob Tyka // Adapted from PASCAL version in // Jonker, R. and Volgenant, A., 1987. A shortest augmenting // path algorithm for dense and sparse linear assignment problems. // Computing 38, pp. 325–340 double LinearAssignment::jonkerVolgenantAssign(){ bool found_unassigned; int i, i1, f0=0, lastfree, f, i0, k, f1; int j, j1, j2, end, last, low, up; int *pred = new int[N]; //predecessor-array for shortest path tree int *free = new int[N]; //unassigned rows (number f0, index f) int *col = new int[N]; //col array of columns, scanned int *match = new int[N]; int *d = new int[N]; // shortest path lengths float min, h, u1, u2, tmp1; float *v = new float[N]; // column cost float *u = new float[N]; // row cost for(i=0; i=0; j--) { min = cm_lazy(0,j); i1 = 0; for(i=1; i=u1) { u2 = h; j2 = j; } else { u2 = u1; u1 = h; j2 = j1; j1 = j; } } } i0 = wellToPart[j1]; if(u1 < u2){ v[j1] = v[j1] - (u2 - u1); }else{ if(i0>=0) { j1 = j2; i0 = wellToPart[j2]; } } partToWell[i] = j1; wellToPart[j1] = i; if(i0 >= 0){ if(u1 < u2){ free[--k] = i0; }else{ free[f0++] = i0; } } } } while(loopcnt < 2); // routine applied twice //augmentation for(f=0; f