struct SubSig *I_NewSubSig(struct SigSet *S, struct ClassSig *C)
{
    struct SubSig *Sp;
    int i;

    if (C->nsubclasses == 0)
	C->SubSig = (struct SubSig *)G_malloc(sizeof(struct SubSig));
    else
	C->SubSig = (struct SubSig *)G_realloc((char *)C->SubSig,
					       sizeof(struct SubSig) *
					       (C->nsubclasses + 1));

    Sp = &C->SubSig[C->nsubclasses++];
    Sp->used = 1;
    Sp->R = (double **)G_calloc(S->nbands, sizeof(double *));
    Sp->R[0] = (double *)G_calloc(S->nbands * S->nbands, sizeof(double));
    for (i = 1; i < S->nbands; i++)
	Sp->R[i] = Sp->R[i - 1] + S->nbands;
    Sp->Rinv = (double **)G_calloc(S->nbands, sizeof(double *));
    Sp->Rinv[0] = (double *)G_calloc(S->nbands * S->nbands, sizeof(double));
    for (i = 1; i < S->nbands; i++)
	Sp->Rinv[i] = Sp->Rinv[i - 1] + S->nbands;
    Sp->means = (double *)G_calloc(S->nbands, sizeof(double));
    Sp->N = 0;
    Sp->pi = 0;
    Sp->cnst = 0;
    return Sp;
}

/*!
   \brief Thin line

   For now, just eliminate points at regular interval

   \param gln line (geoline)
   \param factor

   \return pointer to geoline struct
   \return NULL on failure
 */
static geoline *thin_line(geoline * gln, float factor)
{
    geoline *newln;
    int i, nextp, targp;

    newln = (geoline *) G_malloc(sizeof(geoline));	/* G_fatal_error */
    if (!newln) {
	return (NULL);
    }

    targp = (int)(gln->npts / factor);

    if (targp < 2) {
	targp = 2;
    }

    newln->npts = targp;

    if (2 == (newln->dims = gln->dims)) {
	newln->p2 = (Point2 *) G_calloc(targp, sizeof(Point2));	/* G_fatal_error */
	if (!newln->p2) {
	    return (NULL);
	}

	for (i = 0; i < targp; i++) {
	    if (i == targp - 1) {
		nextp = gln->npts - 1;	/* avoid rounding error */
	    }
	    else {
		nextp = (int)((i * (gln->npts - 1)) / (targp - 1));
	    }

	    newln->p2[i][X] = gln->p2[nextp][X];
	    newln->p2[i][Y] = gln->p2[nextp][Y];
	}
    }
    else {
	newln->p3 = (Point3 *) G_calloc(targp, sizeof(Point3));	/* G_fatal_error */
	if (!newln->p3) {
	    return (NULL);
	}

	for (i = 0; i < targp; i++) {
	    if (i == targp - 1) {
		nextp = gln->npts - 1;	/* avoid rounding error */
	    }
	    else {
		nextp = (int)((i * (gln->npts - 1)) / (targp - 1));
	    }

	    newln->p3[i][X] = gln->p3[nextp][X];
	    newln->p3[i][Y] = gln->p3[nextp][Y];
	    newln->p3[i][Z] = gln->p3[nextp][Z];
	}
    }

    newln->next = NULL;

    return (newln);
}

int o_alloc_bufs(int ncols, int size)
{
    buffer[0] = (void *)G_calloc(ncols, size);
    buffer[1] = (void *)G_calloc(ncols, size);

    return 0;
}

/*!
 * \brief Allocate memory for a quadratic or not quadratic linear equation system
 *
 * The type of the linear equation system must be G_MATH_NORMAL_LES for
 * a regular quadratic matrix or G_MATH_SPARSE_LES for a sparse matrix
 *
 * <p>
 * In case of G_MATH_NORMAL_LES
 * 
 * A quadratic matrix of size rows*rows*sizeof(double) will allocated
 *
 * <p>
 * In case of G_MATH_SPARSE_LES
 *
 * a vector of size row will be allocated, ready to hold additional allocated sparse vectors.
 * each sparse vector may have a different size.
 *
 * Parameter parts defines which parts of the les should be allocated.
 * The number of columns and rows defines if the matrix is quadratic.
 *
 * \param rows int
 * \param cols int
 * \param type int
 * \param parts int -- 2 = A, x and b; 1 = A and x; 0 = A allocated
 * \return G_math_les *
 *
 * */
