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