Added a new library to grass called ccmath (LGPL license) to replace the

NR algorithms of the gmath library.
Moved the linear equation solver code from gpde lib to gmath lib.
Added blas level 1, 2 and 3 algorithm in gmath lib.
Modified the gmath solver to use the grass blas implementation.
Added wrapper for ATLAS blas level 1 algorithms.
Updated the gpde library tests.
Added gmath library tests for the numerical part and ccmath wrapper.
Modified the groundwater flow modules (raster, raster3d)to use the gmath
solver.
Patched i.cca, i.pca and i.smap to use gmath vecotr and matrix functions
and the ccmath wrapper for eigen value computation.
Removed NR svd and eigen value code.



git-svn-id: https://svn.osgeo.org/grass/grass/trunk@39389 15284696-431f-4ddb-bdfa-cd5b030d7da7

Author: huhabla

imagery/i.cca/local_proto.h

@@ -3,22 +3,19 @@
 
 #include <grass/raster.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__ */

imagery/i.cca/main.c

@@ -53,26 +53,26 @@ int main(int argc, char *argv[])
     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;
@@ -129,9 +129,28 @@ int main(int argc, char *argv[])
 
     /* 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
@@ -154,14 +173,26 @@ int main(int argc, char *argv[])
 
     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++) {
@@ -205,6 +236,25 @@ int main(int argc, char *argv[])
 
     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);
 }

imagery/i.cca/matrix.c

@@ -4,25 +4,39 @@
 #include "local_proto.h"
 
 
-int product(double vector[MX], double factor, double matrix1[MX][MX],
+int print_matrix(double **matrix, int bands)
+{
+    int i, j;
+
+    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 +46,37 @@ int setdiag(double eigval[MX], int bands, double l[MX][MX])
 
 
 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;
+
+    tmp = G_alloc_matrix(bands, bands);
 
-    for (i = 1; i <= bands; i++)
+    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);
-    return 0;
-}
 
+    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);
 
-int solveq(double q[MX][MX], int bands, double w[MX][MX], double p[MX][MX])
-{
-    double tmp[MX][MX];
+    G_free_matrix(tmp);
 
-    matmul(tmp, w, p, bands);
-    matmul(q, tmp, w, bands);
     return 0;
 }
 
 
-int matmul(double res[MX][MX], double m1[MX][MX], double m2[MX][MX], int dim)
+int solveq(double **q, int bands, double **w, double **p)
 {
-    int i, j, k;
-    double sum;
-
-    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;
-	}
-    }
+    double **tmp;
+
+    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);
+
+    G_free_matrix(tmp);
 
     return 0;
 }
+

imagery/i.cca/stats.c

@@ -1,76 +1,68 @@
 #include <grass/gis.h>
-#include "local_proto.h"
+#include <grass/gmath.h>
 
+#include "local_proto.h"
 
 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;
 }
 
 
 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];
-
-    for (i = 0; i < MX; i++)
-	newvec[i] = 0.0;
+    double **tmp0, **tmp1, **tmp2;
+    double *newvec;
 
-    for (i = 1; i <= bands; i++)
-	for (j = 1; j <= bands; j++)
-	    tmp1[i][j] = tmp2[i][j] = 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 <= 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;
 }
+