libgmath: move new matrix/vector functions from i.spec.unmix Addon here

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

include/defs/la.h

@@ -27,6 +27,7 @@ double G_matrix_get_element(mat_struct *, int, int);
 
 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 @@ int G_vector_set(vec_struct *, int, int, vtype, int);
 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 @@ double G_vector_norm1(vec_struct *);
 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 */

lib/gmath/la.c

@@ -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,6 +199,47 @@ mat_struct *G_matrix_subtract(mat_struct * mt1, mat_struct * mt2)
     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 @@ int G_matvect_retrieve_matrix(vec_struct * vc)
 
 
 /*!
+ * \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 @@ vec_struct *G_vector_sub(vec_struct * v1, vec_struct * v2, vec_struct * out)
     return out;
 }
 
-
 /*!
  * \fn int G_vector_set (vec_struct *A, int cells, int ldim, vtype vt, int vindx)
  *
@@ -1247,6 +1340,80 @@ double G_vector_norm1(vec_struct * vc)
     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 @@ int G_matrix_read(FILE * fp, mat_struct * out)
 
 
 /*!
+ * \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