#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <assert.h>
#if defined(HAVE_CALLGRIND)
#include <valgrind/callgrind.h>
#else
#define CALLGRIND_START_INSTRUMENTATION do { } while (0)
#define CALLGRIND_STOP_INSTRUMENTATION do { } while (0)
#define CALLGRIND_DUMP_STATS do { } while (0)
#define CALLGRIND_ZERO_STATS do { } while (0)
#define CALLGRIND_TOGGLE_COLLECT do { } while (0)
#endif
#include "smut.h"
#include "timer.h"
#include "auction.h"
#if !defined(HAVE_RESTRICT)
#define restrict
#endif
static int icntl[10];
void mc64ad_ (int* JOB,
int* N,
int* NE,
int* IP,
int* IRN,
double* A,
int* NUM,
int* CPERM,
int* LIW,
int* IW,
int* LDW,
double* DW,
int* ICNTL,
int* INFO);
void mc64wd_ (int* n, int* ne, int* ip, int* irn, double* a,
int* iperm, int* num, int* jperm, int* out,
int* pr, int* q, int* L, double* U, double* D);
int
main (int argc, char **argv)
{
FILE *f, *fout = NULL;
int err;
struct spcsr_t A;
double mu_min;
double *R, *C;
double *expint;
int *match, *invmatch;
double *price;
int errflg, i, j;
double p, d;
volatile double auction_primal = 0.0, mc64_primal = 0.0;
/* gcc-3.4 compiler bug? ap becomes zero after calling mc64
with optimization */
struct Timer timer;
assert(sizeof(int) == 4);
if (argc <= 1) {
printf("No files\n");
return -1;
}
f = fopen(argv[1], "rb");
if (!f) {
perror("Error opening file: ");
return -1;
}
if (argc > 2) {
fout = fopen (argv[2], "wb");
}
err = spcsr_load_binfile (f, &A);
if (err) {
printf ("Error reading: %d\n", err);
return -1;
}
fclose(f);
R = malloc (A.nr * sizeof(double));
C = malloc (A.nc * sizeof(double));
expint = malloc (A.nent * sizeof(double));
match = malloc (A.nc * sizeof(int));
invmatch = malloc (A.nr * sizeof(int));
memset (match, -1, A.nc * sizeof(int));
memset (invmatch, -1, A.nr * sizeof(int));
price = calloc (A.nc, sizeof(double));
printf ("nr = %d nc = %d nent = %d\n", A.nr, A.nc, A.nent);
if (getenv ("PRESCALE")) {
spcsr_lascaling (&A, R, C);
spcsr_apply_scaling (&A, R, C);
}
else
spcsr_lascale (&A, R, C);
if (getenv ("TOINT"))
auction_toexpint (&A, expint);
else
auction_toexp (&A, expint);
auction_shift (&A, expint);
{
int k;
double mx = -HUGE_VAL, mn = HUGE_VAL;
for (k = 0; k < A.nent; ++k) {
double e = expint[k];
if (e < mn) mn = e;
if (e > mx) mx = e;
}
printf ("expint max: %g , min %g\n", mx, mn);
}
cacheblast();
mu_min = 1.0 / A.nr;
/* be approximate */
/*mu_min = 1.0;*/
/*mu_min = 100.0 / A.nr;*/
#if 0
{
int n_unmatched;
double mu_min;
double p, d;
initialize_timer (&timer);
start_timer (&timer);
n_unmatched = auction_simple (A, expint, match, price, mu_min);
stop_timer (&timer);
auction_eval_primal_mdual (A, expint, match, price, &p, &d);
printf ("* Primal: %20g\n* Dual: %20g\n* Time: %20g\n", p, d,
timer_duration(timer));
a_printstats ();
}
#endif
cacheblast();
memset (match, -1, A.nc * sizeof(int));
memset (price, 0, A.nc * sizeof(double));
CALLGRIND_START_INSTRUMENTATION;
initialize_timer (&timer);
start_timer (&timer);
errflg = auction_scaling (A, expint, match, price, mu_min, 0);
stop_timer (&timer);
if (errflg)
fprintf (stderr, "Errflg: %d\n", errflg);
auction_eval_primal_mdual (A, expint, match, price, &auction_primal, &d);
printf ("Primal: %20g\nDual: %20g\nTime: %20g\nmu_min: %20g\n",
auction_primal, d, timer_duration(timer), mu_min);
a_printstats();
{
int i, j;
for (j = 0; j < A.nc; ++j) {
if (match[j] < 0) printf("column %d is unmatched\n", j);
if (match[j] >= A.nr) printf("column %d's match out of range (%d/%d)\n",
j, match[j], A.nr);
}
for (i = 0; i < A.nr; ++i) {
int k;
int found_match = 0;
for (k = A.rowoff[i]; k < A.rowoff[i+1]; ++k) {
const int j = A.colind[k];
if (i == match[j]) ++found_match;
}
if (!found_match) printf ("row %d is unmatched\n", i);
if (found_match > 1) printf ("row %d is multiply matched\n", i);
}
}
if (fout) {
int tag = 1;
fwrite (&tag, sizeof(int), 1, fout);
fwrite (&A.nc, sizeof(int), 1, fout);
fwrite (match, sizeof(int), A.nc, fout);
}
#if 0
cacheblast();
{
int n_unmatched;
double mu_min;
double p, d;
initialize_timer (&timer);
start_timer (&timer);
n_unmatched = auction_scaling_cache (A, expint, match, price, mu_min);
stop_timer (&timer);
auction_eval_primal_mdual (A, expint, match, price, &p, &d);
printf ("_ Primal: %20g\n_ Dual: %20g\n_ Time: %20g\n", p, d,
timer_duration(timer));
a_printstats ();
}
#endif
CALLGRIND_STOP_INSTRUMENTATION;
cacheblast();
/* mc64 */
if (getenv ("RUN_MC64")) {
{
int job = 4; /* sum of diag, already log scale */
int num_on_diag = 0;
int liw;
int *iw;
int ldw;
double *dw;
int info = 0;
int i, k;
for (i = 0; i <= A.nr; ++i)
++A.rowoff[i];
for (k = 0; k < A.nent; ++k)
++A.colind[k];
initialize_timer (&timer);
start_timer (&timer);
liw = 5 * A.nr;
iw = malloc (liw * sizeof(int));
ldw = 2 * A.nr + A.nent;
dw = malloc (ldw * sizeof(double));
memset (icntl, 0, 10 * sizeof(int));
icntl[0] = 6; icntl[1] = 6; icntl[2] = -1;
#if !defined(NDEBUG)
icntl[3] = 0;
#else
icntl[3] = -1;
#endif
mc64ad_ (&job, &A.nr, &A.nent, A.rowoff, A.colind, expint,
&num_on_diag, invmatch, &liw, iw, &ldw, dw, icntl,
&info);
free (iw);
free (dw);
stop_timer (&timer);
for (i = 0; i < A.nr; ++i) --invmatch[i];
for (i = 0; i <= A.nr; ++i) --A.rowoff[i];
for (k = 0; k < A.nent; ++k)
--A.colind[k];
}
for (i = 0; i < A.nr; ++i) {
if (invmatch[i] < 0) printf("mc64: row %d is unmatched\n", i);
}
memset (match, -1, A.nc * sizeof(int));
for (i = 0; i < A.nr; ++i) {
const int j = invmatch[i];
match[j] = i;
}
for (j = 0; j < A.nc; ++j) {
if (match[j] < 0) printf("mc64: col %d is unmatched\n", j);
}
auction_eval_primal_mdual (A, expint, invmatch, price, &mc64_primal, &d);
printf ("mc64 Primal: %20g\nmc64 Time: %20g\n", mc64_primal,
timer_duration(timer));
printf (" primal diff: %20g\n", auction_primal - mc64_primal);
}
free (price);
free (invmatch);
free (match);
free (expint);
free (R);
free (C);
spcsr_free (&A);
return 0;
}