#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <assert.h>
#include <mpi.h>
#include "smut.h"
#include "stripesmut.h"
#include "auction.h"
#include "stripe_auction.h"
int comm_ngives = 0;
void
stripe_auction_shift (const struct gcsr_t *gA, double *b, MPI_Comm comm)
{
int k;
const int nent = gA->local.nent;
double mn = HUGE_VAL, gmn;
for (k = 0; k < nent; ++k) {
const double e = b[k];
if (e < mn) mn = e;
}
MPI_Allreduce (&mn, &gmn, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
for (k = 0; k < nent; ++k) {
b[k] -= gmn;
}
/*
{
int rank;
MPI_Comm_rank (comm, &rank);
fprintf (stderr, "%d: shifting by %g (%g)\n", rank, gmn, mn);
}
*/
}
void
stripe_auction_vars_l2g (const struct gcsr_t *gA,
int *lmatch, double *lprice,
int *match, double *price,
MPI_Comm comm)
{
/* written for debugging, not efficiency */
int rank = gA->rowmap.p, csize = gA->rowmap.np;
int am_root;
int p;
/*
const int lnr = gA->local.nr;
const int gnr = gA->nr;
*/
const int lnc = gA->local.nc;
const int gnc = gA->nc;
const int *l2g = gA->colmap.l2g;
const int lindstart = gA->rowmap.lindstart;
const int lnind = gA->rowmap.lnind;
int *matchbuf_out, *matchbuf;
double *pricebuf_out, *pricebuf;
int j, lj;
am_root = (0 == rank);
matchbuf_out = malloc (gnc * sizeof(int));
pricebuf_out = malloc (gnc * sizeof(double));
#if !defined(NDEBUG)
for (j = 0; j < gnc; ++j) matchbuf_out[j] = -(rank + 5);
for (j = 0; j < gnc; ++j) pricebuf_out[j] = -(rank + 5);
if (am_root) {
memset (match, -1, gnc * sizeof(int));
for (j = 0; j < gnc; ++j) price[j] = -1.0;
}
#else
memset (matchbuf_out, -1, gnc * sizeof(int));
memset (pricebuf_out, 0, gnc * sizeof(double));
if (am_root) {
memset (match, -1, gnc * sizeof(int));
memset (price, 0, gnc * sizeof(double));
}
#endif
for (lj = 0; lj < lnc; ++lj) {
j = (l2g? l2g[lj] : lj);
if (lmatch[lj] >= 0) {
assert(lmatch[lj] < lindstart + lnind);
matchbuf_out[j] = lmatch[lj] + lindstart;
}
pricebuf_out[j] = lprice[lj];
}
#if !defined(NDEBUG)
if (0 == csize) {
for (j = 0; j < gnc; ++j) {
assert(matchbuf_out[j] == lmatch[j]);
assert(pricebuf_out[j] == lprice[j]);
}
}
#endif
if (am_root) {
matchbuf = malloc (csize * gnc * sizeof(int));
pricebuf = malloc (csize * gnc * sizeof(double));
#if !defined(NDEBUG)
for (j = 0; j < gnc; ++j) matchbuf[j] = -(rank + 5)*10;
for (j = 0; j < gnc; ++j) pricebuf[j] = -(rank + 5)*10;
#endif
}
MPI_Gather (matchbuf_out, gnc, MPI_INT, matchbuf, gnc, MPI_INT, 0, comm);
MPI_Gather (pricebuf_out, gnc, MPI_DOUBLE, pricebuf, gnc, MPI_DOUBLE, 0, comm);
if (am_root) { /* check and collapse */
for (j = 0; j < gnc; ++j) {
int set_by = -1;
int price_set_by = -1;
double price_so_far = -1.0;
for (p = 0; p < csize; ++p) {
int i;
double pr;
i = matchbuf[gnc*p + j];
pr = pricebuf[gnc*p + j];
/*
if (1264 == j) {
fprintf (stderr, "col %d, proc %d: pb[%d*%d + %d = %d] = %g\n",
j, p, gnc, p, j, gnc*p + j, pr);
fprintf (stderr, "col %d, proc %d: mb[%d*%d + %d = %d] = %d\n",
j, p, gnc, p, j, gnc*p + j, i);
}
*/
if (i >= 0) {
if (set_by >= 0) {
fprintf (stderr,
"column collision on col %d, %d setting, first set by %d\n",
j, p, set_by);
}
else {
set_by = p;
match[j] = i;
}
}
if (pr >= 0.0) {
if (price_set_by >= 0 && pr != price_so_far) {
fprintf (stderr,
"prices wrong on col %d, %g on %d, first %g by %d\n",
j, pr, p, price_so_far, price_set_by);
}
else {
price_so_far = pr;
price_set_by = p;
price[j] = pr;
}
}
}
}
}
MPI_Barrier(comm);
}
int
stripe_auction_simple (const struct gcsr_t *gA, const double *lb,
int *lmatch, double *lprice,
const double mu_min,
struct comm_state_t *comm_state,
MPI_Comm comm)
{
const struct spcsr_t local = gA->local;
const int lnr = gA->local.nr;
const int lnc = gA->local.nc;
const int gnc = gA->nc;
const int gnr = gA->nr;
const int lnent = local.nent;
int max_lnc = 0;
int rank = gA->rowmap.p, csize = gA->rowmap.np;
int g_n_unmatched = gnr;
int n_unmatched = 0;
int *unmatched;
int n_changes;
int *changed_col_list, *col_changed_flag;
const int * restrict l2g = gA->colmap.l2g;
int * restrict g2l = NULL;
const int lindstart = gA->rowmap.lindstart;
const int lnind = gA->rowmap.lnind;
int j, k;
int pass = 0;
/*
fprintf (stderr, "%d: gA: (%d, %d)\n", rank, gA->nr, gA->nc);
fprintf (stderr, "%d: gA->local: (%d, %d; %d)\n", rank,
gA->local.nr, gA->local.nc, gA->local.nent);
*/
memset (lmatch, -1, lnc * sizeof(int));
memset (lprice, 0, lnc * sizeof(double));
init_comm_state (comm_state, gA, comm);
g2l = comm_state_g2l(comm_state);
unmatched = malloc (lnr * sizeof(int));
for (n_unmatched = 0; n_unmatched < lnr; ++n_unmatched) {
unmatched[n_unmatched] = n_unmatched;
}
g_n_unmatched = comm_g_n_unmatched (comm_state, n_unmatched, 0, comm);
changed_col_list = malloc (lnc * sizeof(int));
#if !defined(NDEBUG)
memset (changed_col_list, -1, lnc * sizeof(int));
#endif
col_changed_flag = malloc (lnc * sizeof(int));
memset (col_changed_flag, -1, lnc * sizeof(int));
/*
fprintf (stderr, "%d: entering, g_n_unmatched = %d, n_unmatched = %d\n",
rank, g_n_unmatched, n_unmatched);
*/
while (g_n_unmatched > 0) {
/*
fprintf (stderr, "%d: pass %d, g_n_unmatched = %d, n_unmatched = %d\n",
rank, pass, g_n_unmatched, n_unmatched);
*/
n_changes = 0;
while (n_unmatched > 0) {
n_unmatched = auction_pass_cl (local, lb, lmatch, lprice,
n_unmatched, unmatched,
&n_changes,
changed_col_list, col_changed_flag,
mu_min, -HUGE_VAL);
}
/*
{
double primal, dual;
auction_eval_primal_mdual (local, lb, lmatch, lprice, &primal, &dual);
fprintf (stderr,
"%d: local primal: %g local dual: %g\n", rank, primal, dual);
}
*/
/*fprintf (stderr, "%d: giving %d changes\n", rank, n_changes);*/
++comm_ngives;
comm_give_changes (comm_state, gA, n_unmatched,
n_changes, changed_col_list, col_changed_flag,
lmatch, lprice, comm);
{
int k;
for (k = 0; k < n_changes; ++k) {
int j, gj;
j = changed_col_list[k];
#if !defined(NDEBUG)
assert(col_changed_flag[j] == k);
changed_col_list[k] = -1;
#endif
col_changed_flag[j] = -1;
}
n_changes = 0;
}
{
int n_changes_in;
const struct cbuf_entry *clist;
while (comm_get_changes (comm_state, comm, &n_changes_in, &clist)) {
/*fprintf (stderr, "%d: merging %d changes\n", rank, n_changes_in);*/
merge_cbuf_into_vars (n_changes_in, clist,
g2l,
gA, lmatch, lprice,
&n_unmatched, unmatched);
}
}
g_n_unmatched = comm_g_n_unmatched (comm_state, n_unmatched, 0, comm);
++pass;
}
infeas:
/* free stuff*/
free_comm_state (comm_state, comm);
if (col_changed_flag) free (col_changed_flag);
if (changed_col_list) free (changed_col_list);
if (unmatched) free (unmatched);
return g_n_unmatched;
}
int
stripe_auction_scaling (const struct gcsr_t *gA, const double *lb_in,
int *lmatch_in, double *lprice_in,
const double mu_final, const int relgap,
struct comm_state_t *comm_state,
MPI_Comm comm)
{
const struct spcsr_t local = gA->local;
const int lnr = gA->local.nr;
const int lnc = gA->local.nc;
const int gnc = gA->nc;
const int gnr = gA->nr;
const int lnent = local.nent;
int max_lnc = 0;
int rank = gA->rowmap.p, csize = gA->rowmap.np;
int g_n_unmatched = gnr;
int n_unmatched = 0;
int *unmatched;
double tmp, tmpv1[3], tmpv2[3];
double mu, pmax, feasbnd;
double brange, bmn, bmx, gap;
int n_changes;
int *changed_col_list, *col_changed_flag;
const double * restrict lb = lb_in;
int * restrict lmatch = lmatch_in;
double * restrict lprice = lprice_in;
const double theta = 1.0/8.0;
const int * restrict l2g = gA->colmap.l2g;
int * restrict g2l = NULL;
const int lindstart = gA->rowmap.lindstart;
const int lnind = gA->rowmap.lnind;
int j, k;
int pass = 0;
/*
fprintf (stderr, "%d: gA: (%d, %d)\n", rank, gA->nr, gA->nc);
fprintf (stderr, "%d: gA->local: (%d, %d; %d)\n", rank,
gA->local.nr, gA->local.nc, gA->local.nent);
*/
memset (lmatch, -1, lnc * sizeof(int));
memset (lprice, 0, lnc * sizeof(double));
init_comm_state (comm_state, gA, comm);
g2l = comm_state_g2l (comm_state);
unmatched = malloc (lnr * sizeof(int));
for (n_unmatched = 0; n_unmatched < lnr; ++n_unmatched) {
unmatched[n_unmatched] = n_unmatched;
}
changed_col_list = malloc (lnc * sizeof(int));
#if !defined(NDEBUG)
memset (changed_col_list, -1, lnc * sizeof(int));
#endif
col_changed_flag = malloc (lnc * sizeof(int));
memset (col_changed_flag, -1, lnc * sizeof(int));
/* initial greedy matching while scanning largest benefit */
n_changes = 0;
n_unmatched = auction_greedy_gather0_cl (local, lb, lmatch, lprice,
unmatched,
&n_changes,
changed_col_list,
col_changed_flag,
&bmn, &bmx, &gap, &pmax);
tmpv1[0] = -bmn; tmpv1[1] = bmx; tmpv1[2] = pmax;
MPI_Allreduce (&tmpv1[0], &tmpv2[0], 3, MPI_DOUBLE, MPI_MAX, comm);
bmn = -tmpv2[0]; bmx = tmpv2[1]; pmax = tmpv2[2];
brange = bmx - bmn;
MPI_Allreduce (&gap, &tmp, 1, MPI_DOUBLE, MPI_SUM, comm);
gap = tmp;
g_n_unmatched = comm_g_n_unmatched (comm_state, n_unmatched, n_changes, comm);
mu = theta * gap / gnc;
/*
fprintf (stderr, "%d: pass: %d gap: %g mu: %g\n", rank, pass, gap, mu);
*/
if (!g_n_unmatched) goto freemem;
if (gap == -HUGE_VAL) { n_unmatched = -1; goto freemem; }
if (mu < 8.0*mu_final) mu = mu_final;
feasbnd = -(2 * gnc - 1) * brange - gnc * mu_final;
for (;;) {
double new_mu;
/* the matchall loop */
while (g_n_unmatched > 0) {
/*
fprintf (stderr, "%d: pass %d, g_n_unmatched = %d, n_unmatched = %d\n",
rank, pass, g_n_unmatched, n_unmatched);
*/
while (n_unmatched > 0) {
/*
fprintf (stderr, "%d: n_unmatched before pass_cl: %d\n", rank,
n_unmatched);
*/
n_unmatched = auction_pass_cl (local, lb, lmatch, lprice,
n_unmatched, unmatched,
&n_changes,
changed_col_list, col_changed_flag,
mu, -HUGE_VAL);
/*
fprintf (stderr, "%d: n_unmatched after pass_cl: %d\n", rank,
n_unmatched);
*/
}
#if 0
{
double primal, dual;
auction_eval_primal_mdual (local, lb, lmatch, lprice, &primal, &dual);
/*
fprintf (stderr,
"%d: local primal: %g local dual: %g gap %g\n",
rank, primal, dual, dual - primal);
*/
}
#endif
/*fprintf (stderr, "%d, %d: giving %d changes\n", rank, pass, n_changes);*/
++comm_ngives;
comm_give_changes (comm_state, gA, n_unmatched,
n_changes, changed_col_list, col_changed_flag,
lmatch, lprice, comm);
{ /* clear changes */
int k;
for (k = 0; k < n_changes; ++k) {
int j, gj;
j = changed_col_list[k];
#if !defined(NDEBUG)
assert(col_changed_flag[j] == k);
changed_col_list[k] = -1;
#endif
col_changed_flag[j] = -1;
}
n_changes = 0;
}
{
int n_changes_in;
const struct cbuf_entry *clist;
while (comm_get_changes (comm_state, comm, &n_changes_in, &clist)) {
/*
fprintf (stderr, "%d, %d: merging %d changes\n", rank, pass, n_changes_in);
*/
merge_cbuf_into_vars (n_changes_in, clist,
g2l,
gA, lmatch, lprice,
&n_unmatched, unmatched);
}
}
g_n_unmatched = comm_g_n_unmatched (comm_state,
n_unmatched, n_changes, comm);
/*fprintf (stderr, "%d: nun %d gnun %d\n", rank, n_unmatched, g_n_unmatched);*/
} /* the matchall loop */
/*
fprintf (stderr, "%d: pass: %d g_n_un: %d (%d) mu %g muf %g\n",
rank, pass, g_n_unmatched, n_unmatched, mu, mu_final);
*/
if (n_unmatched || mu <= mu_final) {
/*
fprintf (stderr, "%d: aauugghh %d, or %g <= %g\n",
rank, n_unmatched, mu, mu_final);
*/
break;
}
{
int j;
double pg[2] = {0.0, 0.0};
double all_pg[2] = {0.0, 0.0};
pmax = -HUGE_VAL;
pg[0] = auction_eval_gap (local, lb, lmatch, lprice, &pg[1]);
/*fprintf (stderr, "%d: pass: %d gap: %g mu: %g\n", rank, pass, gap, mu);*/
MPI_Allreduce (&pg[0], &all_pg[0], 2, MPI_DOUBLE, MPI_SUM, comm);
gap = all_pg[1];
if (!relgap) {
#if 1
if (gap < gnc * mu_final) break; /* before clearing match */
#else
if (gap < mu_final) break; /* before clearing match */
#endif
}
else {
if (all_pg[1] / all_pg[0] < mu_final) break;
}
new_mu = theta * (gap) / gnc;
if (new_mu < 8.0*mu_final) new_mu = mu_final;
if (new_mu > theta * mu) new_mu = theta * mu;
/*fprintf (stderr, "%d: pass: %d gap: %g mu: %g\n", rank, pass, gap, mu);*/
for (j = 0; j < lnc; ++j) {
const int i = lmatch[j];
const double p = lprice[j];
if (p == HUGE_VAL) continue;
if (p > pmax) pmax = p;
if (i < 0) continue;
lmatch[j] = -1;
unmatched[n_unmatched] = i;
++n_unmatched;
}
/* fprintf (stderr, "%d: %g ------- %g\n", rank, mu, new_mu); */
mu = new_mu;
g_n_unmatched = comm_g_n_unmatched (comm_state, n_unmatched, n_changes, comm);
MPI_Allreduce (&pmax, &tmp, 1, MPI_DOUBLE, MPI_MAX, comm);
pmax = tmp;
}
++pass;
}
freemem:
/* free stuff*/
free_comm_state (comm_state, comm);
if (col_changed_flag) free (col_changed_flag);
if (changed_col_list) free (changed_col_list);
if (unmatched) free (unmatched);
return g_n_unmatched;
}