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);
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;
93 memset(A->vals, 0,
sizeof(A->vals));
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));
150 G_warning(_(
"Matrix is not initialised fully."));
155 G_warning(_(
"Unable to allocate space for matrix copy"));
159 memcpy(&B->vals[0], &A->vals[0], A->cols * A->ldim *
sizeof(doublereal));
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)
233 for (i = 0; i < m; i++) {
234 for (j = 0; j < n; j++) {
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"));
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)
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"));
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
414 doublereal *dbo, *dbt, *dbx, *dby;
418 if (mt->cols % 2 == 0)
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)
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"));
502 G_warning(_(
"Could not allocate space for solution matrix"));
508 G_warning(_(
"Could not allocate space for working matrix"));
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)
610 mat_struct *mt0, *res;
613 if (mt->rows != mt->cols) {
614 G_warning(_(
"Matrix is not square. Cannot determine inverse"));
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;
641 G_warning(_(
"Problem in LA procedure."));
652 #warning G_matrix_inverse() not compiled; requires BLAS and LAPACK libraries
690 char buf[64], numbuf[64];
692 for (i = 0; i < mt->rows; i++) {
695 for (j = 0; j < mt->cols; j++) {
699 if (j < mt->cols - 1)
706 fprintf(stderr,
"\n");
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;
766 return (val = (
double)mt->vals[rowval + colval * mt->ldim]);
787 if (col < 0 || col >= mt->cols) {
788 G_warning(_(
"Specified matrix column index is outside range"));
793 G_warning(_(
"Matrix is not initialised"));
798 G_warning(_(
"Could not allocate space for vector structure"));
802 for (i = 0; i < mt->rows; i++)
828 if (row < 0 || row >= mt->cols) {
829 G_warning(_(
"Specified matrix row index is outside range"));
834 G_warning(_(
"Matrix is not initialised"));
839 G_warning(_(
"Could not allocate space for vector structure"));
843 for (i = 0; i < mt->cols; i++)
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"));
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++) {
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;
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++)
1113 for (i = 0; i < v1->rows; i++)
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)
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
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;
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++)
1336 for (i = 0; i < vc->rows; i++)
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++)
1407 for (i = 0; i < v1->rows; i++)
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;
1551 if (!
G_getl(buff,
sizeof(buff), fp))
1557 if (sscanf(buff,
"Matrix: %d by %d", &rows, &cols) != 2) {
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) {
1600 int i, j, p, index = 0;
1601 for (i = 0; i < rows; i++)
1602 for (j = 0; j < cols; 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++)
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++)
1663 if (d->type == ROWVEC_)
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++)
1677 if (d->type == ROWVEC_)
1689 static int egcmp(
const void *pa,
const void *pb)
1691 double a = *(doublereal *
const)pa;
1692 double b = *(doublereal *
const)pb;