| 122 | |
| 123 | |
| 124 | int EMD_wrap_omp(int n1, int n2, double *X, double *Y, double *D, double *G, |
| 125 | double* alpha, double* beta, double *cost, uint64_t maxIter, int numThreads) { |
| 126 | // beware M and C are stored in row major C style!!! |
| 127 | |
| 128 | using namespace lemon_omp; |
| 129 | uint64_t n, m, cur; |
| 130 | |
| 131 | typedef FullBipartiteDigraph Digraph; |
| 132 | DIGRAPH_TYPEDEFS(Digraph); |
| 133 | |
| 134 | // Get the number of non zero coordinates for r and c |
| 135 | n=0; |
| 136 | for (int i=0; i<n1; i++) { |
| 137 | double val=*(X+i); |
| 138 | if (val>0) { |
| 139 | n++; |
| 140 | }else if(val<0){ |
| 141 | return INFEASIBLE; |
| 142 | } |
| 143 | } |
| 144 | m=0; |
| 145 | for (int i=0; i<n2; i++) { |
| 146 | double val=*(Y+i); |
| 147 | if (val>0) { |
| 148 | m++; |
| 149 | }else if(val<0){ |
| 150 | return INFEASIBLE; |
| 151 | } |
| 152 | } |
| 153 | |
| 154 | // Define the graph |
| 155 | |
| 156 | std::vector<uint64_t> indI(n), indJ(m); |
| 157 | std::vector<double> weights1(n), weights2(m); |
| 158 | Digraph di(n, m); |
| 159 | NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, (int) (n + m), n * m, maxIter, numThreads); |
| 160 | |
| 161 | // Set supply and demand, don't account for 0 values (faster) |
| 162 | |
| 163 | cur=0; |
| 164 | for (uint64_t i=0; i<n1; i++) { |
| 165 | double val=*(X+i); |
| 166 | if (val>0) { |
| 167 | weights1[ cur ] = val; |
| 168 | indI[cur++]=i; |
| 169 | } |
| 170 | } |
| 171 | |
| 172 | // Demand is actually negative supply... |
| 173 | |
| 174 | cur=0; |
| 175 | for (uint64_t i=0; i<n2; i++) { |
| 176 | double val=*(Y+i); |
| 177 | if (val>0) { |
| 178 | weights2[ cur ] = -val; |
| 179 | indJ[cur++]=i; |
| 180 | } |
| 181 | } |