#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <errno.h>
#include <assert.h>
#include "smut.h"
#if !defined(SMUT_USE_RESTRICT)
#define restrict
#endif
void
spcsr_init_clear (struct spcsr_t *A)
{
if (!A) return;
A->nr = A->nc = A->nent = 0;
A->rowoff = A->colind = A->coldeg = NULL;
A->entry = NULL;
}
void
spcsr_free (struct spcsr_t *A)
{
if (!A) return;
if (A->coldeg) free (A->coldeg);
if (A->entry) free (A->entry);
if (A->colind) free (A->colind);
if (A->rowoff) free (A->rowoff);
spcsr_init_clear (A);
}
int
spcsr_alloc_ (int nr, int nc, int nent,
int **rowoff, int **colind, double **entry,
int **coldeg)
{
int *ro = NULL, *ci = NULL;
double *e = NULL;
int *cd = NULL;
ro = malloc ((nr+1) * sizeof(int));
ci = malloc (nent * sizeof(int));
e = malloc (nent * sizeof(double));
if (coldeg) cd = malloc (nc * sizeof(int));
if (ro && ci && e && (cd || coldeg == NULL)) {
*rowoff = ro;
*colind = ci;
*entry = e;
if (coldeg) *coldeg = cd;
return 0;
}
if (cd) free(cd);
if (e) free(e);
if (ci) free(ci);
if (ro) free(ro);
return 1;
}
int
spcsr_load (const int nrows, const int ncols,
const int *rowoff, const int *colind, const double *entry,
const int *coldeg,
struct spcsr_t* A)
{
int nent;
if (!A) return 1;
if (!rowoff || !colind || !entry) goto errclear;
nent = rowoff[nrows];
if (nent < 0) goto errclear;
if (spcsr_alloc_ (nrows, ncols, nent, &A->rowoff, &A->colind, &A->entry,
(coldeg? &A->coldeg : NULL)) )
goto errclear;
A->nr = nrows;
A->nc = ncols;
A->nent = nent;
memcpy(A->rowoff, rowoff, (nrows+1)*sizeof(int));
memcpy(A->colind, colind, nent*sizeof(int));
memcpy(A->entry, entry, nent*sizeof(double));
if (coldeg)
memcpy(A->coldeg, coldeg, nent*sizeof(double));
return 0;
errclear:
spcsr_init_clear (A);
return 1;
}
int
spcsr_take (const int nrows, const int ncols,
int **rowoff, int **colind, double **entry,
int **coldeg,
struct spcsr_t *A)
{
int nent;
if (!A) return 1;
if (!rowoff || !colind || !entry) goto errclear;
nent = (*rowoff)[nrows];
if (nent < 0) goto errclear;
A->nr = nrows;
A->nc = ncols;
A->nent = nent;
A->rowoff = *rowoff; *rowoff = NULL;
A->colind = *colind; *colind = NULL;
A->entry = *entry; *entry = NULL;
if (coldeg) {
A->coldeg = *coldeg; *coldeg = NULL;
} else {
A->coldeg = NULL;
}
return 0;
errclear:
spcsr_init_clear (A);
return 1;
}
void
spcsr_get_sizes (const struct spcsr_t *A,
int *nrows, int *ncols, int *nents)
{
if (A && nrows && ncols) {
*nrows = A->nr;
*ncols = A->nc;
*nents = A->nent;
}
}
int
spcsr_copy (const struct spcsr_t *A, struct spcsr_t *Aout)
{
if (!A || !Aout) return 1;
*Aout = *A;
if (spcsr_alloc_ (A->nr, A->nc, A->nent,
&Aout->rowoff, &Aout->colind, &Aout->entry,
(A->coldeg? &Aout->coldeg : NULL)))
goto errclear;
memcpy(Aout->rowoff, A->rowoff, (A->nr+1)*sizeof(int));
memcpy(Aout->colind, A->colind, A->nent*sizeof(int));
memcpy(Aout->entry, A->entry, A->nent*sizeof(double));
if (A->coldeg)
memcpy(Aout->coldeg, A->coldeg, A->nent*sizeof(double));
return 0;
errclear:
spcsr_init_clear (Aout);
return 1;
}
void
spcsr_copyout (const struct spcsr_t *A,
int *rowoff, int *colind, double *entry,
int *coldeg)
{
if (!A) return;
if (rowoff) memcpy (rowoff, A->rowoff, (A->nr+1) * sizeof(int));
if (colind) memcpy (colind, A->colind, A->nent * sizeof(int));
if (entry) memcpy (entry, A->entry, A->nent * sizeof(double));
if (coldeg && A->coldeg)
memcpy (coldeg, A->coldeg, A->nc * sizeof(int));
}
void
spcsr_giveout (struct spcsr_t *A,
int **rowoff, int **colind, double **entry,
int **coldeg)
{
if (rowoff)
*rowoff = A->rowoff;
else
if (A->rowoff) free (A->rowoff);
if (colind)
*colind = A->colind;
else
if (A->colind) free (A->colind);
if (entry)
*entry = A->entry;
else
if (A->entry) free (A->entry);
if (coldeg)
*coldeg = A->coldeg;
else
if (A->coldeg) free (A->coldeg);
spcsr_init_clear (A);
}
static void
bswap32 (void *buf)
{
unsigned int *i = (unsigned int *) buf; /* XXX: Assume aligned */
unsigned int tmp;
if (!i) return;
tmp = *i;
*i = ((tmp & 0xff) << 24)
| ((tmp & 0xff00) << 8)
| ((tmp & 0xff0000) >> 8)
| ((tmp & 0xff000000) >> 24);
}
static void
bswap64 (void *buf)
{
unsigned int *i = (unsigned int *) buf; /* XXX: Assume aligned */
unsigned int tmp1, tmp2;
if (!i) return;
tmp1 = i[0];
tmp2 = i[1];
tmp1 = ((tmp1 & 0xff) << 24)
| ((tmp1 & 0xff00) << 8)
| ((tmp1 & 0xff0000) >> 8)
| ((tmp1 & 0xff000000) >> 24);
tmp2 = ((tmp2 & 0xff) << 24)
| ((tmp2 & 0xff00) << 8)
| ((tmp2 & 0xff0000) >> 8)
| ((tmp2 & 0xff000000) >> 24);
i[1] = tmp1;
i[0] = tmp2;
}
int
spcsr_load_binfile (FILE *f, struct spcsr_t *A)
{
int needswap;
int i, tag, nr, nc, nent;
int *off = NULL, *ind = NULL;
double *ent = NULL;
int errcode = 0;
if (!f || !A) return 0;
if (1 != fread (&tag, sizeof(int), 1, f)) return -1;
needswap = (tag != 1);
if (1 != fread (&nr, sizeof(int), 1, f)) return -2;
if (needswap) bswap32 (&nr);
if (1 != fread (&nc, sizeof(int), 1, f)) return -3;
if (needswap) bswap32 (&nc);
off = malloc ( (nr+1) * sizeof(int) );
if (!off) { errcode = ENOMEM; goto errfree; }
if (nr+1 != fread (off, sizeof(int), nr+1, f)) {
errcode = -5;
goto errfree;
}
if (needswap)
for (i = 0; i <= nr; ++i) bswap32(&off[i]);
nent = off[nr];
if (nent < 0) { errcode = EINVAL; goto errfree; }
ind = malloc ( nent * sizeof(int) );
ent = malloc ( nent * sizeof(double) );
if (!ind || !ent) { errcode = ENOMEM; goto errfree; }
if (nent != fread (ind, sizeof(int), nent, f)) {
errcode = -6;
goto errfree;
}
if (needswap)
for (i = 0; i < nent; ++i) bswap32(&ind[i]);
if (nent != fread (ent, sizeof(double), nent, f)) {
errcode = -7;
goto errfree;
}
if (needswap)
for (i = 0; i < nent; ++i) bswap64(&ent[i]);
spcsr_take (nr, nc, &off, &ind, &ent, NULL, A);
return 0;
errfree:
if (off) free (off);
if (ind) free (ind);
if (ent) free (ent);
return errcode;
}
int
smut_load_bin_vectors (FILE* f, int *nrp, int *ncp, double **datap)
{
int needswap;
int i, tag, nr, nc, nent;
double *data = NULL;
int errcode = 0;
if (!f || !nrp || !ncp || !data) return 0;
if (1 != fread (&tag, sizeof(int), 1, f)) return -1;
needswap = (tag != 1);
if (1 != fread (&nr, sizeof(int), 1, f)) return -2;
if (needswap) bswap32 (&nr);
if (1 != fread (&nc, sizeof(int), 1, f)) return -3;
if (needswap) bswap32 (&nc);
data = malloc (nr * nc * sizeof(double));
if (!data) { errcode = ENOMEM; goto errfree; }
if (nr * nc != fread (data, sizeof(double), nr*nc, f)) {
errcode = -5;
goto errfree;
}
if (needswap)
for (i = 0; i < nr*nc; ++i) bswap64(&data[i]);
*nrp = nr;
*ncp = nc;
*datap = data;
return 0;
errfree:
if (data) free (data);
return errcode;
}
int
smut_load_bin_vectors_bf (FILE* f, int *nrp, int *ncp,
smut_bigfloat_t **datap)
{
int needswap;
int i, tag, nr, nc, nent;
smut_bigfloat_t *data = NULL;
int errcode = 0;
if (!f || !nrp || !ncp || !data) return 0;
if (1 != fread (&tag, sizeof(int), 1, f)) return -1;
needswap = (tag != 1);
if (needswap) {
errcode = EINVAL;
goto errfree;
}
if (1 != fread (&nr, sizeof(int), 1, f)) return -2;
if (needswap) bswap32 (&nr);
if (1 != fread (&nc, sizeof(int), 1, f)) return -3;
if (needswap) bswap32 (&nc);
data = malloc (nr * nc * sizeof(smut_bigfloat_t));
if (!data) { errcode = ENOMEM; goto errfree; }
if (nr * nc != fread (data, sizeof(smut_bigfloat_t), nr*nc, f)) {
errcode = -5;
goto errfree;
}
*nrp = nr;
*ncp = nc;
*datap = data;
return 0;
errfree:
if (data) free (data);
return errcode;
}
int
spcsr_write_binfile (const struct spcsr_t *A, FILE *f)
{
int tag = 1;
if (!A || !f) return 0;
if (1 != fwrite(&tag, sizeof(int), 1, f)) return -1;
if (1 != fwrite(&A->nr, sizeof(int), 1, f)) return -2;
if (1 != fwrite(&A->nc, sizeof(int), 1, f)) return -3;
if (A->nr+1 != fwrite(A->rowoff, sizeof(int), A->nr+1, f))
return -4;
if (A->nent != fwrite(A->colind, sizeof(int), A->nent, f))
return -5;
if (A->nent != fwrite(A->entry, sizeof(double), A->nent, f))
return -6;
return 0;
}
int
smut_write_bin_vectors (const int nr, const int nc,
const double *data, FILE* f)
{
int tag = 1;
if (!nr || !nc || !data || !f) return 0;
if (1 != fwrite (&tag, sizeof(int), 1, f)) return -1;
if (1 != fwrite (&nr, sizeof(int), 1, f)) return -2;
if (1 != fwrite (&nc, sizeof(int), 1, f)) return -3;
if (nr*nc != fwrite (data, sizeof(double), nr*nc, f)) return -4;
return 0;
}
int
smut_write_bin_vectors_bf (const int nr, const int nc,
const smut_bigfloat_t *data, FILE *f)
{
int tag = 1;
if (!nr || !nc || !data || !f) return 0;
if (1 != fwrite (&tag, sizeof(int), 1, f)) return -1;
if (1 != fwrite (&nr, sizeof(int), 1, f)) return -2;
if (1 != fwrite (&nc, sizeof(int), 1, f)) return -3;
if (nr*nc != fwrite (data, sizeof(smut_bigfloat_t), nr*nc, f)) return -4;
return 0;
}
void
spcsr_coldeg_data (const int nr, const int nc, const int nent,
const int *rowoff, const int *colind,
const double *entry,
int *coldeg)
{
int k;
memset (coldeg, 0, nc * sizeof(int));
for (k = 0; k < nent; ++k) {
const int j = colind[k];
assert (j >= 0);
assert (j < nc);
++coldeg[j];
}
}
int
spcsr_update_coldeg (struct spcsr_t *A)
{
int *cd;
if (!A) return 0;
cd = malloc (A->nc * sizeof(int));
if (!cd) return 1;
spcsr_coldeg_data (A->nr, A->nc, A->nent, A->rowoff, A->colind, A->entry,
cd);
A->coldeg = cd;
return 0;
}
void
spcsr_transpose_data (const int nr, const int nc,
const int *rowoff_in, const int *colind_in,
const double *entry_in, const int *coldeg_in,
int *coloff_in, int *rowind_in, double *rowent_in)
{
int i, j, k, cumsum;
int nent;
const int restrict *rowoff = rowoff_in;
const int restrict *colind = colind_in;
const double restrict *entry = entry_in;
const int restrict *coldeg = coldeg_in;
int restrict *coloff = coloff_in;
int restrict *rowind = rowind_in;
double restrict *rowent = rowent_in;
if (!rowoff || !colind || !coloff || !rowind || !rowent) return;
nent = rowoff[nr];
if (!coldeg) {
cumsum = 0;
memset (coloff, 0, (nc+1) * sizeof(int));
for (k = 0; k < nent; ++k) {
++coloff[1+colind[k]];
}
for (j = 0; j < nc; ++j) {
int tmp = coloff[j+1];
coloff[j+1] = cumsum;
cumsum += tmp;
}
} else {
cumsum = 0;
coloff[0] = 0;
for (j = 0; j < nc; ++j) {
coloff[j+1] = cumsum;
cumsum += coldeg[j];
}
}
assert(nent == cumsum);
for (i = 0; i < nr; ++i) {
for (k = rowoff[i]; k < rowoff[i+1]; ++k) {
const int j = colind[k];
int kT = coloff[j+1]++;
rowind[kT] = i;
if (rowent && entry) {
const double e = entry[k];
rowent[kT] = e;
}
}
}
assert(nent == coloff[nc]);
}
int
spcsr_transpose (const struct spcsr_t *A, struct spcsr_t *A_T)
{
int nent;
int do_coldeg = 0;
int *rowoff = NULL, *colind = NULL, *coldeg = NULL;
double *entry = NULL;
if (!A || !A_T) return 0;
spcsr_init_clear (A_T);
A_T->nr = A->nc;
A_T->nc = A->nr;
nent = A_T->nent = A->nent;
if (A->coldeg) do_coldeg = 1;
rowoff = malloc ( (1 + A_T->nr) * sizeof(int) );
colind = malloc ( nent * sizeof(int) );
entry = malloc ( nent * sizeof(double) );
if (do_coldeg) coldeg = malloc (A_T->nc * sizeof(int));
if (!rowoff || !colind || !entry || (do_coldeg && !coldeg))
goto errclear;
spcsr_transpose_data (A->nr, A->nc, A->rowoff, A->colind, A->entry,
A->coldeg,
rowoff, colind, entry);
if (do_coldeg) {
const int nr = A->nr;
int i;
for (i = 0; i < nr; ++i)
coldeg[i] = A->rowoff[i+1] - A->rowoff[i];
}
A_T->rowoff = rowoff;
A_T->colind = colind;
A_T->entry = entry;
A_T->coldeg = coldeg;
return 0;
errclear:
if (coldeg) free(coldeg);
if (entry) free(entry);
if (colind) free(colind);
if (rowoff) free(rowoff);
return 1;
}
#if !defined(HAVE_LOGB)
#if defined(HAVE_LOG2)
#define logb(x) floor(log2(x))
#else
static double logb (double x) { int exp; frexp(x, &exp); return exp; }
#endif
#endif
#if !defined(HAVE_SCALB)
#define scalb(x,e) ldexp(x,(int)e);
#endif
void
spcsr_lascale (struct spcsr_t *A, double *R, double *C)
{
int i, j, k;
const int nr = A->nr;
const int nc = A->nc;
const int nent = A->nent;
int * restrict rowoff = A->rowoff;
int * restrict colind = A->colind;
double * restrict entry = A->entry;
for (j = 0; j < nc; ++j) C[j] = -HUGE_VAL;
for (i = 0; i < nr; ++i) {
const int rowend = rowoff[i+1];
double mx = -HUGE_VAL;
for (k = rowoff[i]; k < rowend; ++k) {
double fe;
fe = fabs(entry[k]);
if (fe > mx) mx = fe;
}
if (mx == -HUGE_VAL)
mx = 0.0;
else if (mx != 0.0)
mx = logb(mx);
R[i] = scalb(1.0, mx);
for (k = rowoff[i]; k < rowend; ++k) {
const int j = colind[k];
double se;
/* se = scalb(entry[k], -mx); */
se = entry[k] / R[i];
entry[k] = se;
se = fabs(se);
if (se > C[j]) C[j] = se;
}
}
for (j = 0; j < nc; ++j) {
double tmp = C[j];
if (tmp == HUGE_VAL) C[j] = 0.0;
else if (tmp != 0.0) tmp = logb(tmp);
C[j] = scalb(1.0, tmp);
}
for (k = 0; k < nent; ++k) {
const int j = colind[k];
double se;
se = entry[k] / C[j];
entry[k] = se;
}
}
void
spcsr_lascaling (struct spcsr_t *A, double *R_in, double *C_in)
{
int i, j, k;
const int nr = A->nr;
const int nc = A->nc;
int * restrict rowoff = A->rowoff;
int * restrict colind = A->colind;
double * restrict entry = A->entry;
double * restrict R = R_in;
double * restrict C = C_in;
for (j = 0; j < nc; ++j) C[j] = -HUGE_VAL;
for (i = 0; i < nr; ++i) {
const int rowend = rowoff[i+1];
double mx = -HUGE_VAL;
for (k = rowoff[i]; k < rowend; ++k) {
double fe;
fe = fabs(entry[k]);
if (fe > mx) mx = fe;
}
if (mx == -HUGE_VAL)
mx = 0.0;
else if (mx != 0.0)
mx = logb(mx);
R[i] = scalb(1.0, mx);
for (k = rowoff[i]; k < rowend; ++k) {
const int j = colind[k];
double se;
/*se = scalb(entry[k], -mx);*/
se = entry[k] / R[i];
se = fabs(se);
if (se > C[j]) C[j] = se;
}
}
for (j = 0; j < nc; ++j) {
double tmp = C[j];
if (tmp == HUGE_VAL) C[j] = 0.0;
else if (tmp != 0.0) tmp = logb(tmp);
C[j] = scalb(1.0,tmp);
}
}
void
spcsr_apply_scaling (struct spcsr_t *A, double *R_in, double *C_in)
{
int i;
const int nr = A->nr;
int * restrict rowoff = A->rowoff;
int * restrict colind = A->colind;
double * restrict entry = A->entry;
double * restrict R = R_in;
double * restrict C = C_in;
for (i = 0; i < nr; ++i) {
int k;
for (k = rowoff[i]; k < rowoff[i+1]; ++k) {
const int j = colind[k];
double e = entry[k];
e /= (R[i] * C[j]);
entry[k] = e;
}
}
}