[GRASS-SVN] r38891 - in grass/branches/develbranch_6: imagery/i.cca
imagery/i.pca imagery/i.smap/bouman raster/r.gwflow
raster3d/r3.gwflow
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 27 16:58:55 EDT 2009
Author: huhabla
Date: 2009-08-27 16:58:55 -0400 (Thu, 27 Aug 2009)
New Revision: 38891
Modified:
grass/branches/develbranch_6/imagery/i.cca/local_proto.h
grass/branches/develbranch_6/imagery/i.cca/main.c
grass/branches/develbranch_6/imagery/i.cca/matrix.c
grass/branches/develbranch_6/imagery/i.cca/stats.c
grass/branches/develbranch_6/imagery/i.cca/transform.c
grass/branches/develbranch_6/imagery/i.pca/main.c
grass/branches/develbranch_6/imagery/i.smap/bouman/model.c
grass/branches/develbranch_6/raster/r.gwflow/Makefile
grass/branches/develbranch_6/raster/r.gwflow/main.c
grass/branches/develbranch_6/raster3d/r3.gwflow/Makefile
grass/branches/develbranch_6/raster3d/r3.gwflow/main.c
Log:
Patching the modules to use new functions from gmath and
ccmath wrapper libraries
Modified: grass/branches/develbranch_6/imagery/i.cca/local_proto.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.cca/local_proto.h 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/imagery/i.cca/local_proto.h 2009-08-27 20:58:55 UTC (rev 38891)
@@ -1,22 +1,19 @@
#ifndef __LOCAL_PROTO_H__
#define __LOCAL_PROTO_H__
-#define MX 9
-#define MC 50
-
/* matrix.c */
-int product(double[MX], double, double[MX][MX], int);
-int setdiag(double[MX], int, double[MX][MX]);
-int getsqrt(double[MX][MX], int, double[MX][MX], double[MX][MX]);
-int solveq(double[MX][MX], int, double[MX][MX], double[MX][MX]);
-int matmul(double[MX][MX], double[MX][MX], double[MX][MX], int);
+int product(double*, double, double**, int);
+int setdiag(double*, int, double**);
+int getsqrt(double**, int, double**, double**);
+int solveq(double**, int, double**, double**);
+int print_matrix(double **matrix, int bands);
/* stats.c */
-int within(int, int, double[MC], double[MC][MX][MX], double[MX][MX], int);
-int between(int, int, double[MC], double[MC][MX], double[MX][MX], int);
+int within(int, int, double*, double***, double**, int);
+int between(int, int, double*, double**, double**, int);
/* transform.c */
-int transform(int[MX], int[MX], int, int, double[MX][MX], int, CELL[MX],
- CELL[MX]);
+int transform(int*, int*, int, int, double**, int, CELL*,
+ CELL*);
#endif /* __LOCAL_PROTO_H__ */
Modified: grass/branches/develbranch_6/imagery/i.cca/main.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.cca/main.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/imagery/i.cca/main.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -53,26 +53,26 @@
int bands; /* Number of image bands */
int nclass; /* Number of classes */
int samptot; /* Total number of sample points */
- double mu[MC][MX]; /* Mean vector for image classes */
- double w[MX][MX]; /* Within Class Covariance Matrix */
- double p[MX][MX]; /* Between class Covariance Matrix */
- double l[MX][MX]; /* Diagonal matrix of eigenvalues */
- double q[MX][MX]; /* Transformation matrix */
- double cov[MC][MX][MX]; /* Individual class Covariance Matrix */
- double nsamp[MC]; /* Number of samples in a given class */
- double eigval[MX]; /* Eigen value vector */
- double eigmat[MX][MX]; /* Eigen Matrix */
- char tempname[50];
+ double **mu; /* Mean vector for image classes */
+ double **w; /* Within Class Covariance Matrix */
+ double **p; /* Between class Covariance Matrix */
+ double **l; /* Diagonal matrix of eigenvalues */
+ double **q; /* Transformation matrix */
+ double ***cov; /* Individual class Covariance Matrix */
+ double *nsamp; /* Number of samples in a given class */
+ double *eigval; /* Eigen value vector */
+ double **eigmat; /* Eigen Matrix */
+ char tempname[1024];
/* used to make the color tables */
- CELL outbandmax[MX]; /* will hold the maximums found in the out maps */
- CELL outbandmin[MX]; /* will hold the minimums found in the out maps */
+ CELL *outbandmax; /* will hold the maximums found in the out maps */
+ CELL *outbandmin; /* will hold the minimums found in the out maps */
struct Colors color_tbl;
struct Signature sigs;
FILE *sigfp;
struct Ref refs;
- int datafds[MX];
- int outfds[MX];
+ int *datafds;
+ int *outfds;
struct GModule *module;
struct Option *grp_opt, *subgrp_opt, *sig_opt, *out_opt;
@@ -141,10 +141,30 @@
/* check the number of input bands */
bands = refs.nfiles;
- if (bands > MX - 1)
- G_fatal_error(_("Subgroup too large. Maximum number of bands is %d\n."),
- MX - 1);
+
+ /*memory allocation*/
+ mu = G_alloc_matrix(nclass, bands);
+ w = G_alloc_matrix(bands, bands);
+ p = G_alloc_matrix(bands, bands);
+ l = G_alloc_matrix(bands, bands);
+ q = G_alloc_matrix(bands, bands);
+ eigmat = G_alloc_matrix(bands, bands);
+ nsamp = G_alloc_vector(nclass);
+ eigval = G_alloc_vector(bands);
+
+ cov = (double***)G_calloc(nclass, sizeof(double**));
+ for(i = 0; i < nclass; i++)
+ {
+ cov[i] = G_alloc_matrix(bands,bands);
+ }
+
+ outbandmax = (CELL*)G_calloc(nclass, sizeof(CELL));
+ outbandmin = (CELL*)G_calloc(nclass, sizeof(CELL));
+ datafds = (int*)G_calloc(nclass, sizeof(int));
+ outfds = (int*)G_calloc(nclass, sizeof(int));
+
+
/*
Here is where the information regarding
a) Number of samples per class
@@ -153,40 +173,53 @@
*/
samptot = 0;
- for (i = 1; i <= nclass; i++) {
- nsamp[i] = sigs.sig[i - 1].npoints;
+ for (i = 0; i < nclass; i++) {
+ nsamp[i] = sigs.sig[i].npoints;
samptot += nsamp[i];
- for (j = 1; j <= bands; j++) {
- mu[i][j] = sigs.sig[i - 1].mean[j - 1];
- for (k = 1; k <= j; k++)
+ for (j = 0; j < bands; j++) {
+ mu[i][j] = sigs.sig[i].mean[j];
+ for (k = 0; k <= j; k++)
+ {
cov[i][j][k] = cov[i][k][j] =
- sigs.sig[i - 1].var[j - 1][k - 1];
+ sigs.sig[i].var[j][k];
+ }
}
}
within(samptot, nclass, nsamp, cov, w, bands);
between(samptot, nclass, nsamp, mu, p, bands);
- jacobi(w, (long)bands, eigval, eigmat);
- egvorder(eigval, eigmat, (long)bands);
+ G_math_d_copy(w[0], eigmat[0], bands*bands);
+ G_math_eigen(eigmat, eigval, bands);
+ G_math_egvorder(eigval, eigmat, bands);
setdiag(eigval, bands, l);
getsqrt(w, bands, l, eigmat);
solveq(q, bands, w, p);
- jacobi(q, (long)bands, eigval, eigmat);
- egvorder(eigval, eigmat, (long)bands);
- matmul(q, eigmat, w, bands);
+ G_math_d_copy(q[0], eigmat[0], bands*bands);
+ G_math_eigen(eigmat, eigval, bands);
+ G_math_egvorder(eigval, eigmat, bands);
+ G_math_d_AB(eigmat, w, q, bands, bands, bands);
+ for(i = 0; i < bands; i++)
+ {
+ G_verbose_message("%i. eigen value: %+6.5f", i, eigval[i]);
+ G_verbose_message("eigen vector:");
+ for(j = 0; j < bands; j++)
+ G_verbose_message("%+6.5f ", eigmat[i][j]);
+
+ }
+
/* open the cell maps */
- for (i = 1; i <= bands; i++) {
+ for (i = 0; i < bands; i++) {
outbandmax[i] = (CELL) 0;
outbandmin[i] = (CELL) 0;
- if ((datafds[i] = G_open_cell_old(refs.file[i - 1].name,
- refs.file[i - 1].mapset)) < 0) {
+ if ((datafds[i] = G_open_cell_old(refs.file[i].name,
+ refs.file[i].mapset)) < 0) {
G_fatal_error(_("Cannot open raster map <%s>"),
- refs.file[i - 1].name);
+ refs.file[i].name);
}
- sprintf(tempname, "%s.%d", out_opt->answer, i);
+ sprintf(tempname, "%s.%d", out_opt->answer, i + 1);
if ((outfds[i] = G_open_cell_new(tempname)) < 0)
G_fatal_error(_("Cannot create raster map <%s>"), tempname);
}
@@ -199,17 +232,18 @@
G_init_colors(&color_tbl);
/* close the cell maps */
- for (i = 1; i <= bands; i++) {
+ for (i = 0; i < bands; i++) {
G_close_cell(datafds[i]);
G_close_cell(outfds[i]);
if (outbandmin[i] < (CELL) 0 || outbandmax[i] > (CELL) 255) {
G_warning(_("The output cell map <%s.%d> has values "
- "outside the 0-255 range."), out_opt->answer, i);
+ "outside the 0-255 range. Min: %d Max: %d"),
+ out_opt->answer, i + 1, outbandmin[i], outbandmax[i]);
}
G_make_grey_scale_colors(&color_tbl, 0, outbandmax[i]);
- sprintf(tempname, "%s.%d", out_opt->answer, i);
+ sprintf(tempname, "%s.%d", out_opt->answer, i + 1);
/* write a color table */
G_write_colors(tempname, G_mapset(), &color_tbl);
@@ -218,5 +252,26 @@
I_free_signatures(&sigs);
I_free_group_ref(&refs);
+
+ /*free memory*/
+ G_free_matrix(mu);
+ G_free_matrix(w);
+ G_free_matrix(p);
+ G_free_matrix(l);
+ G_free_matrix(q);
+ G_free_matrix(eigmat);
+ for(i = 0; i < nclass; i++)
+ G_free_matrix(cov[i]);
+ G_free(cov);
+
+ G_free_vector(nsamp);
+ G_free_vector(eigval);
+
+ G_free(outbandmax);
+ G_free(outbandmin);
+ G_free(datafds);
+ G_free(outfds);
+
+
exit(EXIT_SUCCESS);
}
Modified: grass/branches/develbranch_6/imagery/i.cca/matrix.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.cca/matrix.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/imagery/i.cca/matrix.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -3,26 +3,39 @@
#include <grass/gmath.h>
#include "local_proto.h"
+int print_matrix(double **matrix, int bands)
+{
+ int i, j;
-int product(double vector[MX], double factor, double matrix1[MX][MX],
+ for (i = 0; i < bands; i++)
+ {
+ for (j = 0; j < bands; j++) {
+ printf("%g ", matrix[i][j]);
+ }
+ printf("\n");
+ }
+ return 0;
+}
+
+int product(double *vector, double factor, double **matrix1,
int bands)
{
int i, j;
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++) {
+ for (i = 0; i < bands; i++)
+ for (j = 0; j < bands; j++) {
matrix1[i][j] = (double)factor *(vector[i] * vector[j]);
}
return 0;
}
-int setdiag(double eigval[MX], int bands, double l[MX][MX])
+int setdiag(double *eigval, int bands, double **l)
{
int i, j;
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++)
+ for (i = 0; i < bands; i++)
+ for (j = 0; j < bands; j++)
if (i == j)
l[i][j] = eigval[i];
else
@@ -32,43 +45,36 @@
int
-getsqrt(double w[MX][MX], int bands, double l[MX][MX], double eigmat[MX][MX])
+getsqrt(double **w, int bands, double **l, double **eigmat)
{
int i;
- double tmp[MX][MX];
+ double **tmp;
- for (i = 1; i <= bands; i++)
+ tmp = G_alloc_matrix(bands, bands);
+
+ for (i = 0; i < bands; i++)
l[i][i] = 1.0 / sqrt(l[i][i]);
- matmul(tmp, eigmat, l, bands);
- transpose(eigmat, bands);
- matmul(w, tmp, eigmat, bands);
+
+ G_math_d_AB(eigmat, l, tmp, bands, bands, bands);
+ G_math_d_A_T(eigmat, bands);
+ G_math_d_AB(tmp, eigmat, w, bands, bands, bands);
+
+ G_free_matrix(tmp);
+
return 0;
}
-int solveq(double q[MX][MX], int bands, double w[MX][MX], double p[MX][MX])
+int solveq(double **q, int bands, double **w, double **p)
{
- double tmp[MX][MX];
+ double **tmp;
- matmul(tmp, w, p, bands);
- matmul(q, tmp, w, bands);
- return 0;
-}
+ tmp = G_alloc_matrix(bands, bands);
+ G_math_d_AB(w, p, tmp, bands, bands, bands);
+ G_math_d_AB(tmp, w, q, bands, bands, bands);
-int matmul(double res[MX][MX], double m1[MX][MX], double m2[MX][MX], int dim)
-{
- int i, j, k;
- double sum;
+ G_free_matrix(tmp);
- for (i = 1; i <= dim; i++) {
- for (j = 1; j <= dim; j++) {
- sum = 0.0;
- for (k = 1; k <= dim; k++)
- sum += m1[i][k] * m2[k][j];
- res[i][j] = sum;
- }
- }
-
return 0;
}
Modified: grass/branches/develbranch_6/imagery/i.cca/stats.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.cca/stats.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/imagery/i.cca/stats.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -3,23 +3,23 @@
int
-within(int samptot, int nclass, double nsamp[MC], double cov[MC][MX][MX],
- double w[MX][MX], int bands)
+within(int samptot, int nclass, double *nsamp, double ***cov,
+ double **w, int bands)
{
int i, j, k;
/* Initialize within class covariance matrix */
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++)
+ for (i = 0; i < bands; i++)
+ for (j = 0; j < bands; j++)
w[i][j] = 0.0;
- for (i = 1; i <= nclass; i++)
- for (j = 1; j <= bands; j++)
- for (k = 1; k <= bands; k++)
+ for (i = 0; i < nclass; i++)
+ for (j = 0; j < bands; j++)
+ for (k = 0; k < bands; k++)
w[j][k] += (nsamp[i] - 1) * cov[i][j][k];
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++)
+ for (i = 0; i < bands; i++)
+ for (j = 0; j < bands; j++)
w[i][j] = (1.0 / ((double)(samptot - nclass))) * w[i][j];
return 0;
@@ -27,50 +27,40 @@
int
-between(int samptot, int nclass, double nsamp[MC], double mu[MC][MX],
- double p[MX][MX], int bands)
+between(int samptot, int nclass, double *nsamp, double **mu,
+ double **p, int bands)
{
int i, j, k;
- double tmp0[MX][MX], tmp1[MX][MX], tmp2[MX][MX];
- double newvec[MX];
+ double **tmp0, **tmp1, **tmp2;
+ double *newvec;
- for (i = 0; i < MX; i++)
- newvec[i] = 0.0;
+ tmp0 = G_alloc_matrix(bands, bands);
+ tmp1 = G_alloc_matrix(bands, bands);
+ tmp2 = G_alloc_matrix(bands, bands);
+ newvec = G_alloc_vector(bands);
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++)
- tmp1[i][j] = tmp2[i][j] = 0.0;
-
- /* for (i = 1 ; i <= nclass ; i++)
- product(mu[i],nsamp[i],tmp0,tmp1,bands);
- for (i = 1 ; i <= nclass ; i++)
- for (j = 1 ; j <= bands ; j++)
- newvec[j] += nsamp[i] * mu[i][j];
- for (i = 1 ; i <= bands ; i++)
- for (j = 1 ; i <= bands ; j++)
- tmp2[i][j] = (newvec[i] * newvec[j]) / samptot;
- for (i = 1 ; i <= bands ; i++)
- for (j = 1 ; j <= bands ; j++)
- p[i][j] = (tmp1[i][j] - tmp2[i][j]) / (nclass - 1);
- */
-
- for (i = 1; i <= nclass; i++)
- for (j = 1; j <= bands; j++)
+ for (i = 0; i < nclass; i++)
+ for (j = 0; j < bands; j++)
newvec[j] += nsamp[i] * mu[i][j];
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++)
+ for (i = 0; i < bands; i++)
+ for (j = 0; j < bands; j++)
tmp1[i][j] = (newvec[i] * newvec[j]) / samptot;
- for (k = 1; k <= nclass; k++) {
+ for (k = 0; k < nclass; k++) {
product(mu[k], nsamp[k], tmp0, bands);
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++)
+ for (i = 0; i < bands; i++)
+ for (j = 0; j < bands; j++)
tmp2[i][j] += tmp0[i][j] - tmp1[i][j];
}
- for (i = 1; i <= bands; i++)
- for (j = 1; j <= bands; j++)
+ for (i = 0; i < bands; i++)
+ for (j = 0; j < bands; j++)
p[i][j] = tmp2[i][j] / (nclass - 1);
+ G_free_matrix(tmp0);
+ G_free_matrix(tmp1);
+ G_free_matrix(tmp2);
+ G_free_vector(newvec);
+
return 0;
}
Modified: grass/branches/develbranch_6/imagery/i.cca/transform.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.cca/transform.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/imagery/i.cca/transform.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -5,33 +5,37 @@
int
-transform(int datafds[MX], int outfds[MX], int rows, int cols,
- double eigmat[MX][MX], int bands, CELL mins[MX], CELL maxs[MX])
+transform(int *datafds, int *outfds, int rows, int cols,
+ double **eigmat, int bands, CELL *mins, CELL *maxs)
{
int i, j, k, l;
- double sum[MX];
- CELL *rowbufs[MX];
+ double *sum;
+ CELL **rowbufs;
+ sum = G_alloc_vector(bands);
+ rowbufs = (CELL**)G_calloc(bands, sizeof(CELL*));
+
+
/* allocate row buffers for each band */
- for (i = 1; i <= bands; i++)
+ for (i = 0; i < bands; i++)
if ((rowbufs[i] = G_allocate_cell_buf()) == NULL)
G_fatal_error(_("Unable to allocate cell buffers."));
for (i = 0; i < rows; i++) {
/* get one row of data */
- for (j = 1; j <= bands; j++)
+ for (j = 0; j < bands; j++)
if (G_get_map_row(datafds[j], rowbufs[j], i) < 0)
G_fatal_error(_("Error reading cell map during transform."));
/* transform each cell in the row */
for (l = 0; l < cols; l++) {
- for (j = 1; j <= bands; j++) {
+ for (j = 0; j < bands; j++) {
sum[j] = 0.0;
- for (k = 1; k <= bands; k++) {
+ for (k = 0; k < bands; k++) {
sum[j] += eigmat[j][k] * (double)rowbufs[k][l];
}
}
- for (j = 1; j <= bands; j++) {
+ for (j = 0; j < bands; j++) {
rowbufs[j][l] = (CELL) (sum[j] + 0.5);
if (rowbufs[j][l] > maxs[j])
maxs[j] = rowbufs[j][l];
@@ -41,13 +45,16 @@
}
/* output the row of data */
- for (j = 1; j <= bands; j++)
+ for (j = 0; j < bands; j++)
if (G_put_raster_row(outfds[j], rowbufs[j], CELL_TYPE) < 0)
G_fatal_error(_("Error writing cell map during transform."));
}
- for (i = 1; i <= bands; i++)
+ for (i = 0; i < bands; i++)
G_free(rowbufs[i]);
+ G_free(rowbufs);
+ G_free_vector(sum);
+
G_message(_("Transform completed.\n"));
return 0;
Modified: grass/branches/develbranch_6/imagery/i.pca/main.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.pca/main.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/imagery/i.pca/main.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -101,22 +101,13 @@
set_output_scale(opt_scale, &scale, &scale_min, &scale_max);
/* allocate memory */
- covar = (double **)G_calloc(bands, sizeof(double *));
- mu = (double *)G_malloc(bands * sizeof(double));
- inp_fd = (int *)G_malloc(bands * sizeof(int));
- eigmat = (double **)G_calloc(bands, sizeof(double *));
- eigval = (double *)G_calloc(bands, sizeof(double));
+ covar = G_alloc_matrix(bands, bands);
+ mu = G_alloc_vector(bands);
+ inp_fd = G_alloc_ivector(bands);
+ eigmat = G_alloc_matrix(bands, bands);
+ eigval = G_alloc_vector(bands);
+
- /* allocate memory for matrices */
- for (i = 0; i < bands; i++) {
- covar[i] = (double *)G_malloc(bands * sizeof(double));
- eigmat[i] = (double *)G_calloc(bands, sizeof(double));
-
- /* initialize covariance matrix */
- for (j = 0; j < bands; j++)
- covar[i][j] = 0.;
- }
-
/* open and check input/output files */
for (i = 0; i < bands; i++) {
char tmpbuf[128];
@@ -147,8 +138,9 @@
}
}
+ G_math_d_copy(covar[0], eigmat[0], bands*bands);
G_debug(1, "Calculating eigenvalues and eigenvectors...");
- eigen(covar, eigmat, eigval, bands);
+ G_math_eigen(eigmat, eigval, bands);
#ifdef PCA_DEBUG
/* dump eigen matrix and eigen values */
@@ -156,10 +148,10 @@
#endif
G_debug(1, "Ordering eigenvalues in descending order...");
- egvorder2(eigval, eigmat, bands);
+ G_math_egvorder(eigval, eigmat, bands);
G_debug(1, "Transposing eigen matrix...");
- transpose2(eigmat, bands);
+ G_math_d_A_T(eigmat, bands);
/* write output images */
write_pca(eigmat, inp_fd, opt_out->answer, bands, scale, scale_min,
@@ -177,7 +169,12 @@
/* close output file */
G_unopen_cell(inp_fd[i]);
}
-
+ /* allocate memory */
+ G_free_matrix(covar);
+ G_free_vector(mu);
+ G_free_ivector(inp_fd);
+ G_free_matrix(eigmat);
+ G_free_vector(eigval);
exit(EXIT_SUCCESS);
}
Modified: grass/branches/develbranch_6/imagery/i.smap/bouman/model.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.smap/bouman/model.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/imagery/i.smap/bouman/model.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -43,7 +43,7 @@
}
/* Test for positive definite matrix */
- eigen(SubS->Rinv, NULL, lambda, nbands);
+ G_math_eigval(SubS->Rinv, lambda, nbands);
for (b1 = 0; b1 < nbands; b1++) {
if (lambda[b1] <= 0.0)
G_warning(_("Nonpositive eigenvalues for class %d subclass %d"),
Modified: grass/branches/develbranch_6/raster/r.gwflow/Makefile
===================================================================
--- grass/branches/develbranch_6/raster/r.gwflow/Makefile 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/raster/r.gwflow/Makefile 2009-08-27 20:58:55 UTC (rev 38891)
@@ -2,8 +2,8 @@
PGM=r.gwflow
-LIBES = $(GISLIB) $(GPDELIB)
-DEPENDENCIES = $(GISDEP) $(GPDEDEP)
+LIBES = $(GISLIB) $(GMATHLIB) $(GPDELIB)
+DEPENDENCIES = $(GISDEP) $(GMATHDEP) $(GPDEDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass/branches/develbranch_6/raster/r.gwflow/main.c
===================================================================
--- grass/branches/develbranch_6/raster/r.gwflow/main.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/raster/r.gwflow/main.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -23,13 +23,14 @@
#include <grass/glocale.h>
#include <grass/N_pde.h>
#include <grass/N_gwflow.h>
+#include <grass/gmath.h>
/*- Parameters and global variables -----------------------------------------*/
typedef struct
{
struct Option *output, *phead, *status, *hc_x, *hc_y, *q, *s, *r, *top,
- *bottom, *vector, *type, *dt, *maxit, *error, *solver, *sor,
+ *bottom, *vector, *type, *dt, *maxit, *error, *solver,
*river_head, *river_bed, *river_leak, *drain_bed, *drain_leak;
struct Flag *sparse;
} paramType;
@@ -37,13 +38,12 @@
paramType param; /*Parameters */
/*- prototypes --------------------------------------------------------------*/
-void set_params(void); /*Fill the paramType structure */
-void copy_result(N_array_2d * status, N_array_2d * phead_start,
- double *result, struct Cell_head *region,
- N_array_2d * target);
-N_les *create_solve_les(N_geom_data * geom, N_gwflow_data2d * data,
- N_les_callback_2d * call, const char *solver,
- int maxit, double error, double sor);
+static void set_params(void); /*Fill the paramType structure */
+static void copy_result(N_array_2d * status, N_array_2d * phead_start, double *result,
+ struct Cell_head *region, N_array_2d * target);
+static N_les *create_solve_les(N_geom_data * geom, N_gwflow_data2d * data,
+ N_les_callback_2d * call, const char *solver, int maxit,
+ double error);
/* ************************************************************************* */
/* Set up the arguments we are expecting ********************************** */
@@ -62,8 +62,7 @@
param.status->type = TYPE_STRING;
param.status->required = YES;
param.status->gisprompt = "old,raster,raster";
- param.status->description =
- _("Boundary condition status, 0-inactive, 1-active, 2-dirichlet");
+ param.status->description = _("Boundary condition status, 0-inactive, 1-active, 2-dirichlet");
param.hc_x = G_define_option();
param.hc_x->key = "hc_x";
@@ -122,16 +121,15 @@
param.output->type = TYPE_STRING;
param.output->required = YES;
param.output->gisprompt = "new,raster,raster";
- param.output->description = _("The map storing the numerical result [m]");
+ param.output->description = _("The map storing the numerical result [m]");
param.vector = G_define_option();
param.vector->key = "velocity";
param.vector->type = TYPE_STRING;
param.vector->required = NO;
param.vector->gisprompt = "new,raster,raster";
- param.vector->description =
- _("Calculate the groundwater filter velocity vector field [m/s]\n"
- "and write the x, and y components to maps named name_[xy]");
+ param.vector->description = _("Calculate the groundwater filter velocity vector field [m/s]\n"
+ "and write the x, and y components to maps named name_[xy]");
param.type = G_define_option();
param.type->key = "type";
@@ -147,7 +145,7 @@
param.river_bed->type = TYPE_STRING;
param.river_bed->required = NO;
param.river_bed->gisprompt = "old,raster,raster";
- param.river_bed->description = _("The height of the river bed in [m]");
+ param.river_bed->description = _("The hight of the river bed in [m]");
param.river_head = G_define_option();
param.river_head->key = "river_head";
@@ -170,26 +168,24 @@
param.drain_bed->type = TYPE_STRING;
param.drain_bed->required = NO;
param.drain_bed->gisprompt = "old,raster,raster";
- param.drain_bed->description = _("The height of the drainage bed in [m]");
+ param.drain_bed->description = _("The hight of the drainage bed in [m]");
param.drain_leak = G_define_option();
param.drain_leak->key = "drain_leak";
param.drain_leak->type = TYPE_STRING;
param.drain_leak->required = NO;
param.drain_leak->gisprompt = "old,raster,raster";
- param.drain_leak->description =
- _("The leakage coefficient of the drainage bed in [1/s]");
-
+ param.drain_leak->description = _("The leakage coefficient of the drainage bed in [1/s]");
+
param.dt = N_define_standard_option(N_OPT_CALC_TIME);
param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
param.error = N_define_standard_option(N_OPT_ITERATION_ERROR);
param.solver = N_define_standard_option(N_OPT_SOLVER_SYMM);
- param.sor = N_define_standard_option(N_OPT_SOR_VALUE);
+ param.solver->options = "cg,pcg,cholesky";
param.sparse = G_define_flag();
param.sparse->key = 's';
- param.sparse->description =
- _("Use a sparse matrix, only available with iterative solvers");
+ param.sparse->description = _("Use a sparse matrix, only available with iterative solvers");
}
@@ -205,7 +201,7 @@
N_les_callback_2d *call = NULL;
double *tmp_vect = NULL;
struct Cell_head region;
- double error, sor, max_norm = 0, tmp;
+ double error, max_norm = 0, tmp;
int maxit, i, inner_count = 0;
char *solver;
int x, y, stat;
@@ -221,8 +217,7 @@
module = G_define_module();
module->keywords = _("raster");
- module->description =
- _("Numerical calculation program for transient, confined and unconfined groundwater flow in two dimensions.");
+ module->description = _("Numerical calculation program for transient, confined and unconfined groundwater flow in two dimensions.");
/* Get parameters from user */
set_params();
@@ -236,9 +231,8 @@
param.river_head->answer == NULL) {
with_river = 0;
}
- else if (param.river_leak->answer != NULL &&
- param.river_bed->answer != NULL &&
- param.river_head->answer != NULL) {
+ else if (param.river_leak->answer != NULL && param.river_bed->answer != NULL
+ && param.river_head->answer != NULL) {
with_river = 1;
}
else {
@@ -263,7 +257,6 @@
sscanf(param.maxit->answer, "%i", &(maxit));
/*Set the calculation error break criteria */
sscanf(param.error->answer, "%lf", &(error));
- sscanf(param.sor->answer, "%lf", &(sor));
/*set the solver */
solver = param.solver->answer;
@@ -363,11 +356,10 @@
/*assemble the linear equation system and solve it */
- les = create_solve_les(geom, data, call, solver, maxit, error, sor);
+ les = create_solve_les(geom, data, call, solver, maxit, error);
/* copy the result into the phead array for output or unconfined calculation */
- copy_result(data->status, data->phead_start, les->x, ®ion,
- data->phead);
+ copy_result(data->status, data->phead_start, les->x, ®ion, data->phead);
N_convert_array_2d_null_to_zero(data->phead);
/****************************************************/
@@ -387,7 +379,7 @@
inner_count = 0;
do {
- G_message(_("Calculation of unconfined groundwater flow loop %i\n"),
+ G_message(_("Calculation of unconfined groundwater flow: loop %i\n"),
inner_count + 1);
/* we will allocate a new les for each loop */
@@ -395,8 +387,7 @@
N_free_les(les);
/*assemble the linear equation system and solve it */
- les =
- create_solve_les(geom, data, call, solver, maxit, error, sor);
+ les = create_solve_les(geom, data, call, solver, maxit, error);
/*calculate the maximum norm of the groundwater height difference */
tmp = 0;
@@ -419,7 +410,7 @@
N_convert_array_2d_null_to_zero(data->phead);
/**/ inner_count++;
}
- while (max_norm > 0.01 && inner_count < 50);
+ while (max_norm > 0.01 && inner_count < 500);
if (tmp_vect)
free(tmp_vect);
@@ -514,49 +505,40 @@
/* ***** create and solve the linear equation system ************* */
/* *************************************************************** */
N_les *create_solve_les(N_geom_data * geom, N_gwflow_data2d * data,
- N_les_callback_2d * call, const char *solver,
- int maxit, double error, double sor)
+ N_les_callback_2d * call, const char *solver, int maxit,
+ double error)
{
N_les *les;
/*assemble the linear equation system */
if (param.sparse->answer)
- les =
- N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status,
- data->phead, (void *)data, call);
+ les = N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status, data->phead, (void *)data, call);
else
- les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead, (void *)data, call);
+ les = N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead, (void *)data, call);
N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead);
- /*solve the equation system */
- if (strcmp(solver, N_SOLVER_ITERATIVE_JACOBI) == 0)
- N_solver_jacobi(les, maxit, sor, error);
-
- if (strcmp(solver, N_SOLVER_ITERATIVE_SOR) == 0)
- N_solver_SOR(les, maxit, sor, error);
-
+ /*solve the linear equation system */
+ if(les && les->type == G_MATH_NORMAL_LES)
+ {
if (strcmp(solver, N_SOLVER_ITERATIVE_CG) == 0)
- N_solver_cg(les, maxit, error);
+ G_math_solver_cg(les->A, les->x, les->b, les->rows, maxit, error);
if (strcmp(solver, N_SOLVER_ITERATIVE_PCG) == 0)
- N_solver_pcg(les, maxit, error, N_DIAGONAL_PRECONDITION);
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, maxit, error, N_DIAGONAL_PRECONDITION);
- if (strcmp(solver, N_SOLVER_ITERATIVE_BICGSTAB) == 0)
- N_solver_bicgstab(les, maxit, error);
-
- if (strcmp(solver, N_SOLVER_DIRECT_LU) == 0)
- N_solver_lu(les);
-
if (strcmp(solver, N_SOLVER_DIRECT_CHOLESKY) == 0)
- N_solver_cholesky(les);
+ G_math_solver_cholesky(les->A, les->x, les->b, les->rows, les->rows);
+ } else if (les && les->type == G_MATH_SPARSE_LES)
+ {
+ if (strcmp(solver, N_SOLVER_ITERATIVE_CG) == 0)
+ G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows, maxit, error);
- if (strcmp(solver, N_SOLVER_DIRECT_GAUSS) == 0)
- N_solver_gauss(les);
+ if (strcmp(solver, N_SOLVER_ITERATIVE_PCG) == 0)
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, maxit, error, N_DIAGONAL_PRECONDITION);
+ }
if (les == NULL)
G_fatal_error(_("Unable to create and solve the linear equation system"));
Modified: grass/branches/develbranch_6/raster3d/r3.gwflow/Makefile
===================================================================
--- grass/branches/develbranch_6/raster3d/r3.gwflow/Makefile 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/raster3d/r3.gwflow/Makefile 2009-08-27 20:58:55 UTC (rev 38891)
@@ -2,8 +2,8 @@
PGM=r3.gwflow
-LIBES = $(GPDELIB) $(G3DLIB) $(GISLIB)
-DEPENDENCIES = $(GPDEDEP) $(G3DDEP) $(GISDEP)
+LIBES = $(GMATHLIB) $(GPDELIB) $(G3DLIB) $(GISLIB)
+DEPENDENCIES = $(GMATHDEP) $(GPDEDEP) $(G3DDEP) $(GISDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Modified: grass/branches/develbranch_6/raster3d/r3.gwflow/main.c
===================================================================
--- grass/branches/develbranch_6/raster3d/r3.gwflow/main.c 2009-08-27 20:57:08 UTC (rev 38890)
+++ grass/branches/develbranch_6/raster3d/r3.gwflow/main.c 2009-08-27 20:58:55 UTC (rev 38891)
@@ -20,6 +20,7 @@
#include <string.h>
#include <grass/gis.h>
#include <grass/G3d.h>
+#include <grass/gmath.h>
#include <grass/glocale.h>
#include <grass/N_pde.h>
#include <grass/N_gwflow.h>
@@ -29,7 +30,7 @@
typedef struct
{
struct Option *output, *phead, *status, *hc_x, *hc_y, *hc_z, *q, *s, *r,
- *vector, *dt, *maxit, *error, *solver, *sor;
+ *vector, *dt, *maxit, *error, *solver;
struct Flag *mask;
struct Flag *sparse;
} paramType;
@@ -37,11 +38,12 @@
paramType param; /*Parameters */
/*- prototypes --------------------------------------------------------------*/
-void set_params(void); /*Fill the paramType structure */
-void write_result(N_array_3d * status, N_array_3d * phead_start,
- N_array_3d * phead, double *result, G3D_Region * region,
- char *name);
+static void set_params(void); /*Fill the paramType structure */
+static void write_result(N_array_3d * status, N_array_3d * phead_start,
+ N_array_3d * phead, double *result,
+ G3D_Region * region, char *name);
+
/* ************************************************************************* */
/* Set up the arguments we are expecting ********************************** */
/* ************************************************************************* */
@@ -60,7 +62,8 @@
param.status->required = YES;
param.status->gisprompt = "old,grid3,3d-raster";
param.status->description =
- _("The status for each cell, = 0 - inactive, 1 - active, 2 - dirichlet");
+ _
+ ("The status for each cell, = 0 - inactive, 1 - active, 2 - dirichlet");
param.hc_x = G_define_option();
param.hc_x->key = "hc_x";
@@ -112,23 +115,23 @@
param.output->type = TYPE_STRING;
param.output->required = YES;
param.output->gisprompt = "new,grid3,3d-raster";
- param.output->description =
- _("The piezometric head result of the numerical calculation will be written to this map");
+ param.output->description = _("The piezometric head result of the numerical calculation will be written to this map");
param.vector = G_define_option();
param.vector->key = "velocity";
param.vector->type = TYPE_STRING;
param.vector->required = NO;
param.vector->gisprompt = "new,grid3,3d-raster";
- param.vector->description =
- _("Calculate the groundwater distance velocity vector field and write the x, y, and z components to maps named name_[xyz]. Name is basename for the new raster3d maps");
+ param.vector->description = _("Calculate the groundwater distance velocity vector field \n"
+ "and write the x, y, and z components to maps named name_[xyz].\n"
+ "Name is basename for the new raster3d maps");
param.dt = N_define_standard_option(N_OPT_CALC_TIME);
param.maxit = N_define_standard_option(N_OPT_MAX_ITERATIONS);
param.error = N_define_standard_option(N_OPT_ITERATION_ERROR);
param.solver = N_define_standard_option(N_OPT_SOLVER_SYMM);
- param.sor = N_define_standard_option(N_OPT_SOR_VALUE);
+ param.solver->options = "cg,pcg,cholesky";
param.mask = G_define_flag();
param.mask->key = 'm';
@@ -136,8 +139,7 @@
param.sparse = G_define_flag();
param.sparse->key = 's';
- param.sparse->description =
- _("Use a sparse linear equation system, only available with iterative solvers");
+ param.sparse->description = _("Use a sparse linear equation system, only available with iterative solvers");
}
/* ************************************************************************* */
@@ -146,19 +148,33 @@
int main(int argc, char *argv[])
{
struct GModule *module = NULL;
+
N_gwflow_data3d *data = NULL;
+
N_geom_data *geom = NULL;
+
N_les *les = NULL;
+
N_les_callback_3d *call = NULL;
+
G3D_Region region;
+
N_gradient_field_3d *field = NULL;
+
N_array_3d *xcomp = NULL;
+
N_array_3d *ycomp = NULL;
+
N_array_3d *zcomp = NULL;
- double error, sor;
+
+ double error;
+
int maxit;
+
const char *solver;
+
int x, y, z, stat;
+
char *buff = NULL;
/* Initialize GRASS */
@@ -166,8 +182,7 @@
module = G_define_module();
module->keywords = _("raster3d, voxel");
- module->description =
- _("Numerical calculation program for transient, confined groundwater flow in three dimensions");
+ module->description = _("Numerical calculation program for transient, confined groundwater flow in three dimensions");
/* Get parameters from user */
set_params();
@@ -180,14 +195,9 @@
sscanf(param.maxit->answer, "%i", &(maxit));
/*Set the calculation error break criteria */
sscanf(param.error->answer, "%lf", &(error));
- sscanf(param.sor->answer, "%lf", &(sor));
/*Set the solver */
solver = param.solver->answer;
- if (strcmp(solver, N_SOLVER_DIRECT_LU) == 0 && param.sparse->answer)
- G_fatal_error(_("The direct LU solver do not work with sparse matrices"));
- if (strcmp(solver, N_SOLVER_DIRECT_GAUSS) == 0 && param.sparse->answer)
- G_fatal_error(_("The direct Gauss solver do not work with sparse matrices"));
if (strcmp(solver, N_SOLVER_DIRECT_CHOLESKY) == 0 && param.sparse->answer)
G_fatal_error(_("The direct cholesky solver do not work with sparse matrices"));
@@ -264,32 +274,28 @@
(void *)data, call);
}
+ if (les && les->type == G_MATH_NORMAL_LES) {
+ if (strcmp(solver, N_SOLVER_ITERATIVE_CG) == 0)
+ G_math_solver_cg(les->A, les->x, les->b, les->rows, maxit, error);
- /*solve the equation system */
- if (strcmp(solver, N_SOLVER_ITERATIVE_JACOBI) == 0)
- N_solver_jacobi(les, maxit, sor, error);
+ if (strcmp(solver, N_SOLVER_ITERATIVE_PCG) == 0)
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, maxit, error,
+ N_DIAGONAL_PRECONDITION);
- if (strcmp(solver, N_SOLVER_ITERATIVE_SOR) == 0)
- N_solver_SOR(les, maxit, sor, error);
+ if (strcmp(solver, N_SOLVER_DIRECT_CHOLESKY) == 0)
+ G_math_solver_cholesky(les->A, les->x, les->b, les->rows,
+ les->rows);
+ }
+ else if (les && les->type == G_MATH_SPARSE_LES) {
+ if (strcmp(solver, N_SOLVER_ITERATIVE_CG) == 0)
+ G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows,
+ maxit, error);
- if (strcmp(solver, N_SOLVER_ITERATIVE_CG) == 0)
- N_solver_cg(les, maxit, error);
+ if (strcmp(solver, N_SOLVER_ITERATIVE_PCG) == 0)
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows,
+ maxit, error, N_DIAGONAL_PRECONDITION);
+ }
- if (strcmp(solver, N_SOLVER_ITERATIVE_PCG) == 0)
- N_solver_pcg(les, maxit, error, N_DIAGONAL_PRECONDITION);
-
- if (strcmp(solver, N_SOLVER_ITERATIVE_BICGSTAB) == 0)
- N_solver_bicgstab(les, maxit, error);
-
- if (strcmp(solver, N_SOLVER_DIRECT_LU) == 0)
- N_solver_lu(les);
-
- if (strcmp(solver, N_SOLVER_DIRECT_GAUSS) == 0)
- N_solver_gauss(les);
-
- if (strcmp(solver, N_SOLVER_DIRECT_CHOLESKY) == 0)
- N_solver_cholesky(les);
-
if (les == NULL)
G_fatal_error(_("Unable to create and solve the linear equation system"));
@@ -356,8 +362,11 @@
char *name)
{
void *map = NULL;
+
int changemask = 0;
+
int z, y, x, rows, cols, depths, count, stat;
+
double d1 = 0;
rows = region->rows;
More information about the grass-commit
mailing list