int
read_input(FILE *infile, const char *filename,
           double **datap, int *ncolsp, int **cumsizesp, int *nrunsp)
{
    char b[2];
    double number;

    int newline;                /* last input was newline */
    int datai;                  /* the current element of (*datap) */
    int col;                    /* the column being read */
    int set = *nrunsp;          /* the current data set */
    int ncols = *ncolsp;        /* the number of columns */

    int datasize;
    int setsize;

    datai = (set == 0)
        ? 0
        : ncols * (*cumsizesp)[set-1];

    setsize = (set == 0)
        ? 0
        : ((set-1) / SET_INC + 1) * SET_INC;

    *cumsizesp = realloc(*cumsizesp, setsize * sizeof(int));

    datasize = (datai == 0)
        ? 0
        : ((datai - 1) / DATA_INC + 1) * DATA_INC;

    *datap = realloc(*datap, datasize * sizeof(double));

    /* remove leading whitespace */
    fscanf (infile, "%*[ \t\n]");

    if (feof(infile)) {
        warnprintf("%s: file is empty.", filename);
        return -1;
    }

    do {
        newline = 0;
        if (set == setsize) {
            setsize += SET_INC;
            *cumsizesp = realloc(*cumsizesp,  setsize * sizeof(int));
        }

        (*cumsizesp)[set] = (set == 0) ? 0 : (*cumsizesp)[set-1];       /* beginning of data set */

        while (newline == 0) {
        col = 0;        /* beginning of row */
            while (newline == 0) {
        if (fscanf(infile, "%lf", &number) != 1) {
                    char buffer[64];

                    fscanf(infile, "%60s", buffer);
                    errprintf(
                        "could not convert string `%s' to double, exiting...",
                        buffer);
                }

                if (datai == datasize) {
                    datasize += DATA_INC;
                    *datap = realloc(*datap,  datasize * sizeof(double));
                }
                (*datap)[datai++] = number;

#if DEBUG > 2
                fprintf(stderr, "set %d, row %d, column %d, x = %g\n",
                        set, (*cumsizesp)[set], col, number);
#endif
                col++;    /* new column */

                fscanf(infile, "%*[ \t]");
                newline = fscanf(infile, "%1[\n]", b);
            }

            if (!ncols)
                ncols = col;
            else if (col != ncols) {
                if ((*cumsizesp)[0] == 0)
                    errprintf ("reference point has dimension %d"
                               " while input has dimension %d",
                               ncols, col);
                else
                    errprintf("row %d has different length (%d)"
                              " than previous rows (%d), exiting...",
                              (*cumsizesp)[set], col, ncols);
            }

            (*cumsizesp)[set]++;

            fscanf (infile, "%*[ \t]");
        newline = fscanf (infile, "%1[\n]", b);
    }
#if DEBUG > 1
        fprintf (stderr, "Set %d, %d rows in total\n", set, (*cumsizesp)[set]);
#endif
    set++; /* new data set */

        fscanf (infile, "%*[ \t\n]");

    } while (!feof(infile));

    *ncolsp = ncols;
    *nrunsp = set;

    return 0;
}