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