[GRASS-SVN] r64265 - in grass/trunk: include/defs lib/gmath
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jan 20 14:31:51 PST 2015
Author: neteler
Date: 2015-01-20 14:31:51 -0800 (Tue, 20 Jan 2015)
New Revision: 64265
Modified:
grass/trunk/include/defs/la.h
grass/trunk/lib/gmath/la.c
Log:
libgmath: move new matrix/vector functions from i.spec.unmix Addon here
Modified: grass/trunk/include/defs/la.h
===================================================================
--- grass/trunk/include/defs/la.h 2015-01-20 21:03:39 UTC (rev 64264)
+++ grass/trunk/include/defs/la.h 2015-01-20 22:31:51 UTC (rev 64265)
@@ -27,6 +27,7 @@
vec_struct *G_matvect_get_column(mat_struct *, int);
vec_struct *G_matvect_get_row(mat_struct *, int);
+vec_struct *G_matvect_product(mat_struct *, vec_struct *, vec_struct *);
int G_matvect_extract_vector(mat_struct *, vtype, int);
int G_matvect_retrieve_matrix(vec_struct *);
@@ -38,6 +39,7 @@
double G_vector_norm_euclid(vec_struct *);
double G_vector_norm_maxval(vec_struct *, int);
vec_struct *G_vector_copy(const vec_struct *, int);
+vec_struct *G_vector_product (vec_struct *, vec_struct *, vec_struct *);
/* Matrix and vector routines corresponding to ?? */
@@ -47,5 +49,7 @@
int G_matrix_read(FILE *, mat_struct *);
int G_matrix_stdin(mat_struct *);
int G_matrix_eigen_sort(vec_struct *, mat_struct *);
+mat_struct *G_matrix_scalar_mul(double, mat_struct *, mat_struct *);
+mat_struct *G_matrix_resize(mat_struct *, int, int);
#endif /* GRASS_LADEFS_H */
Modified: grass/trunk/lib/gmath/la.c
===================================================================
--- grass/trunk/lib/gmath/la.c 2015-01-20 21:03:39 UTC (rev 64264)
+++ grass/trunk/lib/gmath/la.c 2015-01-20 22:31:51 UTC (rev 64265)
@@ -8,6 +8,7 @@
* 26th. Sep. 2000
* Last updated:
* 2006-11-23
+ * 2015-01-20
* This file is part of GRASS GIS. It is free software. You can
* redistribute it and/or modify it under the terms of
@@ -198,7 +199,48 @@
return G__matrix_add(mt1, mt2, 1, -1);
}
+/*!
+ * \fn mat_struct *G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
+ *
+ * \brief Calculates the scalar-matrix multiplication
+ *
+ * Calculates the scalar-matrix multiplication
+ *
+ * \param scalar
+ * \param A
+ * \return mat_struct
+ */
+mat_struct *G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
+{
+ int m, n, i, j;
+ int index = 0;
+
+ if (matrix == NULL) {
+ G_warning (_("Input matrix is uninitialized"));
+ return NULL;
+ }
+
+ if (out == NULL)
+ out = G_matrix_init(matrix->rows, matrix->cols, matrix->rows);
+
+ if (out->rows != matrix->rows || out->cols != matrix->cols)
+ out = G_matrix_resize(out, matrix->rows, matrix->cols);
+
+ m = matrix->rows;
+ n = matrix->cols;
+
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < n; j++) {
+ doublereal value = scalar * G_matrix_get_element(matrix, i, j);
+ G_matrix_set_element (out, i,j, value);
+ }
+ }
+
+ return (out);
+}
+
+
/*!
* \fn mat_struct *G_matrix_scale (mat_struct *mt1, const double c)
*
@@ -885,6 +927,58 @@
/*!
+ * \fn vec_struct *G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
+ *
+ * \brief Calculates the matrix-vector product
+ *
+ * Calculates the product of a matrix and a vector
+ *
+ * \param A
+ * \param b
+ * \return vec_struct
+ */
+
+vec_struct *G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
+{
+ unsigned int i, m, n, j;
+ register doublereal sum;
+
+/* G_message("A=%d,%d,%d", A->cols, A->rows, A->ldim); */
+/* G_message("B=%d,%d,%d", b->cols, b->rows, b->ldim); */
+ if (A->cols != b->cols) {
+ G_warning (_("Input matrix and vector have differing dimensions1"));
+
+ return NULL;
+
+ }
+ if (!out) {
+ G_warning (_("Output vector is uninitialized"));
+ return NULL;
+ }
+/* if (out->ldim != A->rows) {*/
+/* G_warning (_("Output vector has incorrect dimension"));*/
+/* exit(1);*/
+/* return NULL;*/
+/* }*/
+
+ m = A->rows;
+ n = A->cols;
+
+ for (i = 0; i < m; i++) {
+ sum = 0.0;
+ int width = A->rows;
+ for (j = 0; j < n; j++) {
+
+ sum +=G_matrix_get_element(A, i, j) * G_matrix_get_element(b, 0, j);
+ /*sum += A->vals[i + j * width] * b->vals[j];*/
+ out->vals[i] = sum;
+ }
+ }
+ return (out);
+}
+
+
+/*!
* \fn vec_struct *G_vector_init (int cells, int ldim, vtype vt)
*
* \brief Initialize a vector structure
@@ -1025,7 +1119,6 @@
return out;
}
-
/*!
* \fn int G_vector_set (vec_struct *A, int cells, int ldim, vtype vt, int vindx)
*
@@ -1247,7 +1340,81 @@
return result;
}
+/*!
+ * \fn vec_struct *G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out)
+ *
+ * \brief Calculates the vector product
+ *
+ * Calculates the vector product of two vectors
+ *
+ * \param v1
+ * \param v2
+ * \return vec_struct
+ */
+vec_struct *G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out)
+{
+ int idx1, idx2, idx0;
+ int i;
+
+ if (!out->is_init) {
+ G_warning (_("Output vector is uninitialized"));
+ return NULL;
+ }
+
+ if (v1->type != v2->type) {
+ G_warning (_("Vectors are not of the same type"));
+ return NULL;
+ }
+
+ if (v1->type != out->type) {
+ G_warning (_("Output vector is of incorrect type"));
+ return NULL;
+ }
+
+ if (v1->type == MATRIX_) {
+ G_warning (_("Matrices not allowed"));
+ return NULL;
+ }
+
+ if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
+ (v1->type == COLVEC_ && v1->rows != v2->rows))
+ {
+ G_warning (_("Vectors have differing dimensions"));
+ return NULL;
+ }
+
+ if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
+ (v1->type == COLVEC_ && v1->rows != out->rows))
+ {
+ G_warning (_("Output vector has incorrect dimension"));
+ return NULL;
+ }
+
+#if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS)
+ f77_dhad (v1->cols, 1.0, v1->vals, 1, v2->vals, 1, 0.0, out->vals, 1.0);
+#else
+ idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
+ idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
+ idx0 = (out->v_indx > 0) ? out->v_indx : 0;
+
+ if (v1->type == ROWVEC_) {
+ for (i = 0; i < v1->cols; i++)
+ G_matrix_set_element(out, idx0, i,
+ G_matrix_get_element(v1, idx1, i) *
+ G_matrix_get_element(v2, idx2, i));
+ } else {
+ for (i = 0; i < v1->rows; i++)
+ G_matrix_set_element(out, i, idx0,
+ G_matrix_get_element(v1, i, idx1) *
+ G_matrix_get_element(v2, i, idx2));
+ }
+#endif
+
+ return out;
+}
+
+
/*!
* \fn vec_struct *G_vector_copy (const vec_struct *vc1, int comp_flag)
*
@@ -1414,6 +1581,40 @@
/*!
+ * \fn mat_struct *G_matrix_resize(mat_struct *in, int rows, int cols)
+ *
+ * \brief Resizes a matrix
+ *
+ * Resizes a matrix
+ *
+ * \param A
+ * \param rows
+ * \param cols
+ * \return mat_struct
+ */
+
+mat_struct *G_matrix_resize(mat_struct *in, int rows, int cols)
+{
+ mat_struct *matrix;
+ matrix = G_matrix_init(rows, cols, rows);
+ int i, j, p, index = 0;
+ for (i = 0; i < rows; i++)
+ for (j = 0; j < cols; j++)
+ G_matrix_set_element(matrix, i, j, G_matrix_get_element(in, i, j));
+/* matrix->vals[index++] = in->vals[i + j * cols];*/
+
+ int old_size = in->rows * in->cols;
+ int new_size = rows * cols;
+
+ if (new_size > old_size)
+ for (p = old_size; p < new_size; p++)
+ G_matrix_set_element(matrix, i, j, 0.0);
+
+ return (matrix);
+}
+
+
+/*!
* \fn int G_matrix_read_stdin (mat_struct *out)
*
* \brief Read a matrix from standard input
More information about the grass-commit
mailing list