G_math_les *G_math_alloc_les_param(int rows, int cols, int type, int parts)
{
	G_math_les *les;

	if (type == G_MATH_SPARSE_LES)
		G_debug(
				2,
				"Allocate memory for a sparse linear equation system with %i rows\n",
				rows);
	else
		G_debug(
				2,
				"Allocate memory for a regular linear equation system with %i rows and %i cols\n",
				rows, cols);

	les = (G_math_les *) G_calloc(1, sizeof(G_math_les));
	les->x = NULL;
	les->b = NULL;

	if (parts > 0)
	{
		les->x = (double *)G_calloc(cols, sizeof(double));
	}

	if (parts > 1)
	{
		les->b = (double *)G_calloc(cols, sizeof(double));
	}

	les->A = NULL;
	les->data = NULL;
	les->Asp = NULL;
	les->rows = rows;
	les->cols = cols;
	les->symm = 0;
	les->bandwith = cols;
	if (rows == cols)
		les->quad = 1;
	else
		les->quad = 0;

	if (type == G_MATH_SPARSE_LES)
	{
		les->Asp = (G_math_spvector **) G_calloc(rows,
				sizeof(G_math_spvector *));
		les->type = G_MATH_SPARSE_LES;
	}
	else
	{
		les->A = G_alloc_matrix(rows, cols);
		/*save the start pointer of the matrix*/
		les->data = les->A[0];
		les->type = G_MATH_NORMAL_LES;
	}

	return les;
}

static void initialize_bins(void)
{
    int cat;

    G_message(_("Computing bins"));

    for (cat = 0; cat < num_cats; cat++) {
	struct basecat *bc = &basecats[cat];
	int slot;
	double next;
	int num_values = 0;
	int bin = 0;
	unsigned long accum = 0;
	int quant = 0;

	bc->bins = G_calloc(num_quants, sizeof(struct bin));
	bc->slot_bins = G_calloc(num_slots, sizeof(unsigned char));

	next = get_quantile(bc, quant);

	for (slot = 0; slot < num_slots; slot++) {
	    unsigned int count = bc->slots[slot];
	    unsigned long accum2 = accum + count;

	    if (accum2 > next) {
		struct bin *b = &bc->bins[bin];

		bc->slot_bins[slot] = ++bin;

		b->origin = accum;
		b->base = num_values;
		b->count = 0;
		b->min = min + slot_size * slot;
		b->max = min + slot_size * (slot + 1);

		while (accum2 > next)
		    next = get_quantile(bc, ++quant);

		num_values += count;
	    }

	    accum = accum2;
	}

	bc->num_values = num_values;
	bc->num_bins = bin;

	G_free(bc->slots);

	bc->values = G_calloc(num_values, sizeof(DCELL));
    }
}

/*!
 * \brief Allocate memory for a N_array_2d data structure.
 *
 * This function allocates memory for an array of type N_array_2d
 * and returns the pointer to the new allocated memory.
 * <br><br>
 * The data type of this array is set by "type" and must be
 * CELL_TYPE, FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
 * The offset sets the number of boundary cols and rows.
 * This option is useful to generate homogeneous Neumann boundary conditions around
 * an array or to establish overlapping boundaries. The array is initialized with 0 by default.
 * <br><br>
 * If the offset is greater then 0, negative indices are possible.
 * <br><br>
 *
 * The data structure of a array with 3 rows and cols and an offset of 1
 * will looks like this:
 * <br><br>
 *
 \verbatim
 0 0 0 0 0
 0 0 1 2 0
 0 3 4 5 0
 0 6 7 8 0
 0 0 0 0 0
 \endverbatim
 *
 * 0 is the boundary.
 * <br><br>
 * Internal a one dimensional array is allocated to save memory and to speed up the memory access.
 * To access the one dimensional array with a two dimensional index use the provided
 * get and put functions. The internal representation of the above data will look like this:
 *
 \verbatim
 0 0 0 0 0 0 0 1 2 0 0 3 4 5 0 0 6 7 8 0 0 0 0 0 0
 \endverbatim
 *
 * \param cols int
 * \param rows int
 * \param offset int
 * \param type int
 * \return N_array_2d *
 *
 * */
