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