GRASS GIS 8 Programmer's Manual  8.5.0dev(2025)-3df7050671
la.c
Go to the documentation of this file.
1 /******************************************************************************
2  * la.c
3  * wrapper modules for linear algebra problems
4  * linking to BLAS / LAPACK (and others?)
5 
6  * @Copyright David D.Gray <ddgray@armadce.demon.co.uk>
7  * 26th. Sep. 2000
8  * Last updated:
9  * 2006-11-23
10  * 2015-01-20
11 
12  * This file is part of GRASS GIS. It is free software. You can
13  * redistribute it and/or modify it under the terms of
14  * the GNU General Public License as published by the Free Software
15  * Foundation; either version 2 of the License, or (at your option)
16  * any later version.
17 
18  * This program is distributed in the hope that it will be useful,
19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21  * GNU General Public License for more details.
22 
23  ******************************************************************************/
24 
25 #include <grass/config.h>
26 
27 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
28 
29 #include <math.h>
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <string.h>
33 
34 #if defined(_MSC_VER)
35 #include <complex.h>
36 #define LAPACK_COMPLEX_CUSTOM
37 #define lapack_complex_float _Fcomplex
38 #define lapack_complex_double _Dcomplex
39 #endif
40 
41 #include <lapacke.h>
42 #if defined(HAVE_CBLAS_ATLAS_H)
43 #include <cblas-atlas.h>
44 #else
45 #include <cblas.h>
46 #endif // HAVE_CBLAS_ATLAS_H
47 
48 #include <grass/gis.h>
49 #include <grass/glocale.h>
50 #include <grass/la.h>
51 
52 static int egcmp(const void *pa, const void *pb);
53 
54 /*!
55  * \fn mat_struct *G_matrix_init(int rows, int cols, int ldim)
56  *
57  * \brief Initialize a matrix structure
58  *
59  * Initialize a matrix structure
60  *
61  * \param rows
62  * \param cols
63  * \param ldim
64  * \return mat_struct
65  */
66 
67 mat_struct *G_matrix_init(int rows, int cols, int ldim)
68 {
69  mat_struct *tmp_arry;
70 
71  if (rows < 1 || cols < 1 || ldim < rows) {
72  G_warning(_("Matrix dimensions out of range"));
73  return NULL;
74  }
75 
76  tmp_arry = (mat_struct *)G_malloc(sizeof(mat_struct));
77  tmp_arry->rows = rows;
78  tmp_arry->cols = cols;
79  tmp_arry->ldim = ldim;
80  tmp_arry->type = MATRIX_;
81  tmp_arry->v_indx = -1;
82 
83  tmp_arry->vals = (double *)G_calloc(ldim * cols, sizeof(double));
84  tmp_arry->is_init = 1;
85 
86  return tmp_arry;
87 }
88 
89 /*!
90  * \fn int G_matrix_zero (mat_struct *A)
91  *
92  * \brief Clears (or resets) the matrix values to 0
93  *
94  * \param A
95  * \return 0 on error; 1 on success
96  */
97 
99 {
100  if (!A->vals)
101  return 0;
102 
103  memset(A->vals, 0, (A->ldim * A->cols) * sizeof(double));
104 
105  return 1;
106 }
107 
108 /*!
109  * \fn int G_matrix_set(mat_struct *A, int rows, int cols, int ldim)
110  *
111  * \brief Set parameters for an initialized matrix
112  *
113  * Set parameters for matrix <b>A</b> that is allocated,
114  * but not yet fully initialized. Is an alternative to G_matrix_init().
115  *
116  * \param A
117  * \param rows
118  * \param cols
119  * \param ldim
120  * \return int
121  */
122 
123 int G_matrix_set(mat_struct *A, int rows, int cols, int ldim)
124 {
125  if (rows < 1 || cols < 1 || ldim < 0) {
126  G_warning(_("Matrix dimensions out of range"));
127  return -1;
128  }
129 
130  A->rows = rows;
131  A->cols = cols;
132  A->ldim = ldim;
133  A->type = MATRIX_;
134  A->v_indx = -1;
135 
136  A->vals = (double *)G_calloc(ldim * cols, sizeof(double));
137  A->is_init = 1;
138 
139  return 0;
140 }
141 
142 /*!
143  * \fn mat_struct *G_matrix_copy (const mat_struct *A)
144  *
145  * \brief Copy a matrix
146  *
147  * Copy matrix <b>A</b> by exactly duplicating its contents.
148  *
149  * \param A
150  * \return mat_struct
151  */
152 
154 {
155  mat_struct *B;
156 
157  if (!A->is_init) {
158  G_warning(_("Matrix is not initialised fully."));
159  return NULL;
160  }
161 
162  if ((B = G_matrix_init(A->rows, A->cols, A->ldim)) == NULL) {
163  G_warning(_("Unable to allocate space for matrix copy"));
164  return NULL;
165  }
166 
167  memcpy(&B->vals[0], &A->vals[0],
168  (size_t)A->cols * A->ldim * sizeof(double));
169 
170  return B;
171 }
172 
173 /*!
174  * \fn mat_struct *G_matrix_add (mat_struct *mt1, mat_struct *mt2)
175  *
176  * \brief Adds two matricies
177  *
178  * Adds two matricies <b>mt1</b> and <b>mt2</b> and returns a
179  * resulting matrix. The return structure is automatically initialized.
180  *
181  * \param mt1
182  * \param mt2
183  * \return mat_struct
184  */
185 
187 {
188  return G__matrix_add(mt1, mt2, 1, 1);
189 }
190 
191 /*!
192  * \fn mat_struct *G_matrix_subtract (mat_struct *mt1, mat_struct *mt2)
193  *
194  * \brief Subtract two matricies
195  *
196  * Subtracts two matricies <b>mt1</b> and <b>mt2</b> and returns
197  * a resulting matrix. The return matrix is automatically initialized.
198  *
199  * \param mt1
200  * \param mt2
201  * \return mat_struct
202  */
203 
205 {
206  return G__matrix_add(mt1, mt2, 1, -1);
207 }
208 
209 /*!
210  * \fn mat_struct *G_matrix_scalar_mul(double scalar, mat_struct *matrix,
211  * mat_struct *out)
212  *
213  * \brief Calculates the scalar-matrix multiplication
214  *
215  * Calculates the scalar-matrix multiplication
216  *
217  * \param scalar
218  * \param A
219  * \return mat_struct
220  */
221 
222 mat_struct *G_matrix_scalar_mul(double scalar, mat_struct *matrix,
223  mat_struct *out)
224 {
225  int m, n, i, j;
226 
227  if (matrix == NULL) {
228  G_warning(_("Input matrix is uninitialized"));
229  return NULL;
230  }
231 
232  if (out == NULL)
233  out = G_matrix_init(matrix->rows, matrix->cols, matrix->rows);
234 
235  if (out->rows != matrix->rows || out->cols != matrix->cols)
236  out = G_matrix_resize(out, matrix->rows, matrix->cols);
237 
238  m = matrix->rows;
239  n = matrix->cols;
240 
241  for (i = 0; i < m; i++) {
242  for (j = 0; j < n; j++) {
243  double value = scalar * G_matrix_get_element(matrix, i, j);
244 
245  G_matrix_set_element(out, i, j, value);
246  }
247  }
248 
249  return (out);
250 }
251 
252 /*!
253  * \fn mat_struct *G_matrix_scale (mat_struct *mt1, const double c)
254  *
255  * \brief Scale a matrix by a scalar value
256  *
257  * Scales matrix <b>mt1</b> by scalar value <b>c</b>. The
258  * resulting matrix is automatically initialized.
259  *
260  * \param mt1
261  * \param c
262  * \return mat_struct
263  */
264 
265 mat_struct *G_matrix_scale(mat_struct *mt1, const double c)
266 {
267  return G__matrix_add(mt1, NULL, c, 0);
268 }
269 
270 /*!
271  * \fn mat_struct *G__matrix_add (mat_struct *mt1, mat_struct *mt2, const double
272  * c1, const double c2)
273  *
274  * \brief General add/subtract/scalar multiply routine
275  *
276  * General add/subtract/scalar multiply routine. <b>c2</b> may be
277  * zero, but <b>c1</b> must be non-zero.
278  *
279  * \param mt1
280  * \param mt2
281  * \param c1
282  * \param c2
283  * \return mat_struct
284  */
285 
286 mat_struct *G__matrix_add(mat_struct *mt1, mat_struct *mt2, const double c1,
287  const double c2)
288 {
289  mat_struct *mt3;
290  int i, j; /* loop variables */
291 
292  if (c1 == 0) {
293  G_warning(_("First scalar multiplier must be non-zero"));
294  return NULL;
295  }
296 
297  if (c2 == 0) {
298  if (!mt1->is_init) {
299  G_warning(_("One or both input matrices uninitialised"));
300  return NULL;
301  }
302  }
303 
304  else {
305 
306  if (!((mt1->is_init) && (mt2->is_init))) {
307  G_warning(_("One or both input matrices uninitialised"));
308  return NULL;
309  }
310 
311  if (mt1->rows != mt2->rows || mt1->cols != mt2->cols) {
312  G_warning(_("Matrix order does not match"));
313  return NULL;
314  }
315  }
316 
317  if ((mt3 = G_matrix_init(mt1->rows, mt1->cols, mt1->ldim)) == NULL) {
318  G_warning(_("Unable to allocate space for matrix sum"));
319  return NULL;
320  }
321 
322  if (c2 == 0) {
323 
324  for (i = 0; i < mt3->rows; i++) {
325  for (j = 0; j < mt3->cols; j++) {
326  mt3->vals[i + mt3->ldim * j] =
327  c1 * mt1->vals[i + mt1->ldim * j];
328  }
329  }
330  }
331 
332  else {
333 
334  for (i = 0; i < mt3->rows; i++) {
335  for (j = 0; j < mt3->cols; j++) {
336  mt3->vals[i + mt3->ldim * j] =
337  c1 * mt1->vals[i + mt1->ldim * j] +
338  c2 * mt2->vals[i + mt2->ldim * j];
339  }
340  }
341  }
342 
343  return mt3;
344 }
345 
346 /*!
347  * \fn mat_struct *G_matrix_product (mat_struct *mt1, mat_struct *mt2)
348  *
349  * \brief Returns product of two matricies
350  *
351  * Returns a matrix with the product of matrix <b>mt1</b> and
352  * <b>mt2</b>. The return matrix is automatically initialized.
353  *
354  * \param mt1
355  * \param mt2
356  * \return mat_struct
357  */
358 
360 {
361  mat_struct *mt3;
362  double unity = 1., zero = 0.;
363  int rows, cols, interdim, lda, ldb;
364 
365  if (!((mt1->is_init) || (mt2->is_init))) {
366  G_warning(_("One or both input matrices uninitialised"));
367  return NULL;
368  }
369 
370  if (mt1->cols != mt2->rows) {
371  G_warning(_("Matrix order does not match"));
372  return NULL;
373  }
374 
375  if ((mt3 = G_matrix_init(mt1->rows, mt2->cols, mt1->ldim)) == NULL) {
376  G_warning(_("Unable to allocate space for matrix product"));
377  return NULL;
378  }
379 
380  /* Call the driver */
381 
382  rows = (int)mt1->rows;
383  interdim = (int)mt1->cols;
384  cols = (int)mt2->cols;
385 
386  lda = (int)mt1->ldim;
387  ldb = (int)mt2->ldim;
388 
389  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, rows, cols, interdim,
390  unity, mt1->vals, lda, mt2->vals, ldb, zero, mt3->vals, lda);
391 
392  return mt3;
393 }
394 
395 /*!
396  * \fn mat_struct *G_matrix_transpose (mat_struct *mt)
397  *
398  * \brief Transpose a matrix
399  *
400  * Transpose matrix <b>m1</b> by creating a new one and
401  * populating with transposed elements. The return matrix is
402  * automatically initialized.
403  *
404  * \param mt
405  * \return mat_struct
406  */
407 
409 {
410  mat_struct *mt1;
411  int ldim, ldo;
412  double *dbo, *dbt, *dbx, *dby;
413  int cnt, cnt2;
414 
415  /* Word align the workspace blocks */
416  if (mt->cols % 2 == 0)
417  ldim = mt->cols;
418  else
419  ldim = mt->cols + 1;
420 
421  mt1 = G_matrix_init(mt->cols, mt->rows, ldim);
422 
423  /* Set initial values for reading arrays */
424  dbo = &mt->vals[0];
425  dbt = &mt1->vals[0];
426  ldo = mt->ldim;
427 
428  for (cnt = 0; cnt < mt->cols; cnt++) {
429  dbx = dbo;
430  dby = dbt;
431 
432  for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
433  *dby = *dbx;
434  dby += ldim;
435  dbx++;
436  }
437 
438  *dby = *dbx;
439 
440  if (cnt < mt->cols - 1) {
441  dbo += ldo;
442  dbt++;
443  }
444  }
445 
446  return mt1;
447 }
448 
449 /*!
450  * \fn int G_matrix_LU_solve (const mat_struct *mt1, mat_struct **xmat0,
451  * const mat_struct *bmat, mat_type mtype)
452  *
453  * \brief Solve a general system A.X = B
454  *
455  * Solve a general system A.X = B, where A is a NxN matrix, X and B are
456  * NxC matrices, and we are to solve for C arrays in X given B. Uses LU
457  * decomposition.<br>
458  * Links to LAPACK function dgesv_() and similar to perform the core routine.
459  * (By default solves for a general non-symmetric matrix.)<br>
460  * mtype is a flag to indicate what kind of matrix (real/complex, Hermitian,
461  * symmetric, general etc.) is used (NONSYM, SYM, HERMITIAN).<br>
462  * <b>Warning:</b> NOT YET COMPLETE: only some solutions' options
463  * available. Now, only general real matrix is supported.
464  *
465  * \param mt1
466  * \param xmat0
467  * \param bmat
468  * \param mtype
469  * \return int
470  */
471 
472 /*** NOT YET COMPLETE: only some solutions' options available ***/
473 
474 int G_matrix_LU_solve(const mat_struct *mt1, mat_struct **xmat0,
475  const mat_struct *bmat, mat_type mtype)
476 {
477  mat_struct *wmat, *xmat, *mtx;
478 
479  if (mt1->is_init == 0 || bmat->is_init == 0) {
480  G_warning(_("Input: one or both data matrices uninitialised"));
481  return -1;
482  }
483 
484  if (mt1->rows != mt1->cols || mt1->rows < 1) {
485  G_warning(_("Principal matrix is not properly dimensioned"));
486  return -1;
487  }
488 
489  if (bmat->cols < 1) {
490  G_warning(_("Input: you must have at least one array to solve"));
491  return -1;
492  }
493 
494  /* Now create solution matrix by copying the original coefficient matrix */
495  if ((xmat = G_matrix_copy(bmat)) == NULL) {
496  G_warning(_("Could not allocate space for solution matrix"));
497  return -1;
498  }
499 
500  /* Create working matrix for the coefficient array */
501  if ((mtx = G_matrix_copy(mt1)) == NULL) {
502  G_warning(_("Could not allocate space for working matrix"));
503  return -1;
504  }
505 
506  /* Copy the contents of the data matrix, to preserve the
507  original information
508  */
509  if ((wmat = G_matrix_copy(bmat)) == NULL) {
510  G_warning(_("Could not allocate space for working matrix"));
511  return -1;
512  }
513 
514  /* Now call appropriate LA driver to solve equations */
515  switch (mtype) {
516 
517  case NONSYM: {
518  int *perm, res_info;
519  int num_eqns, nrhs, lda, ldb;
520 
521  perm = (int *)G_malloc(wmat->rows * sizeof(int));
522 
523  /* Set fields to pass to fortran routine */
524  num_eqns = (int)mt1->rows;
525  nrhs = (int)wmat->cols;
526  lda = (int)mt1->ldim;
527  ldb = (int)wmat->ldim;
528 
529  /* Call LA driver */
530  res_info = LAPACKE_dgesv(LAPACK_COL_MAJOR, num_eqns, nrhs, mtx->vals,
531  lda, perm, wmat->vals, ldb);
532 
533  /* Copy the results from the modified data matrix, taking account
534  of pivot permutations ???
535  */
536 
537  /*
538  for(indx1 = 0; indx1 < num_eqns; indx1++) {
539  iperm = perm[indx1];
540  ptin = &wmat->vals[0] + indx1;
541  ptout = &xmat->vals[0] + iperm;
542 
543  for(indx2 = 0; indx2 < nrhs - 1; indx2++) {
544  *ptout = *ptin;
545  ptin += wmat->ldim;
546  ptout += xmat->ldim;
547  }
548 
549  *ptout = *ptin;
550  }
551  */
552 
553  memcpy(xmat->vals, wmat->vals,
554  (size_t)wmat->cols * wmat->ldim * sizeof(double));
555 
556  /* Free temp arrays */
557  G_free(perm);
558  G_matrix_free(wmat);
559  G_matrix_free(mtx);
560 
561  if (res_info > 0) {
562  G_warning(
563  _("Matrix (or submatrix is singular). Solution undetermined"));
564  return 1;
565  }
566  else if (res_info < 0) {
567  G_warning(_("Problem in LA routine."));
568  return -1;
569  }
570  break;
571  }
572  default: {
573  G_warning(_("Procedure not yet available for selected matrix type"));
574  return -1;
575  }
576  } /* end switch */
577 
578  *xmat0 = xmat;
579 
580  return 0;
581 }
582 
583 /*!
584  * \fn mat_struct *G_matrix_inverse (mat_struct *mt)
585  *
586  * \brief Returns the matrix inverse
587  *
588  * Calls G_matrix_LU_solve() to obtain matrix inverse using LU
589  * decomposition. Returns NULL on failure.
590  *
591  * \param mt
592  * \return mat_struct
593  */
594 
596 {
597  mat_struct *mt0, *res;
598  int i, j, k; /* loop */
599 
600  if (mt->rows != mt->cols) {
601  G_warning(_("Matrix is not square. Cannot determine inverse"));
602  return NULL;
603  }
604 
605  if ((mt0 = G_matrix_init(mt->rows, mt->rows, mt->ldim)) == NULL) {
606  G_warning(_("Unable to allocate space for matrix"));
607  return NULL;
608  }
609 
610  /* Set `B' matrix to unit matrix */
611  for (i = 0; i < mt0->rows - 1; i++) {
612  mt0->vals[i + i * mt0->ldim] = 1.0;
613 
614  for (j = i + 1; j < mt0->cols; j++) {
615  mt0->vals[i + j * mt0->ldim] = mt0->vals[j + i * mt0->ldim] = 0.0;
616  }
617  }
618 
619  mt0->vals[mt0->rows - 1 + (mt0->rows - 1) * mt0->ldim] = 1.0;
620 
621  /* Solve system */
622  if ((k = G_matrix_LU_solve(mt, &res, mt0, NONSYM)) == 1) {
623  G_warning(_("Matrix is singular"));
624  G_matrix_free(mt0);
625  return NULL;
626  }
627  else if (k < 0) {
628  G_warning(_("Problem in LA procedure."));
629  G_matrix_free(mt0);
630  return NULL;
631  }
632  else {
633  G_matrix_free(mt0);
634  return res;
635  }
636 }
637 
638 /*!
639  * \fn void G_matrix_free (mat_struct *mt)
640  *
641  * \brief Free up allocated matrix
642  *
643  * Free up allocated matrix.
644  *
645  * \param mt
646  * \return void
647  */
648 
650 {
651  if (mt->is_init)
652  G_free(mt->vals);
653 
654  G_free(mt);
655 }
656 
657 /*!
658  * \fn void G_matrix_print (mat_struct *mt)
659  *
660  * \brief Print out a matrix
661  *
662  * Print out a representation of the matrix to standard output.
663  *
664  * \param mt
665  * \return void
666  */
667 
669 {
670  int i, j;
671  char buf[2048], numbuf[64];
672 
673  for (i = 0; i < mt->rows; i++) {
674  G_strlcpy(buf, "", sizeof(buf));
675 
676  for (j = 0; j < mt->cols; j++) {
677  snprintf(numbuf, sizeof(numbuf), "%14.6f",
678  G_matrix_get_element(mt, i, j));
679  strcat(buf, numbuf);
680  if (j < mt->cols - 1)
681  strcat(buf, ", ");
682  }
683 
684  G_message("%s", buf);
685  }
686 
687  fprintf(stderr, "\n");
688 }
689 
690 /*!
691  * \fn int G_matrix_set_element (mat_struct *mt, int rowval, int colval, double
692  * val)
693  *
694  * \brief Set the value of the (i, j)th element
695  *
696  * Set the value of the (i, j)th element to a double value. Index values
697  * are C-like ie. zero-based. The row number is given first as is
698  * conventional. Returns -1 if the accessed cell is outside the bounds.
699  *
700  * \param mt
701  * \param rowval
702  * \param colval
703  * \param val
704  * \return int
705  */
706 
707 int G_matrix_set_element(mat_struct *mt, int rowval, int colval, double val)
708 {
709  if (!mt->is_init) {
710  G_warning(_("Element array has not been allocated"));
711  return -1;
712  }
713 
714  if (rowval >= mt->rows || colval >= mt->cols || rowval < 0 || colval < 0) {
715  G_warning(_("Specified element is outside array bounds"));
716  return -1;
717  }
718 
719  mt->vals[rowval + colval * mt->ldim] = (double)val;
720 
721  return 0;
722 }
723 
724 /*!
725  * \fn double G_matrix_get_element (mat_struct *mt, int rowval, int colval)
726  *
727  * \brief Retrieve value of the (i,j)th element
728  *
729  * Retrieve the value of the (i, j)th element to a double value. Index
730  * values are C-like ie. zero-based.
731  * <b>Note:</b> Does currently not set an error flag for bounds checking.
732  *
733  * \param mt
734  * \param rowval
735  * \param colval
736  * \return double
737  */
738 
739 double G_matrix_get_element(mat_struct *mt, int rowval, int colval)
740 {
741  double val;
742 
743  /* Should do some checks, but this would require an error control
744  system: later? */
745 
746  return (val = (double)mt->vals[rowval + colval * mt->ldim]);
747 }
748 
749 /*!
750  * \fn vec_struct *G_matvect_get_column (mat_struct *mt, int col)
751  *
752  * \brief Retrieve a column of the matrix to a vector structure
753  *
754  * Retrieve a column of matrix <b>mt</b> to a returning vector structure
755  *
756  * \param mt
757  * \param col
758  * \return vec_struct
759  */
760 
762 {
763  int i; /* loop */
764  vec_struct *vc1;
765 
766  if (col < 0 || col >= mt->cols) {
767  G_warning(_("Specified matrix column index is outside range"));
768  return NULL;
769  }
770 
771  if (!mt->is_init) {
772  G_warning(_("Matrix is not initialised"));
773  return NULL;
774  }
775 
776  if ((vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) == NULL) {
777  G_warning(_("Could not allocate space for vector structure"));
778  return NULL;
779  }
780 
781  for (i = 0; i < mt->rows; i++)
782  G_matrix_set_element((mat_struct *)vc1, i, 0,
783  G_matrix_get_element(mt, i, col));
784 
785  return vc1;
786 }
787 
788 /*!
789  * \fn vec_struct *G_matvect_get_row (mat_struct *mt, int row)
790  *
791  * \brief Retrieve a row of the matrix to a vector structure
792  *
793  * Retrieves a row from matrix <b>mt</b> and returns it in a vector
794  * structure.
795  *
796  * \param mt
797  * \param row
798  * \return vec_struct
799  */
800 
802 {
803  int i; /* loop */
804  vec_struct *vc1;
805 
806  if (row < 0 || row >= mt->cols) {
807  G_warning(_("Specified matrix row index is outside range"));
808  return NULL;
809  }
810 
811  if (!mt->is_init) {
812  G_warning(_("Matrix is not initialised"));
813  return NULL;
814  }
815 
816  if ((vc1 = G_vector_init(mt->cols, mt->ldim, RVEC)) == NULL) {
817  G_warning(_("Could not allocate space for vector structure"));
818  return NULL;
819  }
820 
821  for (i = 0; i < mt->cols; i++)
822  G_matrix_set_element((mat_struct *)vc1, 0, i,
823  G_matrix_get_element(mt, row, i));
824 
825  return vc1;
826 }
827 
828 /*!
829  * \fn int G_matvect_extract_vector (mat_struct *mt, vtype vt, int indx)
830  *
831  * \brief Convert matrix to vector
832  *
833  * Convert the matrix <b>mt</b> to a vector structure. The vtype,
834  * <b>vt</b>, is RVEC or CVEC which specifies a row vector or column
835  * vector. The index, <b>indx</b>, indicates the row/column number (zero based).
836  *
837  * \param mt
838  * \param vt
839  * \param indx
840  * \return int
841  */
842 
844 {
845  if (vt == RVEC && indx >= mt->rows) {
846  G_warning(_("Specified row index is outside range"));
847  return -1;
848  }
849 
850  else if (vt == CVEC && indx >= mt->cols) {
851  G_warning(_("Specified column index is outside range"));
852  return -1;
853  }
854 
855  switch (vt) {
856 
857  case RVEC: {
858  mt->type = ROWVEC_;
859  mt->v_indx = indx;
860  break;
861  }
862 
863  case CVEC: {
864  mt->type = COLVEC_;
865  mt->v_indx = indx;
866  break;
867  }
868 
869  default: {
870  G_warning(_("Unknown vector type."));
871  return -1;
872  }
873  }
874 
875  return 0;
876 }
877 
878 /*!
879  * \fn int G_matvect_retrieve_matrix (vec_struct *vc)
880  *
881  * \brief Revert a vector to matrix
882  *
883  * Revert vector <b>vc</b> to a matrix.
884  *
885  * \param vc
886  * \return int
887  */
888 
890 {
891  /* We have to take the integrity of the vector structure
892  largely on trust
893  */
894 
895  vc->type = MATRIX_;
896  vc->v_indx = -1;
897 
898  return 0;
899 }
900 
901 /*!
902  * \fn vec_struct *G_matvect_product(mat_struct *A, vec_struct *b, vec_struct
903  * *out)
904  *
905  * \brief Calculates the matrix-vector product
906  *
907  * Calculates the product of a matrix and a vector
908  *
909  * \param A
910  * \param b
911  * \return vec_struct
912  */
913 
915 {
916  unsigned int i, m, n, j;
917  register double sum;
918 
919  /* G_message("A=%d,%d,%d", A->cols, A->rows, A->ldim); */
920  /* G_message("B=%d,%d,%d", b->cols, b->rows, b->ldim); */
921  if (A->cols != b->cols) {
922  G_warning(_("Input matrix and vector have differing dimensions1"));
923 
924  return NULL;
925  }
926  if (!out) {
927  G_warning(_("Output vector is uninitialized"));
928  return NULL;
929  }
930  /* if (out->ldim != A->rows) { */
931  /* G_warning (_("Output vector has incorrect dimension")); */
932  /* exit(1); */
933  /* return NULL; */
934  /* } */
935 
936  m = A->rows;
937  n = A->cols;
938 
939  for (i = 0; i < m; i++) {
940  sum = 0.0;
941  /* int width = A->rows; */
942 
943  for (j = 0; j < n; j++) {
944 
945  sum +=
946  G_matrix_get_element(A, i, j) * G_matrix_get_element(b, 0, j);
947  /*sum += A->vals[i + j * width] * b->vals[j]; */
948  out->vals[i] = sum;
949  }
950  }
951  return (out);
952 }
953 
954 /*!
955  * \fn vec_struct *G_vector_init (int cells, int ldim, vtype vt)
956  *
957  * \brief Initialize a vector structure
958  *
959  * Returns an initialized vector structure with <b>cell</b> cells,
960  * of dimension <b>ldim</b>, and of type <b>vt</b>.
961  *
962  * \param cells
963  * \param ldim
964  * \param vt
965  * \return vec_struct
966  */
967 
968 vec_struct *G_vector_init(int cells, int ldim, vtype vt)
969 {
970  vec_struct *tmp_arry;
971 
972  if ((cells < 1) || (vt == RVEC && ldim < 1) ||
973  (vt == CVEC && ldim < cells) || ldim < 0) {
974  G_warning("Vector dimensions out of range.");
975  return NULL;
976  }
977 
978  tmp_arry = (vec_struct *)G_malloc(sizeof(vec_struct));
979 
980  if (vt == RVEC) {
981  tmp_arry->rows = 1;
982  tmp_arry->cols = cells;
983  tmp_arry->ldim = ldim;
984  tmp_arry->type = ROWVEC_;
985  }
986 
987  else if (vt == CVEC) {
988  tmp_arry->rows = cells;
989  tmp_arry->cols = 1;
990  tmp_arry->ldim = ldim;
991  tmp_arry->type = COLVEC_;
992  }
993 
994  tmp_arry->v_indx = 0;
995 
996  tmp_arry->vals = (double *)G_calloc(ldim * tmp_arry->cols, sizeof(double));
997  tmp_arry->is_init = 1;
998 
999  return tmp_arry;
1000 }
1001 
1002 /*!
1003  * \fn void G_vector_free (vec_struct *v)
1004  *
1005  * \brief Free an allocated vector structure
1006  *
1007  * Free an allocated vector structure.
1008  *
1009  * \param v
1010  * \return void
1011  */
1012 
1014 {
1015  if (v->is_init)
1016  G_free(v->vals);
1017 
1018  G_free(v);
1019 }
1020 
1021 /*!
1022  * \fn vec_struct *G_vector_sub (vec_struct *v1, vec_struct *v2, vec_struct
1023  * *out)
1024  *
1025  * \brief Subtract two vectors
1026  *
1027  * Subtracts two vectors, <b>v1</b> and <b>v2</b>, and returns and
1028  * populates vector <b>out</b>.
1029  *
1030  * \param v1
1031  * \param v2
1032  * \param out
1033  * \return vec_struct
1034  */
1035 
1037 {
1038  int idx1, idx2, idx0;
1039  int i;
1040 
1041  if (!out->is_init) {
1042  G_warning(_("Output vector is uninitialized"));
1043  return NULL;
1044  }
1045 
1046  if (v1->type != v2->type) {
1047  G_warning(_("Vectors are not of the same type"));
1048  return NULL;
1049  }
1050 
1051  if (v1->type != out->type) {
1052  G_warning(_("Output vector is of incorrect type"));
1053  return NULL;
1054  }
1055 
1056  if (v1->type == MATRIX_) {
1057  G_warning(_("Matrices not allowed"));
1058  return NULL;
1059  }
1060 
1061  if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1062  (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1063  G_warning(_("Vectors have differing dimensions"));
1064  return NULL;
1065  }
1066 
1067  if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1068  (v1->type == COLVEC_ && v1->rows != out->rows)) {
1069  G_warning(_("Output vector has incorrect dimension"));
1070  return NULL;
1071  }
1072 
1073  idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
1074  idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
1075  idx0 = (out->v_indx > 0) ? out->v_indx : 0;
1076 
1077  if (v1->type == ROWVEC_) {
1078  for (i = 0; i < v1->cols; i++)
1079  G_matrix_set_element(out, idx0, i,
1080  G_matrix_get_element(v1, idx1, i) -
1081  G_matrix_get_element(v2, idx2, i));
1082  }
1083  else {
1084  for (i = 0; i < v1->rows; i++)
1085  G_matrix_set_element(out, i, idx0,
1086  G_matrix_get_element(v1, i, idx1) -
1087  G_matrix_get_element(v2, i, idx2));
1088  }
1089 
1090  return out;
1091 }
1092 
1093 /*!
1094  * \fn int G_vector_set (vec_struct *A, int cells, int ldim, vtype vt, int
1095  * vindx)
1096  *
1097  * \brief Set parameters for vector structure
1098  *
1099  * Set parameters for a vector structure that is
1100  * allocated but not yet initialised fully. The vtype is RVEC or
1101  * CVEC which specifies a row vector or column vector.
1102  *
1103  * \param A
1104  * \param cells
1105  * \param ldim
1106  * \param vt
1107  * \param vindx
1108  * \return int
1109  */
1110 
1111 int G_vector_set(vec_struct *A, int cells, int ldim, vtype vt, int vindx)
1112 {
1113  if ((cells < 1) || (vt == RVEC && ldim < 1) ||
1114  (vt == CVEC && ldim < cells) || ldim < 0) {
1115  G_warning(_("Vector dimensions out of range"));
1116  return -1;
1117  }
1118 
1119  if ((vt == RVEC && vindx >= A->cols) || (vt == CVEC && vindx >= A->rows)) {
1120  G_warning(_("Row/column out of range"));
1121  return -1;
1122  }
1123 
1124  if (vt == RVEC) {
1125  A->rows = 1;
1126  A->cols = cells;
1127  A->ldim = ldim;
1128  A->type = ROWVEC_;
1129  }
1130  else {
1131  A->rows = cells;
1132  A->cols = 1;
1133  A->ldim = ldim;
1134  A->type = COLVEC_;
1135  }
1136 
1137  if (vindx < 0)
1138  A->v_indx = 0;
1139  else
1140  A->v_indx = vindx;
1141 
1142  A->vals = (double *)G_calloc(ldim * A->cols, sizeof(double));
1143  A->is_init = 1;
1144 
1145  return 0;
1146 }
1147 
1148 /*!
1149  * \fn double G_vector_norm_euclid (vec_struct *vc)
1150  *
1151  * \brief Calculates euclidean norm
1152  *
1153  * Calculates the euclidean norm of a row or column vector, using BLAS
1154  * routine dnrm2_().
1155  *
1156  * \param vc
1157  * \return double
1158  */
1159 
1161 {
1162  int incr, Nval;
1163  double *startpt;
1164 
1165  if (!vc->is_init)
1166  G_fatal_error(_("Matrix is not initialised"));
1167 
1168  if (vc->type == ROWVEC_) {
1169  Nval = (int)vc->cols;
1170  incr = (int)vc->ldim;
1171  if (vc->v_indx < 0)
1172  startpt = vc->vals;
1173  else
1174  startpt = vc->vals + vc->v_indx;
1175  }
1176  else {
1177  Nval = (int)vc->rows;
1178  incr = 1;
1179  if (vc->v_indx < 0)
1180  startpt = vc->vals;
1181  else
1182  startpt = vc->vals + vc->v_indx * vc->ldim;
1183  }
1184 
1185  /* Call the BLAS routine dnrm2_() */
1186  return cblas_dnrm2(Nval, startpt, incr);
1187 }
1188 
1189 /*!
1190  * \fn double G_vector_norm_maxval (vec_struct *vc, int vflag)
1191  *
1192  * \brief Calculates maximum value
1193  *
1194  * Calculates the maximum value of a row or column vector.
1195  * The vflag setting defines which value to be calculated:
1196  * vflag:
1197  * 1 Indicates maximum value<br>
1198  * -1 Indicates minimum value<br>
1199  * 0 Indicates absolute value [???]
1200  *
1201  * \param vc
1202  * \param vflag
1203  * \return double
1204  */
1205 
1206 double G_vector_norm_maxval(vec_struct *vc, int vflag)
1207 {
1208  double xval, *startpt, *curpt;
1209  double cellval;
1210  int ncells, incr;
1211 
1212  if (!vc->is_init)
1213  G_fatal_error(_("Matrix is not initialised"));
1214 
1215  if (vc->type == ROWVEC_) {
1216  ncells = (int)vc->cols;
1217  incr = (int)vc->ldim;
1218  if (vc->v_indx < 0)
1219  startpt = vc->vals;
1220  else
1221  startpt = vc->vals + vc->v_indx;
1222  }
1223  else {
1224  ncells = (int)vc->rows;
1225  incr = 1;
1226  if (vc->v_indx < 0)
1227  startpt = vc->vals;
1228  else
1229  startpt = vc->vals + vc->v_indx * vc->ldim;
1230  }
1231 
1232  xval = *startpt;
1233  curpt = startpt;
1234 
1235  while (ncells > 0) {
1236  if (curpt != startpt) {
1237  switch (vflag) {
1238 
1239  case MAX_POS: {
1240  if (*curpt > xval)
1241  xval = *curpt;
1242  break;
1243  }
1244 
1245  case MAX_NEG: {
1246  if (*curpt < xval)
1247  xval = *curpt;
1248  break;
1249  }
1250 
1251  case MAX_ABS: {
1252  cellval = (double)(*curpt);
1253  if (hypot(cellval, cellval) > (double)xval)
1254  xval = *curpt;
1255  }
1256  } /* switch */
1257  } /* if(curpt != startpt) */
1258 
1259  curpt += incr;
1260  ncells--;
1261  }
1262 
1263  return (double)xval;
1264 }
1265 
1266 /*!
1267  * \fn double G_vector_norm1 (vec_struct *vc)
1268  *
1269  * \brief Calculates the 1-norm of a vector
1270  *
1271  * Calculates the 1-norm of a vector
1272  *
1273  * \param vc
1274  * \return double
1275  */
1276 
1278 {
1279  double result = 0.0;
1280  int idx;
1281  int i;
1282 
1283  if (!vc->is_init) {
1284  G_warning(_("Matrix is not initialised"));
1285  return NAN;
1286  }
1287 
1288  idx = (vc->v_indx > 0) ? vc->v_indx : 0;
1289 
1290  if (vc->type == ROWVEC_) {
1291  for (i = 0; i < vc->cols; i++)
1292  result += fabs(G_matrix_get_element(vc, idx, i));
1293  }
1294  else {
1295  for (i = 0; i < vc->rows; i++)
1296  result += fabs(G_matrix_get_element(vc, i, idx));
1297  }
1298 
1299  return result;
1300 }
1301 
1302 /*!
1303  * \fn vec_struct *G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct
1304  * *out)
1305  *
1306  * \brief Calculates the vector product
1307  *
1308  * Calculates the vector product of two vectors
1309  *
1310  * \param v1
1311  * \param v2
1312  * \return vec_struct
1313  */
1314 
1316 {
1317  int idx1, idx2, idx0;
1318  int i;
1319 
1320  if (!out->is_init) {
1321  G_warning(_("Output vector is uninitialized"));
1322  return NULL;
1323  }
1324 
1325  if (v1->type != v2->type) {
1326  G_warning(_("Vectors are not of the same type"));
1327  return NULL;
1328  }
1329 
1330  if (v1->type != out->type) {
1331  G_warning(_("Output vector is not the same type as others"));
1332  return NULL;
1333  }
1334 
1335  if (v1->type == MATRIX_) {
1336  G_warning(_("Matrices not allowed"));
1337  return NULL;
1338  }
1339 
1340  if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1341  (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1342  G_warning(_("Vectors have differing dimensions"));
1343  return NULL;
1344  }
1345 
1346  if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1347  (v1->type == COLVEC_ && v1->rows != out->rows)) {
1348  G_warning(_("Output vector has incorrect dimension"));
1349  return NULL;
1350  }
1351 
1352  idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
1353  idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
1354  idx0 = (out->v_indx > 0) ? out->v_indx : 0;
1355 
1356  if (v1->type == ROWVEC_) {
1357  for (i = 0; i < v1->cols; i++)
1358  G_matrix_set_element(out, idx0, i,
1359  G_matrix_get_element(v1, idx1, i) *
1360  G_matrix_get_element(v2, idx2, i));
1361  }
1362  else {
1363  for (i = 0; i < v1->rows; i++)
1364  G_matrix_set_element(out, i, idx0,
1365  G_matrix_get_element(v1, i, idx1) *
1366  G_matrix_get_element(v2, i, idx2));
1367  }
1368 
1369  return out;
1370 }
1371 
1372 /*!
1373  * \fn vec_struct *G_vector_copy (const vec_struct *vc1, int comp_flag)
1374  *
1375  * \brief Returns a vector copied from <b>vc1</b>. Underlying structure
1376  * is preserved unless DO_COMPACT flag.
1377  *
1378  * \param vc1
1379  * \param comp_flag
1380  * \return vec_struct
1381  */
1382 
1383 vec_struct *G_vector_copy(const vec_struct *vc1, int comp_flag)
1384 {
1385  vec_struct *tmp_arry;
1386  int incr1, incr2;
1387  double *startpt1, *startpt2, *curpt1, *curpt2;
1388  int cnt;
1389 
1390  if (!vc1->is_init) {
1391  G_warning(_("Vector structure is not initialised"));
1392  return NULL;
1393  }
1394 
1395  tmp_arry = (vec_struct *)G_malloc(sizeof(vec_struct));
1396 
1397  if (comp_flag == DO_COMPACT) {
1398  if (vc1->type == ROWVEC_) {
1399  tmp_arry->rows = 1;
1400  tmp_arry->cols = vc1->cols;
1401  tmp_arry->ldim = 1;
1402  tmp_arry->type = ROWVEC_;
1403  tmp_arry->v_indx = 0;
1404  }
1405  else if (vc1->type == COLVEC_) {
1406  tmp_arry->rows = vc1->rows;
1407  tmp_arry->cols = 1;
1408  tmp_arry->ldim = vc1->ldim;
1409  tmp_arry->type = COLVEC_;
1410  tmp_arry->v_indx = 0;
1411  }
1412  else {
1413  G_warning("Type is not vector.");
1414  return NULL;
1415  }
1416  }
1417  else if (comp_flag == NO_COMPACT) {
1418  tmp_arry->v_indx = vc1->v_indx;
1419  tmp_arry->rows = vc1->rows;
1420  tmp_arry->cols = vc1->cols;
1421  tmp_arry->ldim = vc1->ldim;
1422  tmp_arry->type = vc1->type;
1423  }
1424  else {
1425  G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n");
1426  return NULL;
1427  }
1428 
1429  tmp_arry->vals =
1430  (double *)G_calloc(tmp_arry->ldim * tmp_arry->cols, sizeof(double));
1431  if (comp_flag == DO_COMPACT) {
1432  if (tmp_arry->type == ROWVEC_) {
1433  startpt1 = tmp_arry->vals;
1434  startpt2 = vc1->vals + vc1->v_indx;
1435  curpt1 = startpt1;
1436  curpt2 = startpt2;
1437  incr1 = 1;
1438  incr2 = vc1->ldim;
1439  cnt = vc1->cols;
1440  }
1441  else if (tmp_arry->type == COLVEC_) {
1442  startpt1 = tmp_arry->vals;
1443  startpt2 = vc1->vals + vc1->v_indx * vc1->ldim;
1444  curpt1 = startpt1;
1445  curpt2 = startpt2;
1446  incr1 = 1;
1447  incr2 = 1;
1448  cnt = vc1->rows;
1449  }
1450  else {
1451  G_warning("Structure type is not vector.");
1452  return NULL;
1453  }
1454  }
1455  else if (comp_flag == NO_COMPACT) {
1456  startpt1 = tmp_arry->vals;
1457  startpt2 = vc1->vals;
1458  curpt1 = startpt1;
1459  curpt2 = startpt2;
1460  incr1 = 1;
1461  incr2 = 1;
1462  cnt = vc1->ldim * vc1->cols;
1463  }
1464  else {
1465  G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n");
1466  return NULL;
1467  }
1468 
1469  while (cnt > 0) {
1470  memcpy(curpt1, curpt2, sizeof(double));
1471  curpt1 += incr1;
1472  curpt2 += incr2;
1473  cnt--;
1474  }
1475 
1476  tmp_arry->is_init = 1;
1477 
1478  return tmp_arry;
1479 }
1480 
1481 /*!
1482  * \fn int G_matrix_read (FILE *fp, mat_struct *out)
1483  *
1484  * \brief Read a matrix from a file stream
1485  *
1486  * Populates matrix structure <b>out</b> with matrix read from file
1487  * stream <b>fp</b>. Matrix <b>out</b> is automatically initialized.
1488  * Returns -1 on error and 0 on success.
1489  *
1490  * \param fp
1491  * \param out
1492  * \return int
1493  */
1494 
1495 int G_matrix_read(FILE *fp, mat_struct *out)
1496 {
1497  char buff[100];
1498  int rows, cols;
1499  int i, j, row;
1500  double val;
1501 
1502  /* skip comments */
1503  for (;;) {
1504  if (!G_getl(buff, sizeof(buff), fp))
1505  return -1;
1506  if (buff[0] != '#')
1507  break;
1508  }
1509 
1510  if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) {
1511  G_warning(_("Input format error"));
1512  return -1;
1513  }
1514 
1515  G_matrix_set(out, rows, cols, rows);
1516 
1517  for (i = 0; i < rows; i++) {
1518  if (fscanf(fp, "row%d:", &row) != 1 || row != i) {
1519  G_warning(_("Input format error"));
1520  return -1;
1521  }
1522  for (j = 0; j < cols; j++) {
1523  if (fscanf(fp, "%lf:", &val) != 1) {
1524  G_warning(_("Input format error"));
1525  return -1;
1526  }
1527 
1528  G_matrix_set_element(out, i, j, val);
1529  }
1530  }
1531 
1532  return 0;
1533 }
1534 
1535 /*!
1536  * \fn mat_struct *G_matrix_resize(mat_struct *in, int rows, int cols)
1537  *
1538  * \brief Resizes a matrix
1539  *
1540  * Resizes a matrix
1541  *
1542  * \param A
1543  * \param rows
1544  * \param cols
1545  * \return mat_struct
1546  */
1547 
1548 mat_struct *G_matrix_resize(mat_struct *in, int rows, int cols)
1549 {
1550  mat_struct *matrix;
1551 
1552  matrix = G_matrix_init(rows, cols, rows);
1553  int i, j, p /*, index = 0 */;
1554 
1555  for (i = 0; i < rows; i++)
1556  for (j = 0; j < cols; j++)
1557  G_matrix_set_element(matrix, i, j, G_matrix_get_element(in, i, j));
1558  /* matrix->vals[index++] = in->vals[i + j * cols]; */
1559 
1560  int old_size = in->rows * in->cols;
1561  int new_size = rows * cols;
1562 
1563  if (new_size > old_size)
1564  for (p = old_size; p < new_size; p++)
1565  G_matrix_set_element(matrix, i, j, 0.0);
1566 
1567  return (matrix);
1568 }
1569 
1570 /*!
1571  * \fn int G_matrix_read_stdin (mat_struct *out)
1572  *
1573  * \brief Read a matrix from standard input
1574  *
1575  * Populates matrix <b>out</b> with matrix read from stdin. Matrix
1576  * <b>out</b> is automatically initialized. Returns -1 on failure or 0
1577  * on success.
1578  *
1579  * \param out
1580  * \return int
1581  */
1582 
1584 {
1585  return G_matrix_read(stdin, out);
1586 }
1587 
1588 /*!
1589  * \fn int G_matrix_eigen_sort (vec_struct *d, mat_struct *m)
1590  *
1591  * \brief Sort eigenvectors according to eigenvalues
1592  *
1593  * Sort eigenvectors according to eigenvalues. Returns 0.
1594  *
1595  * \param d
1596  * \param m
1597  * \return int
1598  */
1599 
1601 {
1602  mat_struct tmp;
1603  int i, j;
1604  int idx;
1605 
1606  G_matrix_set(&tmp, m->rows + 1, m->cols, m->ldim + 1);
1607 
1608  idx = (d->v_indx > 0) ? d->v_indx : 0;
1609 
1610  /* concatenate (vertically) m and d into tmp */
1611  for (i = 0; i < m->cols; i++) {
1612  for (j = 0; j < m->rows; j++)
1613  G_matrix_set_element(&tmp, j + 1, i, G_matrix_get_element(m, j, i));
1614  if (d->type == ROWVEC_)
1615  G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, idx, i));
1616  else
1617  G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, i, idx));
1618  }
1619 
1620  /* sort the combined matrix */
1621  qsort(tmp.vals, tmp.cols, tmp.ldim * sizeof(double), egcmp);
1622 
1623  /* split tmp into m and d */
1624  for (i = 0; i < m->cols; i++) {
1625  for (j = 0; j < m->rows; j++)
1626  G_matrix_set_element(m, j, i, G_matrix_get_element(&tmp, j + 1, i));
1627  if (d->type == ROWVEC_)
1628  G_matrix_set_element(d, idx, i, G_matrix_get_element(&tmp, 0, i));
1629  else
1630  G_matrix_set_element(d, i, idx, G_matrix_get_element(&tmp, 0, i));
1631  }
1632 
1633  G_free(tmp.vals);
1634 
1635  return 0;
1636 }
1637 
1638 static int egcmp(const void *pa, const void *pb)
1639 {
1640  double a = *(double *const)pa;
1641  double b = *(double *const)pb;
1642 
1643  if (a > b)
1644  return 1;
1645  if (a < b)
1646  return -1;
1647 
1648  return 0;
1649 }
1650 
1651 #endif // HAVE_LIBLAPACK HAVE_LIBBLAS
1652 
#define NULL
Definition: ccmath.h:32
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_calloc(m, n)
Definition: defs/gis.h:95
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
#define G_malloc(n)
Definition: defs/gis.h:94
void G_message(const char *,...) __attribute__((format(printf
int G_getl(char *, int, FILE *)
Gets a line of text from a file.
Definition: getl.c:33
size_t G_strlcpy(char *, const char *, size_t)
Safe string copy function.
Definition: strlcpy.c:52
#define _(str)
Definition: glocale.h:10
vec_struct * G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
Calculates the matrix-vector product.
Definition: la.c:914
double G_vector_norm_maxval(vec_struct *vc, int vflag)
Calculates maximum value.
Definition: la.c:1206
vec_struct * G_matvect_get_column(mat_struct *mt, int col)
Retrieve a column of the matrix to a vector structure.
Definition: la.c:761
vec_struct * G_vector_sub(vec_struct *v1, vec_struct *v2, vec_struct *out)
Subtract two vectors.
Definition: la.c:1036
int G_matvect_retrieve_matrix(vec_struct *vc)
Revert a vector to matrix.
Definition: la.c:889
double G_vector_norm_euclid(vec_struct *vc)
Calculates euclidean norm.
Definition: la.c:1160
mat_struct * G_matrix_subtract(mat_struct *mt1, mat_struct *mt2)
Subtract two matricies.
Definition: la.c:204
mat_struct * G_matrix_add(mat_struct *mt1, mat_struct *mt2)
Adds two matricies.
Definition: la.c:186
double G_matrix_get_element(mat_struct *mt, int rowval, int colval)
Retrieve value of the (i,j)th element.
Definition: la.c:739
mat_struct * G__matrix_add(mat_struct *mt1, mat_struct *mt2, const double c1, const double c2)
Definition: la.c:286
mat_struct * G_matrix_inverse(mat_struct *mt)
Returns the matrix inverse.
Definition: la.c:595
int G_matrix_stdin(mat_struct *out)
Definition: la.c:1583
void G_matrix_print(mat_struct *mt)
Print out a matrix.
Definition: la.c:668
mat_struct * G_matrix_copy(const mat_struct *A)
Copy a matrix.
Definition: la.c:153
mat_struct * G_matrix_resize(mat_struct *in, int rows, int cols)
Resizes a matrix.
Definition: la.c:1548
int G_matvect_extract_vector(mat_struct *mt, vtype vt, int indx)
Convert matrix to vector.
Definition: la.c:843
vec_struct * G_vector_init(int cells, int ldim, vtype vt)
Initialize a vector structure.
Definition: la.c:968
mat_struct * G_matrix_product(mat_struct *mt1, mat_struct *mt2)
Returns product of two matricies.
Definition: la.c:359
int G_matrix_LU_solve(const mat_struct *mt1, mat_struct **xmat0, const mat_struct *bmat, mat_type mtype)
Solve a general system A.X = B.
Definition: la.c:474
vec_struct * G_vector_product(vec_struct *v1, vec_struct *v2, vec_struct *out)
Calculates the vector product.
Definition: la.c:1315
mat_struct * G_matrix_init(int rows, int cols, int ldim)
Initialize a matrix structure.
Definition: la.c:67
int suppress_empty_translation_unit_compiler_warning
Definition: la.c:1653
double G_vector_norm1(vec_struct *vc)
Calculates the 1-norm of a vector.
Definition: la.c:1277
int G_matrix_read(FILE *fp, mat_struct *out)
Read a matrix from a file stream.
Definition: la.c:1495
int G_matrix_eigen_sort(vec_struct *d, mat_struct *m)
Sort eigenvectors according to eigenvalues.
Definition: la.c:1600
mat_struct * G_matrix_scale(mat_struct *mt1, const double c)
Scale a matrix by a scalar value.
Definition: la.c:265
void G_vector_free(vec_struct *v)
Free an allocated vector structure.
Definition: la.c:1013
vec_struct * G_matvect_get_row(mat_struct *mt, int row)
Retrieve a row of the matrix to a vector structure.
Definition: la.c:801
void G_matrix_free(mat_struct *mt)
Free up allocated matrix.
Definition: la.c:649
vec_struct * G_vector_copy(const vec_struct *vc1, int comp_flag)
Returns a vector copied from vc1. Underlying structure is preserved unless DO_COMPACT flag.
Definition: la.c:1383
int G_vector_set(vec_struct *A, int cells, int ldim, vtype vt, int vindx)
Definition: la.c:1111
int G_matrix_set_element(mat_struct *mt, int rowval, int colval, double val)
Definition: la.c:707
mat_struct * G_matrix_transpose(mat_struct *mt)
Transpose a matrix.
Definition: la.c:408
mat_struct * G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
Calculates the scalar-matrix multiplication.
Definition: la.c:222
int G_matrix_set(mat_struct *A, int rows, int cols, int ldim)
Set parameters for an initialized matrix.
Definition: la.c:123
int G_matrix_zero(mat_struct *A)
Clears (or resets) the matrix values to 0.
Definition: la.c:98
Wrapper headers for BLAS/LAPACK.
#define DO_COMPACT
Definition: la.h:35
vtype
Definition: la.h:44
@ CVEC
Definition: la.h:44
@ RVEC
Definition: la.h:44
#define MAX_NEG
Definition: la.h:32
#define MAX_ABS
Definition: la.h:33
mat_type
Definition: la.h:42
@ NONSYM
Definition: la.h:42
#define NO_COMPACT
Definition: la.h:36
#define MAX_POS
Definition: la.h:31
@ COLVEC_
Definition: la.h:43
@ MATRIX_
Definition: la.h:43
@ ROWVEC_
Definition: la.h:43
double b
Definition: r_raster.c:39
if(!(yy_init))
Definition: sqlp.yy.c:775
Definition: la.h:53
int v_indx
Definition: la.h:56
mat_spec type
Definition: la.h:55
int rows
Definition: la.h:59
int is_init
Definition: la.h:63
int ldim
Definition: la.h:60
double * vals
Definition: la.h:62
int cols
Definition: la.h:59