N_array_2d *N_alloc_array_2d(int cols, int rows, int offset, int type)
{
    N_array_2d *data = NULL;

    if (rows < 1 || cols < 1)
	G_fatal_error("N_alloc_array_2d: cols and rows should be > 0");

    if (type != CELL_TYPE && type != FCELL_TYPE && type != DCELL_TYPE)
	G_fatal_error
	    ("N_alloc_array_2d: Wrong data type, should be CELL_TYPE, FCELL_TYPE or DCELL_TYPE");

    data = (N_array_2d *) G_calloc(1, sizeof(N_array_2d));

    data->cols = cols;
    data->rows = rows;
    data->type = type;
    data->offset = offset;
    data->rows_intern = rows + 2 * offset;	/*offset position at booth sides */
    data->cols_intern = cols + 2 * offset;	/*offset position at booth sides */
    data->cell_array = NULL;
    data->fcell_array = NULL;
    data->dcell_array = NULL;

    if (data->type == CELL_TYPE) {
	data->cell_array =
	    (CELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
			      sizeof(CELL));
	G_debug(3,
		"N_alloc_array_2d: CELL array allocated rows_intern %i cols_intern %i offset %i",
		data->rows_intern, data->cols_intern, data->offset = offset);
    }
    else if (data->type == FCELL_TYPE) {
	data->fcell_array =
	    (FCELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
			       sizeof(FCELL));
	G_debug(3,
		"N_alloc_array_2d: FCELL array allocated rows_intern %i cols_intern %i offset %i",
		data->rows_intern, data->cols_intern, data->offset = offset);

    }
    else if (data->type == DCELL_TYPE) {
	data->dcell_array =
	    (DCELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
			       sizeof(DCELL));
	G_debug(3,
		"N_alloc_array_2d: DCELL array allocated rows_intern %i cols_intern %i offset %i",
		data->rows_intern, data->cols_intern, data->offset = offset);
    }

    return data;
}

struct Stats alloc_Stats(int n)
{
    double *err, *stm;
    struct Stats stat;

    stat.n_points = n;
    err = (double *)G_calloc(n, sizeof(double));
    stm = (double *)G_calloc(n, sizeof(double));

    stat.error = err;
    stat.estima = stm;

    return stat;
}

void filldir(int fe, int fd, int nl, struct band3 *bnd, struct metrics *m)
{
    int i, bufsz;
    CELL *dir;

    /* determine the flow direction in each cell.  On outer rows and columns
     * the flow direction is always directly out of the map */

    dir = G_calloc(bnd->ns, sizeof(CELL));
    bufsz = bnd->ns * sizeof(CELL);

    lseek(fe, 0, SEEK_SET);
    lseek(fd, 0, SEEK_SET);
    advance_band3(fe, bnd);
    for (i = 0; i < nl; i += 1) {
	advance_band3(fe, bnd);
	build_one_row(i, nl, bnd->ns, bnd, dir, m[i]);
	write(fd, dir, bufsz);
    }
    advance_band3(fe, bnd);
    build_one_row(i, nl, bnd->ns, bnd, dir, m[i]);
    write(fd, dir, bufsz);

    G_free(dir);

    return;
}

FLAG *FlagCreate(int nrows, int ncols)
{
    unsigned char *temp;
    FLAG *new_flag;
    register int i;

    new_flag = (FLAG *) G_malloc(sizeof(FLAG));

    new_flag->nrows = nrows;
    new_flag->ncols = ncols;
    new_flag->leng = (ncols + 7) / 8;
    new_flag->array =
	(unsigned char **)G_malloc(nrows * sizeof(unsigned char *));

    temp =
	(unsigned char *)G_calloc(nrows * new_flag->leng,
				  sizeof(unsigned char));

    for (i = 0; i < nrows; i++) {
	new_flag->array[i] = temp;
	temp += new_flag->leng;
    }

    return (new_flag);
}

vec_struct *G_vector_init(int cells, int ldim, vtype vt)
{
    vec_struct *tmp_arry;

    if ((cells < 1) || (vt == RVEC && ldim < 1)
	|| (vt == CVEC && ldim < cells) || ldim < 0) {
	G_warning("Vector dimensions out of range.");
	return NULL;
    }

    tmp_arry = (vec_struct *) G_malloc(sizeof(vec_struct));

    if (vt == RVEC) {
	tmp_arry->rows = 1;
	tmp_arry->cols = cells;
	tmp_arry->ldim = ldim;
	tmp_arry->type = ROWVEC_;
    }

    else if (vt == CVEC) {
	tmp_arry->rows = cells;
	tmp_arry->cols = 1;
	tmp_arry->ldim = ldim;
	tmp_arry->type = COLVEC_;
    }

    tmp_arry->v_indx = 0;

    tmp_arry->vals = (doublereal *) G_calloc(ldim * tmp_arry->cols,
					     sizeof(doublereal));
    tmp_arry->is_init = 1;

    return tmp_arry;
}

