GRASS GIS 7 Programmer's Manual  7.0.5(2016)-r00000
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 
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 
88 int G_matrix_zero(mat_struct * A)
89 {
90  if (!A->vals)
91  return 0;
92 
93  memset(A->vals, 0, sizeof(A->vals));
94 
95  return 1;
96 }
97 
98 
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 
145 mat_struct *G_matrix_copy(const mat_struct * A)
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 
178 mat_struct *G_matrix_add(mat_struct * mt1, mat_struct * mt2)
179 {
180  return G__matrix_add(mt1, mt2, 1, 1);
181 }
182 
183 
197 mat_struct *G_matrix_subtract(mat_struct * mt1, mat_struct * mt2)
198 {
199  return G__matrix_add(mt1, mt2, 1, -1);
200 }
201 
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 
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 
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 
356 mat_struct *G_matrix_product(mat_struct * mt1, mat_struct * mt2)
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 
410 mat_struct *G_matrix_transpose(mat_struct * mt)
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 
477 /*** NOT YET COMPLETE: only some solutions' options available ***/
478 
479 int
480 G_matrix_LU_solve(const mat_struct * mt1, mat_struct ** xmat0,
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 
608 mat_struct *G_matrix_inverse(mat_struct * mt)
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 
667 void G_matrix_free(mat_struct * mt)
668 {
669  if (mt->is_init)
670  G_free(mt->vals);
671 
672  G_free(mt);
673 }
674 
675 
687 void G_matrix_print(mat_struct * mt)
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 
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 
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 
782 vec_struct *G_matvect_get_column(mat_struct * mt, int col)
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 
823 vec_struct *G_matvect_get_row(mat_struct * mt, int row)
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 
866 int G_matvect_extract_vector(mat_struct * mt, vtype vt, int indx)
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 
916 int G_matvect_retrieve_matrix(vec_struct * vc)
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 
941 vec_struct *G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
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 
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 
1042 void G_vector_free(vec_struct * v)
1043 {
1044  if (v->is_init)
1045  G_free(v->vals);
1046 
1047  G_free(v);
1048 }
1049 
1050 
1065 vec_struct *G_vector_sub(vec_struct * v1, vec_struct * v2, vec_struct * out)
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 
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 
1192 double G_vector_norm_euclid(vec_struct * vc)
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 
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 
1318 double G_vector_norm1(vec_struct * vc)
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 
1355 vec_struct *G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out)
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 
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 
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 
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 
1630 int G_matrix_stdin(mat_struct * out)
1631 {
1632  return G_matrix_read(stdin, out);
1633 }
1634 
1635 
1648 int G_matrix_eigen_sort(vec_struct * d, mat_struct * m)
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 */
#define NULL
Definition: ccmath.h:32
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
double b
if(!DBFLoadRecord(psDBF, hEntity)) return NULL
int G_getl(char *buf, int n, FILE *fd)
Gets a line of text from a file.
Definition: getl.c:31
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:89
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203