GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
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 
31 /********
32  ******** only compile this LAPACK/BLAS wrapper file if g2c.h is present!
33  ********/
34 #include <grass/config.h>
35 
36 #if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS) && defined(HAVE_G2C_H)
37 
38 #include <grass/gis.h>
39 #include <grass/glocale.h>
40 #include <grass/la.h>
41 
42 
43 static int egcmp(const void *pa, const void *pb);
44 
45 
59 mat_struct *G_matrix_init(int rows, int cols, int ldim)
60 {
61  mat_struct *tmp_arry;
62 
63  if (rows < 1 || cols < 1 || ldim < rows) {
64  G_warning(_("Matrix dimensions out of range"));
65  return NULL;
66  }
67 
68  tmp_arry = (mat_struct *) G_malloc(sizeof(mat_struct));
69  tmp_arry->rows = rows;
70  tmp_arry->cols = cols;
71  tmp_arry->ldim = ldim;
72  tmp_arry->type = MATRIX_;
73  tmp_arry->v_indx = -1;
74 
75  tmp_arry->vals = (doublereal *) G_calloc(ldim * cols, sizeof(doublereal));
76  tmp_arry->is_init = 1;
77 
78  return tmp_arry;
79 }
80 
81 
91 int G_matrix_zero(mat_struct * A)
92 {
93  if (!A->vals)
94  return 0;
95 
96  memset(A->vals, 0, sizeof(A->vals));
97 
98  return 1;
99 }
100 
101 
117 int G_matrix_set(mat_struct * A, int rows, int cols, int ldim)
118 {
119  if (rows < 1 || cols < 1 || ldim < 0) {
120  G_warning(_("Matrix dimensions out of range"));
121  return -1;
122  }
123 
124  A->rows = rows;
125  A->cols = cols;
126  A->ldim = ldim;
127  A->type = MATRIX_;
128  A->v_indx = -1;
129 
130  A->vals = (doublereal *) G_calloc(ldim * cols, sizeof(doublereal));
131  A->is_init = 1;
132 
133  return 0;
134 }
135 
136 
148 mat_struct *G_matrix_copy(const mat_struct * A)
149 {
150  mat_struct *B;
151 
152  if (!A->is_init) {
153  G_warning(_("Matrix is not initialised fully."));
154  return NULL;
155  }
156 
157  if ((B = G_matrix_init(A->rows, A->cols, A->ldim)) == NULL) {
158  G_warning(_("Unable to allocate space for matrix copy"));
159  return NULL;
160  }
161 
162  memcpy(&B->vals[0], &A->vals[0], A->cols * A->ldim * sizeof(doublereal));
163 
164  return B;
165 }
166 
167 
181 mat_struct *G_matrix_add(mat_struct * mt1, mat_struct * mt2)
182 {
183  return G__matrix_add(mt1, mt2, 1, 1);
184 }
185 
186 
200 mat_struct *G_matrix_subtract(mat_struct * mt1, mat_struct * mt2)
201 {
202  return G__matrix_add(mt1, mt2, 1, -1);
203 }
204 
205 
219 mat_struct *G_matrix_scale(mat_struct * mt1, const double c)
220 {
221  return G__matrix_add(mt1, NULL, c, 0);
222 }
223 
224 
240 mat_struct *G__matrix_add(mat_struct * mt1, mat_struct * mt2, const double c1,
241  const double c2)
242 {
243  mat_struct *mt3;
244  int i, j; /* loop variables */
245 
246  if (c1 == 0) {
247  G_warning(_("First scalar multiplier must be non-zero"));
248  return NULL;
249  }
250 
251  if (c2 == 0) {
252  if (!mt1->is_init) {
253  G_warning(_("One or both input matrices uninitialised"));
254  return NULL;
255  }
256  }
257 
258  else {
259 
260  if (!((mt1->is_init) && (mt2->is_init))) {
261  G_warning(_("One or both input matrices uninitialised"));
262  return NULL;
263  }
264 
265  if (mt1->rows != mt2->rows || mt1->cols != mt2->cols) {
266  G_warning(_("Matrix order does not match"));
267  return NULL;
268  }
269  }
270 
271  if ((mt3 = G_matrix_init(mt1->rows, mt1->cols, mt1->ldim)) == NULL) {
272  G_warning(_("Unable to allocate space for matrix sum"));
273  return NULL;
274  }
275 
276  if (c2 == 0) {
277 
278  for (i = 0; i < mt3->rows; i++) {
279  for (j = 0; j < mt3->cols; j++) {
280  mt3->vals[i + mt3->ldim * j] =
281  c1 * mt1->vals[i + mt1->ldim * j];
282  }
283  }
284  }
285 
286  else {
287 
288  for (i = 0; i < mt3->rows; i++) {
289  for (j = 0; j < mt3->cols; j++) {
290  mt3->vals[i + mt3->ldim * j] =
291  c1 * mt1->vals[i + mt1->ldim * j] + c2 * mt2->vals[i +
292  mt2->
293  ldim *
294  j];
295  }
296  }
297  }
298 
299  return mt3;
300 }
301 
302 
303 #if defined(HAVE_LIBBLAS)
304 
318 mat_struct *G_matrix_product(mat_struct * mt1, mat_struct * mt2)
319 {
320  mat_struct *mt3;
321  doublereal unity = 1, zero = 0;
322  integer rows, cols, interdim, lda, ldb;
323  integer1 no_trans = 'n';
324 
325  if (!((mt1->is_init) || (mt2->is_init))) {
326  G_warning(_("One or both input matrices uninitialised"));
327  return NULL;
328  }
329 
330  if (mt1->cols != mt2->rows) {
331  G_warning(_("Matrix order does not match"));
332  return NULL;
333  }
334 
335  if ((mt3 = G_matrix_init(mt1->rows, mt2->cols, mt1->ldim)) == NULL) {
336  G_warning(_("Unable to allocate space for matrix product"));
337  return NULL;
338  }
339 
340  /* Call the driver */
341 
342  rows = (integer) mt1->rows;
343  interdim = (integer) mt1->cols;
344  cols = (integer) mt2->cols;
345 
346  lda = (integer) mt1->ldim;
347  ldb = (integer) mt2->ldim;
348 
349  f77_dgemm(&no_trans, &no_trans, &rows, &cols, &interdim, &unity,
350  mt1->vals, &lda, mt2->vals, &ldb, &zero, mt3->vals, &lda);
351 
352  return mt3;
353 }
354 
355 #else /* defined(HAVE_LIBBLAS) */
356 #warning G_matrix_product() not compiled; requires BLAS library
357 #endif /* defined(HAVE_LIBBLAS) */
358 
372 mat_struct *G_matrix_transpose(mat_struct * mt)
373 {
374  mat_struct *mt1;
375  int ldim, ldo;
376  doublereal *dbo, *dbt, *dbx, *dby;
377  int cnt, cnt2;
378 
379  /* Word align the workspace blocks */
380  if (mt->cols % 2 == 0)
381  ldim = mt->cols;
382  else
383  ldim = mt->cols + 1;
384 
385  mt1 = G_matrix_init(mt->cols, mt->rows, ldim);
386 
387  /* Set initial values for reading arrays */
388  dbo = &mt->vals[0];
389  dbt = &mt1->vals[0];
390  ldo = mt->ldim;
391 
392  for (cnt = 0; cnt < mt->cols; cnt++) {
393  dbx = dbo;
394  dby = dbt;
395 
396  for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
397  *dby = *dbx;
398  dby += ldim;
399  dbx++;
400  }
401 
402  *dby = *dbx;
403 
404  if (cnt < mt->cols - 1) {
405  dbo += ldo;
406  dbt++;
407  }
408  }
409 
410  return mt1;
411 }
412 
413 
414 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
415 
439 /*** NOT YET COMPLETE: only some solutions' options available ***/
440 
441 int
442 G_matrix_LU_solve(const mat_struct * mt1, mat_struct ** xmat0,
443  const mat_struct * bmat, mat_type mtype)
444 {
445  mat_struct *wmat, *xmat, *mtx;
446 
447  if (mt1->is_init == 0 || bmat->is_init == 0) {
448  G_warning(_("Input: one or both data matrices uninitialised"));
449  return -1;
450  }
451 
452  if (mt1->rows != mt1->cols || mt1->rows < 1) {
453  G_warning(_("Principal matrix is not properly dimensioned"));
454  return -1;
455  }
456 
457  if (bmat->cols < 1) {
458  G_warning(_("Input: you must have at least one array to solve"));
459  return -1;
460  }
461 
462  /* Now create solution matrix by copying the original coefficient matrix */
463  if ((xmat = G_matrix_copy(bmat)) == NULL) {
464  G_warning(_("Could not allocate space for solution matrix"));
465  return -1;
466  }
467 
468  /* Create working matrix for the coefficient array */
469  if ((mtx = G_matrix_copy(mt1)) == NULL) {
470  G_warning(_("Could not allocate space for working matrix"));
471  return -1;
472  }
473 
474  /* Copy the contents of the data matrix, to preserve the
475  original information
476  */
477  if ((wmat = G_matrix_copy(bmat)) == NULL) {
478  G_warning(_("Could not allocate space for working matrix"));
479  return -1;
480  }
481 
482  /* Now call appropriate LA driver to solve equations */
483  switch (mtype) {
484 
485  case NONSYM:
486  {
487  integer *perm, res_info;
488  integer num_eqns, nrhs, lda, ldb;
489 
490  perm = (integer *) G_malloc(wmat->rows);
491 
492  /* Set fields to pass to fortran routine */
493  num_eqns = (integer) mt1->rows;
494  nrhs = (integer) wmat->cols;
495  lda = (integer) mt1->ldim;
496  ldb = (integer) wmat->ldim;
497 
498  /* Call LA driver */
499  f77_dgesv(&num_eqns, &nrhs, mtx->vals, &lda, perm, wmat->vals,
500  &ldb, &res_info);
501 
502  /* Copy the results from the modified data matrix, taking account
503  of pivot permutations ???
504  */
505 
506  /*
507  for(indx1 = 0; indx1 < num_eqns; indx1++) {
508  iperm = perm[indx1];
509  ptin = &wmat->vals[0] + indx1;
510  ptout = &xmat->vals[0] + iperm;
511 
512  for(indx2 = 0; indx2 < nrhs - 1; indx2++) {
513  *ptout = *ptin;
514  ptin += wmat->ldim;
515  ptout += xmat->ldim;
516  }
517 
518  *ptout = *ptin;
519  }
520  */
521 
522  memcpy(xmat->vals, wmat->vals,
523  wmat->cols * wmat->ldim * sizeof(doublereal));
524 
525  /* Free temp arrays */
526  G_free(perm);
527  G_matrix_free(wmat);
528  G_matrix_free(mtx);
529 
530  if (res_info > 0) {
531  G_warning(_("Matrix (or submatrix is singular). Solution undetermined"));
532  return 1;
533  }
534  else if (res_info < 0) {
535  G_warning(_("Problem in LA routine."));
536  return -1;
537  }
538  break;
539  }
540  default:
541  {
542  G_warning(_("Procedure not yet available for selected matrix type"));
543  return -1;
544  }
545  } /* end switch */
546 
547  *xmat0 = xmat;
548 
549  return 0;
550 }
551 
552 #else /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
553 #warning G_matrix_LU_solve() not compiled; requires BLAS and LAPACK libraries
554 #endif /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
555 
556 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
557 
570 mat_struct *G_matrix_inverse(mat_struct * mt)
571 {
572  mat_struct *mt0, *res;
573  int i, j, k; /* loop */
574 
575  if (mt->rows != mt->cols) {
576  G_warning(_("Matrix is not square. Cannot determine inverse"));
577  return NULL;
578  }
579 
580  if ((mt0 = G_matrix_init(mt->rows, mt->rows, mt->ldim)) == NULL) {
581  G_warning(_("Unable to allocate space for matrix"));
582  return NULL;
583  }
584 
585  /* Set `B' matrix to unit matrix */
586  for (i = 0; i < mt0->rows - 1; i++) {
587  mt0->vals[i + i * mt0->ldim] = 1.0;
588 
589  for (j = i + 1; j < mt0->cols; j++) {
590  mt0->vals[i + j * mt0->ldim] = mt0->vals[j + i * mt0->ldim] = 0.0;
591  }
592  }
593 
594  mt0->vals[mt0->rows - 1 + (mt0->rows - 1) * mt0->ldim] = 1.0;
595 
596  /* Solve system */
597  if ((k = G_matrix_LU_solve(mt, &res, mt0, NONSYM)) == 1) {
598  G_warning(_("Matrix is singular"));
599  G_matrix_free(mt0);
600  return NULL;
601  }
602  else if (k < 0) {
603  G_warning(_("Problem in LA procedure."));
604  G_matrix_free(mt0);
605  return NULL;
606  }
607  else {
608  G_matrix_free(mt0);
609  return res;
610  }
611 }
612 
613 #else /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
614 #warning G_matrix_inverse() not compiled; requires BLAS and LAPACK libraries
615 #endif /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
616 
617 
629 void G_matrix_free(mat_struct * mt)
630 {
631  if (mt->is_init)
632  G_free(mt->vals);
633 
634  G_free(mt);
635 }
636 
637 
649 void G_matrix_print(mat_struct * mt)
650 {
651  int i, j;
652  char buf[64], numbuf[64];
653 
654  for (i = 0; i < mt->rows; i++) {
655  strcpy(buf, "");
656 
657  for (j = 0; j < mt->cols; j++) {
658 
659  sprintf(numbuf, "%14.6f", G_matrix_get_element(mt, i, j));
660  strcat(buf, numbuf);
661  if (j < mt->cols - 1)
662  strcat(buf, ", ");
663  }
664 
665  G_message("%s", buf);
666  }
667 
668  fprintf(stderr, "\n");
669 }
670 
671 
688 int G_matrix_set_element(mat_struct * mt, int rowval, int colval, double val)
689 {
690  if (!mt->is_init) {
691  G_warning(_("Element array has not been allocated"));
692  return -1;
693  }
694 
695  if (rowval >= mt->rows || colval >= mt->cols || rowval < 0 || colval < 0) {
696  G_warning(_("Specified element is outside array bounds"));
697  return -1;
698  }
699 
700  mt->vals[rowval + colval * mt->ldim] = (doublereal) val;
701 
702  return 0;
703 }
704 
705 
721 double G_matrix_get_element(mat_struct * mt, int rowval, int colval)
722 {
723  double val;
724 
725  /* Should do some checks, but this would require an error control
726  system: later? */
727 
728  return (val = (double)mt->vals[rowval + colval * mt->ldim]);
729 }
730 
731 
744 vec_struct *G_matvect_get_column(mat_struct * mt, int col)
745 {
746  int i; /* loop */
747  vec_struct *vc1;
748 
749  if (col < 0 || col >= mt->cols) {
750  G_warning(_("Specified matrix column index is outside range"));
751  return NULL;
752  }
753 
754  if (!mt->is_init) {
755  G_warning(_("Matrix is not initialised"));
756  return NULL;
757  }
758 
759  if ((vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) == NULL) {
760  G_warning(_("Could not allocate space for vector structure"));
761  return NULL;
762  }
763 
764  for (i = 0; i < mt->rows; i++)
765  G_matrix_set_element((mat_struct *) vc1, i, 0,
766  G_matrix_get_element(mt, i, col));
767 
768  return vc1;
769 }
770 
771 
785 vec_struct *G_matvect_get_row(mat_struct * mt, int row)
786 {
787  int i; /* loop */
788  vec_struct *vc1;
789 
790  if (row < 0 || row >= mt->cols) {
791  G_warning(_("Specified matrix row index is outside range"));
792  return NULL;
793  }
794 
795  if (!mt->is_init) {
796  G_warning(_("Matrix is not initialised"));
797  return NULL;
798  }
799 
800  if ((vc1 = G_vector_init(mt->cols, mt->ldim, RVEC)) == NULL) {
801  G_warning(_("Could not allocate space for vector structure"));
802  return NULL;
803  }
804 
805  for (i = 0; i < mt->cols; i++)
806  G_matrix_set_element((mat_struct *) vc1, 0, i,
807  G_matrix_get_element(mt, row, i));
808 
809  return vc1;
810 }
811 
812 
828 int G_matvect_extract_vector(mat_struct * mt, vtype vt, int indx)
829 {
830  if (vt == RVEC && indx >= mt->rows) {
831  G_warning(_("Specified row index is outside range"));
832  return -1;
833  }
834 
835  else if (vt == CVEC && indx >= mt->cols) {
836  G_warning(_("Specified column index is outside range"));
837  return -1;
838  }
839 
840  switch (vt) {
841 
842  case RVEC:
843  {
844  mt->type = ROWVEC_;
845  mt->v_indx = indx;
846  }
847 
848  case CVEC:
849  {
850  mt->type = COLVEC_;
851  mt->v_indx = indx;
852  }
853 
854  default:
855  {
856  G_warning(_("Unknown vector type."));
857  return -1;
858  }
859 
860  }
861 
862  return 0;
863 
864 }
865 
866 
878 int G_matvect_retrieve_matrix(vec_struct * vc)
879 {
880  /* We have to take the integrity of the vector structure
881  largely on trust
882  */
883 
884  vc->type = MATRIX_;
885  vc->v_indx = -1;
886 
887  return 0;
888 }
889 
890 
905 vec_struct *G_vector_init(int cells, int ldim, vtype vt)
906 {
907  vec_struct *tmp_arry;
908 
909  if ((cells < 1) || (vt == RVEC && ldim < 1)
910  || (vt == CVEC && ldim < cells) || ldim < 0) {
911  G_warning("Vector dimensions out of range.");
912  return NULL;
913  }
914 
915  tmp_arry = (vec_struct *) G_malloc(sizeof(vec_struct));
916 
917  if (vt == RVEC) {
918  tmp_arry->rows = 1;
919  tmp_arry->cols = cells;
920  tmp_arry->ldim = ldim;
921  tmp_arry->type = ROWVEC_;
922  }
923 
924  else if (vt == CVEC) {
925  tmp_arry->rows = cells;
926  tmp_arry->cols = 1;
927  tmp_arry->ldim = ldim;
928  tmp_arry->type = COLVEC_;
929  }
930 
931  tmp_arry->v_indx = 0;
932 
933  tmp_arry->vals = (doublereal *) G_calloc(ldim * tmp_arry->cols,
934  sizeof(doublereal));
935  tmp_arry->is_init = 1;
936 
937  return tmp_arry;
938 }
939 
940 
952 void G_vector_free(vec_struct * v)
953 {
954  if (v->is_init)
955  G_free(v->vals);
956 
957  G_free(v);
958 }
959 
960 
975 vec_struct *G_vector_sub(vec_struct * v1, vec_struct * v2, vec_struct * out)
976 {
977  int idx1, idx2, idx0;
978  int i;
979 
980  if (!out->is_init) {
981  G_warning(_("Output vector is uninitialized"));
982  return NULL;
983  }
984 
985  if (v1->type != v2->type) {
986  G_warning(_("Vectors are not of the same type"));
987  return NULL;
988  }
989 
990  if (v1->type != out->type) {
991  G_warning(_("Output vector is of incorrect type"));
992  return NULL;
993  }
994 
995  if (v1->type == MATRIX_) {
996  G_warning(_("Matrices not allowed"));
997  return NULL;
998  }
999 
1000  if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1001  (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1002  G_warning(_("Vectors have differing dimensions"));
1003  return NULL;
1004  }
1005 
1006  if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1007  (v1->type == COLVEC_ && v1->rows != out->rows)) {
1008  G_warning(_("Output vector has incorrect dimension"));
1009  return NULL;
1010  }
1011 
1012  idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
1013  idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
1014  idx0 = (out->v_indx > 0) ? out->v_indx : 0;
1015 
1016  if (v1->type == ROWVEC_) {
1017  for (i = 0; i < v1->cols; i++)
1018  G_matrix_set_element(out, idx0, i,
1019  G_matrix_get_element(v1, idx1, i) -
1020  G_matrix_get_element(v2, idx2, i));
1021  }
1022  else {
1023  for (i = 0; i < v1->rows; i++)
1024  G_matrix_set_element(out, i, idx0,
1025  G_matrix_get_element(v1, i, idx1) -
1026  G_matrix_get_element(v2, i, idx2));
1027  }
1028 
1029  return out;
1030 }
1031 
1032 
1050 int G_vector_set(vec_struct * A, int cells, int ldim, vtype vt, int vindx)
1051 {
1052  if ((cells < 1) || (vt == RVEC && ldim < 1)
1053  || (vt == CVEC && ldim < cells) || ldim < 0) {
1054  G_warning(_("Vector dimensions out of range"));
1055  return -1;
1056  }
1057 
1058  if ((vt == RVEC && vindx >= A->cols) || (vt == CVEC && vindx >= A->rows)) {
1059  G_warning(_("Row/column out of range"));
1060  return -1;
1061  }
1062 
1063  if (vt == RVEC) {
1064  A->rows = 1;
1065  A->cols = cells;
1066  A->ldim = ldim;
1067  A->type = ROWVEC_;
1068  }
1069  else {
1070  A->rows = cells;
1071  A->cols = 1;
1072  A->ldim = ldim;
1073  A->type = COLVEC_;
1074  }
1075 
1076  if (vindx < 0)
1077  A->v_indx = 0;
1078  else
1079  A->v_indx = vindx;
1080 
1081  A->vals = (doublereal *) G_calloc(ldim * A->cols, sizeof(doublereal));
1082  A->is_init = 1;
1083 
1084  return 0;
1085 
1086 }
1087 
1088 
1089 #if defined(HAVE_LIBBLAS)
1090 
1103 double G_vector_norm_euclid(vec_struct * vc)
1104 {
1105  integer incr, Nval;
1106  doublereal *startpt;
1107 
1108  if (!vc->is_init)
1109  G_fatal_error(_("Matrix is not initialised"));
1110 
1111  if (vc->type == ROWVEC_) {
1112  Nval = (integer) vc->cols;
1113  incr = (integer) vc->ldim;
1114  if (vc->v_indx < 0)
1115  startpt = vc->vals;
1116  else
1117  startpt = vc->vals + vc->v_indx;
1118  }
1119  else {
1120  Nval = (integer) vc->rows;
1121  incr = 1;
1122  if (vc->v_indx < 0)
1123  startpt = vc->vals;
1124  else
1125  startpt = vc->vals + vc->v_indx * vc->ldim;
1126  }
1127 
1128  /* Call the BLAS routine dnrm2_() */
1129  return (double)f77_dnrm2(&Nval, startpt, &incr);
1130 }
1131 
1132 #else /* defined(HAVE_LIBBLAS) */
1133 #warning G_vector_norm_euclid() not compiled; requires BLAS library
1134 #endif /* defined(HAVE_LIBBLAS) */
1135 
1136 
1154 double G_vector_norm_maxval(vec_struct * vc, int vflag)
1155 {
1156  doublereal xval, *startpt, *curpt;
1157  double cellval;
1158  int ncells, incr;
1159 
1160  if (!vc->is_init)
1161  G_fatal_error(_("Matrix is not initialised"));
1162 
1163  if (vc->type == ROWVEC_) {
1164  ncells = (integer) vc->cols;
1165  incr = (integer) vc->ldim;
1166  if (vc->v_indx < 0)
1167  startpt = vc->vals;
1168  else
1169  startpt = vc->vals + vc->v_indx;
1170  }
1171  else {
1172  ncells = (integer) vc->rows;
1173  incr = 1;
1174  if (vc->v_indx < 0)
1175  startpt = vc->vals;
1176  else
1177  startpt = vc->vals + vc->v_indx * vc->ldim;
1178  }
1179 
1180  xval = *startpt;
1181  curpt = startpt;
1182 
1183  while (ncells > 0) {
1184  if (curpt != startpt) {
1185  switch (vflag) {
1186 
1187  case MAX_POS:
1188  {
1189  if (*curpt > xval)
1190  xval = *curpt;
1191  break;
1192  }
1193 
1194  case MAX_NEG:
1195  {
1196  if (*curpt < xval)
1197  xval = *curpt;
1198  break;
1199  }
1200 
1201  case MAX_ABS:
1202  {
1203  cellval = (double)(*curpt);
1204  if (hypot(cellval, cellval) > (double)xval)
1205  xval = *curpt;
1206  }
1207  } /* switch */
1208  } /* if(curpt != startpt) */
1209 
1210  curpt += incr;
1211  ncells--;
1212  }
1213 
1214  return (double)xval;
1215 }
1216 
1217 
1229 double G_vector_norm1(vec_struct * vc)
1230 {
1231  double result = 0.0;
1232  int idx;
1233  int i;
1234 
1235  if (!vc->is_init) {
1236  G_warning(_("Matrix is not initialised"));
1237  return 0.0 / 0.0; /* NaN */
1238  }
1239 
1240  idx = (vc->v_indx > 0) ? vc->v_indx : 0;
1241 
1242  if (vc->type == ROWVEC_) {
1243  for (i = 0; i < vc->cols; i++)
1244  result += fabs(G_matrix_get_element(vc, idx, i));
1245  }
1246  else {
1247  for (i = 0; i < vc->rows; i++)
1248  result += fabs(G_matrix_get_element(vc, i, idx));
1249  }
1250 
1251  return result;
1252 }
1253 
1254 
1266 vec_struct *G_vector_copy(const vec_struct * vc1, int comp_flag)
1267 {
1268  vec_struct *tmp_arry;
1269  int incr1, incr2;
1270  doublereal *startpt1, *startpt2, *curpt1, *curpt2;
1271  int cnt;
1272 
1273  if (!vc1->is_init) {
1274  G_warning(_("Vector structure is not initialised"));
1275  return NULL;
1276  }
1277 
1278  tmp_arry = (vec_struct *) G_malloc(sizeof(vec_struct));
1279 
1280  if (comp_flag == DO_COMPACT) {
1281  if (vc1->type == ROWVEC_) {
1282  tmp_arry->rows = 1;
1283  tmp_arry->cols = vc1->cols;
1284  tmp_arry->ldim = 1;
1285  tmp_arry->type = ROWVEC_;
1286  tmp_arry->v_indx = 0;
1287  }
1288  else if (vc1->type == COLVEC_) {
1289  tmp_arry->rows = vc1->rows;
1290  tmp_arry->cols = 1;
1291  tmp_arry->ldim = vc1->ldim;
1292  tmp_arry->type = COLVEC_;
1293  tmp_arry->v_indx = 0;
1294  }
1295  else {
1296  G_warning("Type is not vector.");
1297  return NULL;
1298  }
1299  }
1300  else if (comp_flag == NO_COMPACT) {
1301  tmp_arry->v_indx = vc1->v_indx;
1302  tmp_arry->rows = vc1->rows;
1303  tmp_arry->cols = vc1->cols;
1304  tmp_arry->ldim = vc1->ldim;
1305  tmp_arry->type = vc1->type;
1306  }
1307  else {
1308  G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n");
1309  return NULL;
1310  }
1311 
1312  tmp_arry->vals = (doublereal *) G_calloc(tmp_arry->ldim * tmp_arry->cols,
1313  sizeof(doublereal));
1314  if (comp_flag == DO_COMPACT) {
1315  if (tmp_arry->type == ROWVEC_) {
1316  startpt1 = tmp_arry->vals;
1317  startpt2 = vc1->vals + vc1->v_indx;
1318  curpt1 = startpt1;
1319  curpt2 = startpt2;
1320  incr1 = 1;
1321  incr2 = vc1->ldim;
1322  cnt = vc1->cols;
1323  }
1324  else if (tmp_arry->type == COLVEC_) {
1325  startpt1 = tmp_arry->vals;
1326  startpt2 = vc1->vals + vc1->v_indx * vc1->ldim;
1327  curpt1 = startpt1;
1328  curpt2 = startpt2;
1329  incr1 = 1;
1330  incr2 = 1;
1331  cnt = vc1->rows;
1332  }
1333  else {
1334  G_warning("Structure type is not vector.");
1335  return NULL;
1336  }
1337  }
1338  else if (comp_flag == NO_COMPACT) {
1339  startpt1 = tmp_arry->vals;
1340  startpt2 = vc1->vals;
1341  curpt1 = startpt1;
1342  curpt2 = startpt2;
1343  incr1 = 1;
1344  incr2 = 1;
1345  cnt = vc1->ldim * vc1->cols;
1346  }
1347  else {
1348  G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n");
1349  return NULL;
1350  }
1351 
1352  while (cnt > 0) {
1353  memcpy(curpt1, curpt2, sizeof(doublereal));
1354  curpt1 += incr1;
1355  curpt2 += incr2;
1356  cnt--;
1357  }
1358 
1359  tmp_arry->is_init = 1;
1360 
1361  return tmp_arry;
1362 }
1363 
1364 
1379 int G_matrix_read(FILE * fp, mat_struct * out)
1380 {
1381  char buff[100];
1382  int rows, cols;
1383  int i, j, row;
1384  double val;
1385 
1386  /* skip comments */
1387  for (;;) {
1388  if (!G_getl(buff, sizeof(buff), fp))
1389  return -1;
1390  if (buff[0] != '#')
1391  break;
1392  }
1393 
1394  if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) {
1395  G_warning(_("Input format error"));
1396  return -1;
1397  }
1398 
1399  G_matrix_set(out, rows, cols, rows);
1400 
1401  for (i = 0; i < rows; i++) {
1402  if (fscanf(fp, "row%d:", &row) != 1 || row != i) {
1403  G_warning(_("Input format error"));
1404  return -1;
1405  }
1406  for (j = 0; j < cols; j++) {
1407  if (fscanf(fp, "%lf:", &val) != 1) {
1408  G_warning(_("Input format error"));
1409  return -1;
1410  }
1411 
1412  G_matrix_set_element(out, i, j, val);
1413  }
1414  }
1415 
1416  return 0;
1417 }
1418 
1419 
1433 int G_matrix_stdin(mat_struct * out)
1434 {
1435  return G_matrix_read(stdin, out);
1436 }
1437 
1438 
1451 int G_matrix_eigen_sort(vec_struct * d, mat_struct * m)
1452 {
1453  mat_struct tmp;
1454  int i, j;
1455  int idx;
1456 
1457  G_matrix_set(&tmp, m->rows + 1, m->cols, m->ldim + 1);
1458 
1459  idx = (d->v_indx > 0) ? d->v_indx : 0;
1460 
1461  /* concatenate (vertically) m and d into tmp */
1462  for (i = 0; i < m->cols; i++) {
1463  for (j = 0; j < m->rows; j++)
1464  G_matrix_set_element(&tmp, j + 1, i,
1465  G_matrix_get_element(m, j, i));
1466  if (d->type == ROWVEC_)
1467  G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, idx, i));
1468  else
1469  G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, i, idx));
1470  }
1471 
1472  /* sort the combined matrix */
1473  qsort(tmp.vals, tmp.cols, tmp.ldim * sizeof(doublereal), egcmp);
1474 
1475  /* split tmp into m and d */
1476  for (i = 0; i < m->cols; i++) {
1477  for (j = 0; j < m->rows; j++)
1478  G_matrix_set_element(m, j, i,
1479  G_matrix_get_element(&tmp, j + 1, i));
1480  if (d->type == ROWVEC_)
1481  G_matrix_set_element(d, idx, i, G_matrix_get_element(&tmp, 0, i));
1482  else
1483  G_matrix_set_element(d, i, idx, G_matrix_get_element(&tmp, 0, i));
1484  }
1485 
1486  G_free(tmp.vals);
1487 
1488  return 0;
1489 }
1490 
1491 
1492 static int egcmp(const void *pa, const void *pb)
1493 {
1494  double a = *(doublereal * const)pa;
1495  double b = *(doublereal * const)pb;
1496 
1497  if (a > b)
1498  return 1;
1499  if (a < b)
1500  return -1;
1501 
1502  return 0;
1503 }
1504 
1505 
1506 #endif /* HAVE_BLAS && HAVE_LAPACK && HAVE_G2C */