void *dig__frealloc(void *oldptr, int nelem, int elsize, int oldnelem)
{
    char *ptr;

    if (elsize == 0) {
	elsize = 4;
    }
    if (nelem == 0) {
	nelem = 1;
    }

    ptr = G_calloc(nelem, elsize);

    /*  out of memory  */
    if (!ptr)
	return (ptr);

    {
	register char *a;
	register char *b;
	register size_t n;

	n = oldnelem * elsize;
	a = ptr;
	b = oldptr;
	while (n--)
	    *a++ = *b++;
    }

    G_free(oldptr);
    return (ptr);
}

N_solute_transport_data2d *N_alloc_solute_transport_data2d(int cols, int rows)
{
    N_solute_transport_data2d *data = NULL;

    data =
	(N_solute_transport_data2d *) G_calloc(1,
					       sizeof
					       (N_solute_transport_data2d));

    data->c = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->c_start = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->status = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->diff_x = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->diff_y = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->q = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->cs = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->R = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->nf = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->cin = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->top = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->bottom = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);

    /*Allocate the dispersivity tensor */
    data->disp_xx = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->disp_yy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    data->disp_xy = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);

    data->grad = N_alloc_gradient_field_2d(cols, rows);
    data->stab = N_UPWIND_EXP;

    return data;
}

int sort_areas(struct Map_info *Map, struct line_pnts *Points,
               int field, struct cat_list *cat_list)
{
    int i, centroid, nareas_selected;
    struct line_cats *Cats;
    CELL cat;

    G_begin_polygon_area_calculations();
    Cats = Vect_new_cats_struct();

    /* first count valid areas */
    nareas = Vect_get_num_areas(Map);
    if (nareas == 0)
	return 0;

    /* allocate list to hold valid area info */
    list =
	(struct list *)G_calloc(nareas * sizeof(char), sizeof(struct list));

    /* store area size,cat,index in list */
    nareas_selected = 0;
    for (i = 0; i < nareas; i++) {

	centroid = Vect_get_area_centroid(Map, i + 1);
	SETNULL(&cat);
	if (centroid <= 0) {
	    G_debug(2,_("Area without centroid (OK for island)"));
	}
	else {
	    Vect_read_line(Map, NULL, Cats, centroid);
	    if (field > 0) {
		if (Vect_cats_in_constraint(Cats, field, cat_list)) {
		    Vect_cat_get(Cats, field, &cat);
		    nareas_selected++;
		}
		else {
		    G_debug(2, _("Area centroid without category"));
		}
	    }
	    else {
		/* field < 1, process all areas with centroid */
		cat = 0;
		nareas_selected++;
	    }
	}

	list[i].index = i + 1;
	Vect_get_area_points(Map, i + 1, Points);
	list[i].size =
	    G_area_of_polygon(Points->x, Points->y, Points->n_points);

	list[i].cat = cat;
    }
    if (nareas_selected > 0) {
	/* sort the list by size */
	qsort(list, nareas * sizeof(char), sizeof(struct list), compare);
    }

    return nareas_selected;
}

int
G_bz2_compress(unsigned char *src, int src_sz, unsigned char *dst,
		int dst_sz)
{
    int err;
    unsigned int i, nbytes, buf_sz;
    unsigned char *buf;

#ifndef HAVE_BZLIB_H
    G_fatal_error(_("GRASS needs to be compiled with BZIP2 for BZIP2 compression"));
    return -1;
#else

    /* Catch errors early */
    if (src == NULL || dst == NULL)
	return -1;

    /* Don't do anything if src is empty */
    if (src_sz <= 0)
	return 0;

    /* Output buffer has to be 1% + 600 bytes bigger for single pass compression */
    buf_sz = (unsigned int)((double)dst_sz * 1.01 + (double)600);

    if (NULL == (buf = (unsigned char *)
		 G_calloc(buf_sz, sizeof(unsigned char))))
	return -1;

    /* Do single pass compression */
    nbytes = buf_sz;
    err = BZ2_bzBuffToBuffCompress((char *)buf, &nbytes, /* destination */
                                   (char *)src, src_sz,  /* source */
				   9,			 /* blockSize100k */ 
				   0,                    /* verbosity */
				   0);                   /* workFactor */
    if (err != BZ_OK) {
	G_free(buf);
	return -1;
    }

    /* updated buf_sz is bytes of compressed data */
    if (nbytes >= (unsigned int)src_sz) {
	/* compression not possible */
	G_free(buf);
	return -2;
    }

    /* dst too small */
    if ((unsigned int)dst_sz < nbytes)
	return -2;

    /* Copy the data from buf to dst */
    for (i = 0; i < nbytes; i++)
	dst[i] = buf[i];

    G_free(buf);

    return nbytes;
#endif
}				/* G_bz2_compress() */

