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