31 #include <grass/config.h>
33 #if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS)
35 #include <grass/gis.h>
36 #include <grass/glocale.h>
40 static int egcmp(
const void *pa,
const void *pb);
56 mat_struct *G_matrix_init(
int rows,
int cols,
int ldim)
60 if (rows < 1 || cols < 1 || ldim < rows) {
61 G_warning(_(
"Matrix dimensions out of range"));
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;
72 tmp_arry->vals = (doublereal *) G_calloc(ldim * cols,
sizeof(doublereal));
73 tmp_arry->is_init = 1;
88 int G_matrix_zero(mat_struct * A)
93 memset(A->vals, 0,
sizeof(A->vals));
114 int G_matrix_set(mat_struct * A,
int rows,
int cols,
int ldim)
116 if (rows < 1 || cols < 1 || ldim < 0) {
117 G_warning(_(
"Matrix dimensions out of range"));
127 A->vals = (doublereal *) G_calloc(ldim * cols,
sizeof(doublereal));
145 mat_struct *G_matrix_copy(
const mat_struct * A)
150 G_warning(_(
"Matrix is not initialised fully."));
154 if ((B = G_matrix_init(A->rows, A->cols, A->ldim)) ==
NULL) {
155 G_warning(_(
"Unable to allocate space for matrix copy"));
159 memcpy(&B->vals[0], &A->vals[0], A->cols * A->ldim *
sizeof(doublereal));
178 mat_struct *G_matrix_add(mat_struct * mt1, mat_struct * mt2)
180 return G__matrix_add(mt1, mt2, 1, 1);
197 mat_struct *G_matrix_subtract(mat_struct * mt1, mat_struct * mt2)
199 return G__matrix_add(mt1, mt2, 1, -1);
214 mat_struct *G_matrix_scalar_mul(
double scalar, mat_struct *matrix, mat_struct *out)
219 if (matrix ==
NULL) {
220 G_warning (_(
"Input matrix is uninitialized"));
225 out = G_matrix_init(matrix->rows, matrix->cols, matrix->rows);
227 if (out->rows != matrix->rows || out->cols != matrix->cols)
228 out = G_matrix_resize(out, matrix->rows, matrix->cols);
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);
257 mat_struct *G_matrix_scale(mat_struct * mt1,
const double c)
259 return G__matrix_add(mt1,
NULL, c, 0);
278 mat_struct *G__matrix_add(mat_struct * mt1, mat_struct * mt2,
const double c1,
285 G_warning(_(
"First scalar multiplier must be non-zero"));
291 G_warning(_(
"One or both input matrices uninitialised"));
298 if (!((mt1->is_init) && (mt2->is_init))) {
299 G_warning(_(
"One or both input matrices uninitialised"));
303 if (mt1->rows != mt2->rows || mt1->cols != mt2->cols) {
304 G_warning(_(
"Matrix order does not match"));
309 if ((mt3 = G_matrix_init(mt1->rows, mt1->cols, mt1->ldim)) ==
NULL) {
310 G_warning(_(
"Unable to allocate space for matrix sum"));
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];
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 +
341 #if defined(HAVE_LIBBLAS)
356 mat_struct *G_matrix_product(mat_struct * mt1, mat_struct * mt2)
359 doublereal unity = 1, zero = 0;
360 integer rows, cols, interdim, lda, ldb;
361 integer1 no_trans =
'n';
363 if (!((mt1->is_init) || (mt2->is_init))) {
364 G_warning(_(
"One or both input matrices uninitialised"));
368 if (mt1->cols != mt2->rows) {
369 G_warning(_(
"Matrix order does not match"));
373 if ((mt3 = G_matrix_init(mt1->rows, mt2->cols, mt1->ldim)) ==
NULL) {
374 G_warning(_(
"Unable to allocate space for matrix product"));
380 rows = (integer) mt1->rows;
381 interdim = (integer) mt1->cols;
382 cols = (integer) mt2->cols;
384 lda = (integer) mt1->ldim;
385 ldb = (integer) mt2->ldim;
387 f77_dgemm(&no_trans, &no_trans, &rows, &cols, &interdim, &unity,
388 mt1->vals, &lda, mt2->vals, &ldb, &zero, mt3->vals, &lda);
394 #warning G_matrix_product() not compiled; requires BLAS library
410 mat_struct *G_matrix_transpose(mat_struct * mt)
414 doublereal *dbo, *dbt, *dbx, *dby;
418 if (mt->cols % 2 == 0)
423 mt1 = G_matrix_init(mt->cols, mt->rows, ldim);
430 for (cnt = 0; cnt < mt->cols; cnt++) {
434 for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
442 if (cnt < mt->cols - 1) {
452 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
480 G_matrix_LU_solve(
const mat_struct * mt1, mat_struct ** xmat0,
481 const mat_struct * bmat, mat_type mtype)
483 mat_struct *wmat, *xmat, *mtx;
485 if (mt1->is_init == 0 || bmat->is_init == 0) {
486 G_warning(_(
"Input: one or both data matrices uninitialised"));
490 if (mt1->rows != mt1->cols || mt1->rows < 1) {
491 G_warning(_(
"Principal matrix is not properly dimensioned"));
495 if (bmat->cols < 1) {
496 G_warning(_(
"Input: you must have at least one array to solve"));
501 if ((xmat = G_matrix_copy(bmat)) ==
NULL) {
502 G_warning(_(
"Could not allocate space for solution matrix"));
507 if ((mtx = G_matrix_copy(mt1)) ==
NULL) {
508 G_warning(_(
"Could not allocate space for working matrix"));
515 if ((wmat = G_matrix_copy(bmat)) ==
NULL) {
516 G_warning(_(
"Could not allocate space for working matrix"));
525 integer *perm, res_info;
526 integer num_eqns, nrhs, lda, ldb;
528 perm = (integer *) G_malloc(wmat->rows *
sizeof(integer));
531 num_eqns = (integer) mt1->rows;
532 nrhs = (integer) wmat->cols;
533 lda = (integer) mt1->ldim;
534 ldb = (integer) wmat->ldim;
537 f77_dgesv(&num_eqns, &nrhs, mtx->vals, &lda, perm, wmat->vals,
560 memcpy(xmat->vals, wmat->vals,
561 wmat->cols * wmat->ldim *
sizeof(doublereal));
569 G_warning(_(
"Matrix (or submatrix is singular). Solution undetermined"));
572 else if (res_info < 0) {
580 G_warning(_(
"Procedure not yet available for selected matrix type"));
591 #warning G_matrix_LU_solve() not compiled; requires BLAS and LAPACK libraries
594 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
608 mat_struct *G_matrix_inverse(mat_struct * mt)
610 mat_struct *mt0, *res;
613 if (mt->rows != mt->cols) {
614 G_warning(_(
"Matrix is not square. Cannot determine inverse"));
618 if ((mt0 = G_matrix_init(mt->rows, mt->rows, mt->ldim)) ==
NULL) {
619 G_warning(_(
"Unable to allocate space for matrix"));
624 for (i = 0; i < mt0->rows - 1; i++) {
625 mt0->vals[i + i * mt0->ldim] = 1.0;
627 for (j = i + 1; j < mt0->cols; j++) {
628 mt0->vals[i + j * mt0->ldim] = mt0->vals[j + i * mt0->ldim] = 0.0;
632 mt0->vals[mt0->rows - 1 + (mt0->rows - 1) * mt0->ldim] = 1.0;
635 if ((k = G_matrix_LU_solve(mt, &res, mt0, NONSYM)) == 1) {
641 G_warning(_(
"Problem in LA procedure."));
652 #warning G_matrix_inverse() not compiled; requires BLAS and LAPACK libraries
667 void G_matrix_free(mat_struct * mt)
687 void G_matrix_print(mat_struct * mt)
690 char buf[64], numbuf[64];
692 for (i = 0; i < mt->rows; i++) {
695 for (j = 0; j < mt->cols; j++) {
697 sprintf(numbuf,
"%14.6f", G_matrix_get_element(mt, i, j));
699 if (j < mt->cols - 1)
706 fprintf(stderr,
"\n");
726 int G_matrix_set_element(mat_struct * mt,
int rowval,
int colval,
double val)
729 G_warning(_(
"Element array has not been allocated"));
733 if (rowval >= mt->rows || colval >= mt->cols || rowval < 0 || colval < 0) {
734 G_warning(_(
"Specified element is outside array bounds"));
738 mt->vals[rowval + colval * mt->ldim] = (doublereal) val;
759 double G_matrix_get_element(mat_struct * mt,
int rowval,
int colval)
766 return (val = (
double)mt->vals[rowval + colval * mt->ldim]);
782 vec_struct *G_matvect_get_column(mat_struct * mt,
int col)
787 if (col < 0 || col >= mt->cols) {
788 G_warning(_(
"Specified matrix column index is outside range"));
793 G_warning(_(
"Matrix is not initialised"));
797 if ((vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) ==
NULL) {
798 G_warning(_(
"Could not allocate space for vector structure"));
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));
823 vec_struct *G_matvect_get_row(mat_struct * mt,
int row)
828 if (row < 0 || row >= mt->cols) {
829 G_warning(_(
"Specified matrix row index is outside range"));
834 G_warning(_(
"Matrix is not initialised"));
838 if ((vc1 = G_vector_init(mt->cols, mt->ldim, RVEC)) ==
NULL) {
839 G_warning(_(
"Could not allocate space for vector structure"));
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));
866 int G_matvect_extract_vector(mat_struct * mt, vtype vt,
int indx)
868 if (vt == RVEC && indx >= mt->rows) {
869 G_warning(_(
"Specified row index is outside range"));
873 else if (vt == CVEC && indx >= mt->cols) {
874 G_warning(_(
"Specified column index is outside range"));
916 int G_matvect_retrieve_matrix(vec_struct * vc)
941 vec_struct *G_matvect_product(mat_struct *A, vec_struct *
b, vec_struct *out)
943 unsigned int i, m, n, j;
944 register doublereal sum;
948 if (A->cols != b->cols) {
949 G_warning (_(
"Input matrix and vector have differing dimensions1"));
955 G_warning (_(
"Output vector is uninitialized"));
967 for (i = 0; i < m; i++) {
970 for (j = 0; j < n; j++) {
972 sum +=G_matrix_get_element(A, i, j) * G_matrix_get_element(b, 0, j);
995 vec_struct *G_vector_init(
int cells,
int ldim, vtype vt)
997 vec_struct *tmp_arry;
999 if ((cells < 1) || (vt == RVEC && ldim < 1)
1000 || (vt == CVEC && ldim < cells) || ldim < 0) {
1001 G_warning(
"Vector dimensions out of range.");
1005 tmp_arry = (vec_struct *) G_malloc(
sizeof(vec_struct));
1009 tmp_arry->cols = cells;
1010 tmp_arry->ldim = ldim;
1011 tmp_arry->type = ROWVEC_;
1014 else if (vt == CVEC) {
1015 tmp_arry->rows = cells;
1017 tmp_arry->ldim = ldim;
1018 tmp_arry->type = COLVEC_;
1021 tmp_arry->v_indx = 0;
1023 tmp_arry->vals = (doublereal *) G_calloc(ldim * tmp_arry->cols,
1024 sizeof(doublereal));
1025 tmp_arry->is_init = 1;
1042 void G_vector_free(vec_struct * v)
1065 vec_struct *G_vector_sub(vec_struct * v1, vec_struct * v2, vec_struct * out)
1067 int idx1, idx2, idx0;
1070 if (!out->is_init) {
1071 G_warning(_(
"Output vector is uninitialized"));
1075 if (v1->type != v2->type) {
1076 G_warning(_(
"Vectors are not of the same type"));
1080 if (v1->type != out->type) {
1081 G_warning(_(
"Output vector is of incorrect type"));
1085 if (v1->type == MATRIX_) {
1090 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1091 (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1092 G_warning(_(
"Vectors have differing dimensions"));
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"));
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;
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));
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));
1139 int G_vector_set(vec_struct * A,
int cells,
int ldim, vtype vt,
int vindx)
1141 if ((cells < 1) || (vt == RVEC && ldim < 1)
1142 || (vt == CVEC && ldim < cells) || ldim < 0) {
1143 G_warning(_(
"Vector dimensions out of range"));
1147 if ((vt == RVEC && vindx >= A->cols) || (vt == CVEC && vindx >= A->rows)) {
1148 G_warning(_(
"Row/column out of range"));
1170 A->vals = (doublereal *) G_calloc(ldim * A->cols,
sizeof(doublereal));
1178 #if defined(HAVE_LIBBLAS)
1192 double G_vector_norm_euclid(vec_struct * vc)
1195 doublereal *startpt;
1200 if (vc->type == ROWVEC_) {
1201 Nval = (integer) vc->cols;
1202 incr = (integer) vc->ldim;
1206 startpt = vc->vals + vc->v_indx;
1209 Nval = (integer) vc->rows;
1214 startpt = vc->vals + vc->v_indx * vc->ldim;
1218 return (
double)f77_dnrm2(&Nval, startpt, &incr);
1222 #warning G_vector_norm_euclid() not compiled; requires BLAS library
1243 double G_vector_norm_maxval(vec_struct * vc,
int vflag)
1245 doublereal xval, *startpt, *curpt;
1252 if (vc->type == ROWVEC_) {
1253 ncells = (integer) vc->cols;
1254 incr = (integer) vc->ldim;
1258 startpt = vc->vals + vc->v_indx;
1261 ncells = (integer) vc->rows;
1266 startpt = vc->vals + vc->v_indx * vc->ldim;
1272 while (ncells > 0) {
1273 if (curpt != startpt) {
1292 cellval = (double)(*curpt);
1293 if (hypot(cellval, cellval) > (double)xval)
1303 return (
double)xval;
1318 double G_vector_norm1(vec_struct * vc)
1320 double result = 0.0;
1325 G_warning(_(
"Matrix is not initialised"));
1329 idx = (vc->v_indx > 0) ? vc->v_indx : 0;
1331 if (vc->type == ROWVEC_) {
1332 for (i = 0; i < vc->cols; i++)
1333 result += fabs(G_matrix_get_element(vc, idx, i));
1336 for (i = 0; i < vc->rows; i++)
1337 result += fabs(G_matrix_get_element(vc, i, idx));
1355 vec_struct *G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out)
1357 int idx1, idx2, idx0;
1360 if (!out->is_init) {
1361 G_warning (_(
"Output vector is uninitialized"));
1365 if (v1->type != v2->type) {
1366 G_warning (_(
"Vectors are not of the same type"));
1370 if (v1->type != out->type) {
1371 G_warning (_(
"Output vector is not the same type as others"));
1375 if (v1->type == MATRIX_) {
1380 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1381 (v1->type == COLVEC_ && v1->rows != v2->rows))
1383 G_warning (_(
"Vectors have differing dimensions"));
1387 if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1388 (v1->type == COLVEC_ && v1->rows != out->rows))
1390 G_warning (_(
"Output vector has incorrect dimension"));
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);
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;
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));
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));
1429 vec_struct *G_vector_copy(
const vec_struct * vc1,
int comp_flag)
1431 vec_struct *tmp_arry;
1433 doublereal *startpt1, *startpt2, *curpt1, *curpt2;
1436 if (!vc1->is_init) {
1437 G_warning(_(
"Vector structure is not initialised"));
1441 tmp_arry = (vec_struct *) G_malloc(
sizeof(vec_struct));
1443 if (comp_flag == DO_COMPACT) {
1444 if (vc1->type == ROWVEC_) {
1446 tmp_arry->cols = vc1->cols;
1448 tmp_arry->type = ROWVEC_;
1449 tmp_arry->v_indx = 0;
1451 else if (vc1->type == COLVEC_) {
1452 tmp_arry->rows = vc1->rows;
1454 tmp_arry->ldim = vc1->ldim;
1455 tmp_arry->type = COLVEC_;
1456 tmp_arry->v_indx = 0;
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;
1471 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
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;
1487 else if (tmp_arry->type == COLVEC_) {
1488 startpt1 = tmp_arry->vals;
1489 startpt2 = vc1->vals + vc1->v_indx * vc1->ldim;
1497 G_warning(
"Structure type is not vector.");
1501 else if (comp_flag == NO_COMPACT) {
1502 startpt1 = tmp_arry->vals;
1503 startpt2 = vc1->vals;
1508 cnt = vc1->ldim * vc1->cols;
1511 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1516 memcpy(curpt1, curpt2,
sizeof(doublereal));
1522 tmp_arry->is_init = 1;
1542 int G_matrix_read(FILE * fp, mat_struct * out)
1551 if (!
G_getl(buff,
sizeof(buff), fp))
1557 if (sscanf(buff,
"Matrix: %d by %d", &rows, &cols) != 2) {
1562 G_matrix_set(out, rows, cols, rows);
1564 for (i = 0; i < rows; i++) {
1565 if (fscanf(fp,
"row%d:", &row) != 1 || row != i) {
1569 for (j = 0; j < cols; j++) {
1570 if (fscanf(fp,
"%lf:", &val) != 1) {
1575 G_matrix_set_element(out, i, j, val);
1596 mat_struct *G_matrix_resize(mat_struct *in,
int rows,
int cols)
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));
1606 int old_size = in->rows * in->cols;
1607 int new_size = rows * cols;
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);
1630 int G_matrix_stdin(mat_struct * out)
1632 return G_matrix_read(stdin, out);
1648 int G_matrix_eigen_sort(vec_struct * d, mat_struct * m)
1654 G_matrix_set(&tmp, m->rows + 1, m->cols, m->ldim + 1);
1656 idx = (d->v_indx > 0) ? d->v_indx : 0;
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));
1666 G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, i, idx));
1670 qsort(tmp.vals, tmp.cols, tmp.ldim *
sizeof(doublereal), egcmp);
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));
1680 G_matrix_set_element(d, i, idx, G_matrix_get_element(&tmp, 0, i));
1689 static int egcmp(
const void *pa,
const void *pb)
1691 double a = *(doublereal *
const)pa;
1692 double b = *(doublereal *
const)pb;
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
if(!DBFLoadRecord(psDBF, hEntity)) return NULL
int G_getl(char *buf, int n, FILE *fd)
Gets a line of text from a file.
void G_message(const char *msg,...)
Print a message to stderr.
void G_free(void *buf)
Free allocated memory.
void G_warning(const char *msg,...)
Print a warning message to stderr.