/*!
 * \brief Allocate memory for a N_array_3d data structure.
 *
 * This functions allocates an array of type N_array_3d and returns a pointer
 * to the new allocated memory.
 * <br><br>
 * The data type of this array set by "type" must be
 * FCELL_TYPE or DCELL_TYPE accordingly to the raster3d map data types.
 * The offsets sets the number of boundary cols, rows and depths.
 * This option is useful to generate homogeneous Neumann boundary conditions around
 * an array or to establish overlapping boundaries. The arrays are initialized with 0 by default.
 * <br><br>
 * If the offset is greater then 0, negative indices are possible.
 * The data structure of a array with 3 depths, rows and cols and an offset of 1
 * will looks like this:
 *
 \verbatim
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0

 0  0  0  0  0
 0  0  1  2  0
 0  3  4  5  0
 0  6  7  8  0
 0  0  0  0  0

 0  0  0  0  0
 0  9 10 11  0
 0 12 13 14  0
 0 15 16 17  0
 0  0  0  0  0

 0  0  0  0  0
 0 18 19 20  0
 0 21 22 23  0
 0 24 25 26  0
 0  0  0  0  0

 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0

 \endverbatim

 The depth counts from the bottom to the top.

 * <br><br>
 * Internal a one dimensional array is allocated to speed up the memory access.
 * To access the dimensional array with a three dimensional indexing use the provided
 * get and put functions.
 *
 * \param cols int
 * \param rows int
 * \param depths int
 * \param offset int
 * \param type int
 * \return N_array_3d *
 *
 * */
N_array_3d *N_alloc_array_3d(int cols, int rows, int depths, int offset,
			     int type)
{
    N_array_3d *data = NULL;

    if (rows < 1 || cols < 1 || depths < 1)
	G_fatal_error
	    ("N_alloc_array_3d: depths, cols and rows should be > 0");

    if (type != DCELL_TYPE && type != FCELL_TYPE)
	G_fatal_error
	    ("N_alloc_array_3d: Wrong data type, should be FCELL_TYPE or DCELL_TYPE");

    data = (N_array_3d *) G_calloc(1, sizeof(N_array_3d));

    data->cols = cols;
    data->rows = rows;
    data->depths = depths;
    data->type = type;
    data->offset = offset;
    data->rows_intern = rows + 2 * offset;
    data->cols_intern = cols + 2 * offset;
    data->depths_intern = depths + 2 * offset;
    data->fcell_array = NULL;
    data->dcell_array = NULL;

    if (data->type == FCELL_TYPE) {
	data->fcell_array =
	    (float *)G_calloc((size_t) data->depths_intern * data->rows_intern *
			      data->cols_intern, sizeof(float));
	G_debug(3,
		"N_alloc_array_3d: float array allocated rows_intern %i cols_intern %i depths_intern %i offset %i",
		data->rows_intern, data->cols_intern, data->depths_intern,
		data->offset = offset);
    }
    else if (data->type == DCELL_TYPE) {
//.........这里部分代码省略.........

/*!
 * \brief allocate a 9 point star data structure
 *
 * \return N_data_star *
 *
 * \attention The 9 point start is not yet implemented in the matrix assembling function
 *
 * */
N_data_star *N_alloc_9star(void)
{
    N_data_star *star = (N_data_star *) G_calloc(1, sizeof(N_data_star));

    star->type = N_9_POINT_STAR;
    star->count = 9;
    return star;
}

int G_zlib_read(int fd, int rbytes, unsigned char *dst, int nbytes)
{
    int bsize, nread, err;
    unsigned char *b;

    if (dst == NULL || nbytes < 0)
	return -2;

    bsize = rbytes;

    /* Our temporary input buffer for read */
    if (NULL == (b = (unsigned char *)
		 G_calloc(bsize, sizeof(unsigned char))))
	return -1;

    /* Read from the file until we get our bsize or an error */
    nread = 0;
    do {
	err = read(fd, b + nread, bsize - nread);
	if (err >= 0)
	    nread += err;
    } while (err > 0 && nread < bsize);

    /* If the bsize if less than rbytes and we didn't get an error.. */
    if (nread < rbytes && err > 0) {
	G_free(b);
	return -1;
    }

    /* Test if row is compressed */
    if (b[0] == G_ZLIB_COMPRESSED_NO) {
	/* Then just copy it to dst */
	for (err = 0; err < nread - 1 && err < nbytes; err++)
	    dst[err] = b[err + 1];

	G_free(b);
	return (nread - 1);
    }
    else if (b[0] != G_ZLIB_COMPRESSED_YES) {
	/* We're not at the start of a row */
	G_free(b);
	return -1;
    }
    /* Okay it's a compressed row */

    /* Just call G_zlib_expand() with the buffer we read,
     * Account for first byte being a flag
     */
    err = G_zlib_expand(b + 1, bsize - 1, dst, nbytes);

    /* We're done with b */
    G_free(b);

    /* Return whatever G_zlib_expand() returned */
    return err;

}				/* G_zlib_read() */

int can_invert(double **a, int n)
{
    int i, imax = 0, j, k;
    double big, dum, sum, temp;
    double *vv;

    vv = (double *)G_calloc(n, sizeof(double));
    for (i = 0; i < n; i++) {
	big = 0.0;
	for (j = 0; j < n; j++)
	    if ((temp = fabs(a[i][j])) > big)
		big = temp;
	if (big == 0.0)
	    goto singular;
	vv[i] = 1.0 / big;
    }
    for (j = 0; j < n; j++) {
	for (i = 0; i < j; i++) {
	    sum = a[i][j];
	    for (k = 0; k < i; k++)
		sum -= a[i][k] * a[k][j];
	    a[i][j] = sum;
	}
	big = 0.0;
	for (i = j; i < n; i++) {
	    sum = a[i][j];
	    for (k = 0; k < j; k++)
		sum -= a[i][k] * a[k][j];
	    a[i][j] = sum;
	    if ((dum = vv[i] * fabs(sum)) >= big) {
		big = dum;
		imax = i;
	    }
	}
	if (j != imax) {
	    for (k = 0; k < n; k++) {
		dum = a[imax][k];
		a[imax][k] = a[j][k];
		a[j][k] = dum;
	    }
	    vv[imax] = vv[j];
	}
	if (a[j][j] == 0.0)
	    goto singular;
	if (j != n - 1) {
	    dum = 1.0 / (a[j][j]);
	    for (i = j; i < n; i++)
		a[i][j] *= dum;
	}
    }
    G_free(vv);
    return 1;
  singular:
    G_free(vv);
    return 0;
}

void new_stats(char *name, struct Reclass *reclass)
{
    struct Histogram histo, histo2;
    struct Range range;
    CELL cat, cat2;
    int i;
    CELL min, max;

    min = reclass->min;
    max = reclass->max;

    /* read histogram for original file */
    G_suppress_warnings(1);
    i = G_read_histogram(reclass->name, reclass->mapset, &histo);
    G_suppress_warnings(0);
    if (i <= 0)
	return;

    /* compute data rage for reclass */
    G_init_range(&range);

    for (i = 0; i < histo.num; i++) {
	cat = histo.list[i].cat;
	if (cat < min || cat > max)
	    continue;
	cat2 = reclass->table[cat - min];
	G_update_range(cat2, &range);
    }
    G_write_range(name, &range);

    /* now generate a histogram from the original */

    /* allocate histogram list */
    histo2.num += range.max - range.min + 1;

    histo2.list = (LIST *) G_calloc(histo2.num, sizeof(LIST));

    /* set all counts to 0 */
    i = 0;
    for (cat = range.min; cat <= range.max; cat++) {
	histo2.list[i].cat = cat;
	histo2.list[i++].count = 0;
    }

    /* go thru original histogram and add into histo2 */
    for (i = 0; i < histo.num; i++) {
	cat = histo.list[i].cat;
	if (cat < min || cat > max)
	    G_set_c_null_value(&cat, 1);
	else
	    cat2 = reclass->table[cat - min];
	if (!G_is_c_null_value(&cat))
	    histo2.list[cat2 - range.min].count += histo.list[i].count;
    }
    G_write_histogram(name, &histo2);
}

/*!
 * \brief Allocate the pde geometry data structure and return a pointer to the new allocated structure
 *
 * \return N_geom_data *
 * */
N_geom_data *N_alloc_geom_data(void)
{
    N_geom_data *geom = (N_geom_data *) G_calloc(1, sizeof(N_geom_data));

    geom->area = NULL;
    geom->planimetric = 1;
    geom->dim = 0;

    return geom;
}