34 #include <grass/config.h>
36 #if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS) && defined(HAVE_G2C_H)
38 #include <grass/gis.h>
39 #include <grass/glocale.h>
43 static int egcmp(
const void *pa,
const void *pb);
59 mat_struct *G_matrix_init(
int rows,
int cols,
int ldim)
63 if (rows < 1 || cols < 1 || ldim < rows) {
64 G_warning(_(
"Matrix dimensions out of range"));
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;
75 tmp_arry->vals = (doublereal *) G_calloc(ldim * cols,
sizeof(doublereal));
76 tmp_arry->is_init = 1;
91 int G_matrix_zero(mat_struct * A)
96 memset(A->vals, 0,
sizeof(A->vals));
117 int G_matrix_set(mat_struct * A,
int rows,
int cols,
int ldim)
119 if (rows < 1 || cols < 1 || ldim < 0) {
120 G_warning(_(
"Matrix dimensions out of range"));
130 A->vals = (doublereal *) G_calloc(ldim * cols,
sizeof(doublereal));
148 mat_struct *G_matrix_copy(
const mat_struct * A)
153 G_warning(_(
"Matrix is not initialised fully."));
157 if ((B = G_matrix_init(A->rows, A->cols, A->ldim)) ==
NULL) {
158 G_warning(_(
"Unable to allocate space for matrix copy"));
162 memcpy(&B->vals[0], &A->vals[0], A->cols * A->ldim *
sizeof(doublereal));
181 mat_struct *G_matrix_add(mat_struct * mt1, mat_struct * mt2)
183 return G__matrix_add(mt1, mt2, 1, 1);
200 mat_struct *G_matrix_subtract(mat_struct * mt1, mat_struct * mt2)
202 return G__matrix_add(mt1, mt2, 1, -1);
219 mat_struct *G_matrix_scale(mat_struct * mt1,
const double c)
221 return G__matrix_add(mt1,
NULL, c, 0);
240 mat_struct *G__matrix_add(mat_struct * mt1, mat_struct * mt2,
const double c1,
247 G_warning(_(
"First scalar multiplier must be non-zero"));
253 G_warning(_(
"One or both input matrices uninitialised"));
260 if (!((mt1->is_init) && (mt2->is_init))) {
261 G_warning(_(
"One or both input matrices uninitialised"));
265 if (mt1->rows != mt2->rows || mt1->cols != mt2->cols) {
266 G_warning(_(
"Matrix order does not match"));
271 if ((mt3 = G_matrix_init(mt1->rows, mt1->cols, mt1->ldim)) ==
NULL) {
272 G_warning(_(
"Unable to allocate space for matrix sum"));
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];
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 +
303 #if defined(HAVE_LIBBLAS)
318 mat_struct *G_matrix_product(mat_struct * mt1, mat_struct * mt2)
321 doublereal unity = 1, zero = 0;
322 integer rows,
cols, interdim, lda, ldb;
323 integer1 no_trans =
'n';
325 if (!((mt1->is_init) || (mt2->is_init))) {
326 G_warning(_(
"One or both input matrices uninitialised"));
330 if (mt1->cols != mt2->rows) {
331 G_warning(_(
"Matrix order does not match"));
335 if ((mt3 = G_matrix_init(mt1->rows, mt2->cols, mt1->ldim)) ==
NULL) {
336 G_warning(_(
"Unable to allocate space for matrix product"));
342 rows = (integer) mt1->rows;
343 interdim = (integer) mt1->cols;
344 cols = (integer) mt2->cols;
346 lda = (integer) mt1->ldim;
347 ldb = (integer) mt2->ldim;
349 f77_dgemm(&no_trans, &no_trans, &rows, &cols, &interdim, &unity,
350 mt1->vals, &lda, mt2->vals, &ldb, &zero, mt3->vals, &lda);
356 #warning G_matrix_product() not compiled; requires BLAS library
372 mat_struct *G_matrix_transpose(mat_struct * mt)
376 doublereal *dbo, *dbt, *dbx, *dby;
380 if (mt->cols % 2 == 0)
385 mt1 = G_matrix_init(mt->cols, mt->rows, ldim);
392 for (cnt = 0; cnt < mt->cols; cnt++) {
396 for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
404 if (cnt < mt->cols - 1) {
414 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
442 G_matrix_LU_solve(
const mat_struct * mt1, mat_struct ** xmat0,
443 const mat_struct * bmat, mat_type mtype)
445 mat_struct *wmat, *xmat, *mtx;
447 if (mt1->is_init == 0 || bmat->is_init == 0) {
448 G_warning(_(
"Input: one or both data matrices uninitialised"));
452 if (mt1->rows != mt1->cols || mt1->rows < 1) {
453 G_warning(_(
"Principal matrix is not properly dimensioned"));
457 if (bmat->cols < 1) {
458 G_warning(_(
"Input: you must have at least one array to solve"));
463 if ((xmat = G_matrix_copy(bmat)) ==
NULL) {
464 G_warning(_(
"Could not allocate space for solution matrix"));
469 if ((mtx = G_matrix_copy(mt1)) ==
NULL) {
470 G_warning(_(
"Could not allocate space for working matrix"));
477 if ((wmat = G_matrix_copy(bmat)) ==
NULL) {
478 G_warning(_(
"Could not allocate space for working matrix"));
487 integer *perm, res_info;
488 integer num_eqns, nrhs, lda, ldb;
490 perm = (integer *) G_malloc(wmat->rows);
493 num_eqns = (integer) mt1->rows;
494 nrhs = (integer) wmat->cols;
495 lda = (integer) mt1->ldim;
496 ldb = (integer) wmat->ldim;
499 f77_dgesv(&num_eqns, &nrhs, mtx->vals, &lda, perm, wmat->vals,
522 memcpy(xmat->vals, wmat->vals,
523 wmat->cols * wmat->ldim *
sizeof(doublereal));
531 G_warning(_(
"Matrix (or submatrix is singular). Solution undetermined"));
534 else if (res_info < 0) {
542 G_warning(_(
"Procedure not yet available for selected matrix type"));
553 #warning G_matrix_LU_solve() not compiled; requires BLAS and LAPACK libraries
556 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
570 mat_struct *G_matrix_inverse(mat_struct * mt)
572 mat_struct *mt0, *res;
575 if (mt->rows != mt->cols) {
576 G_warning(_(
"Matrix is not square. Cannot determine inverse"));
580 if ((mt0 = G_matrix_init(mt->rows, mt->rows, mt->ldim)) ==
NULL) {
581 G_warning(_(
"Unable to allocate space for matrix"));
586 for (i = 0; i < mt0->rows - 1; i++) {
587 mt0->vals[i + i * mt0->ldim] = 1.0;
589 for (j = i + 1; j < mt0->cols; j++) {
590 mt0->vals[i + j * mt0->ldim] = mt0->vals[j + i * mt0->ldim] = 0.0;
594 mt0->vals[mt0->rows - 1 + (mt0->rows - 1) * mt0->ldim] = 1.0;
597 if ((k = G_matrix_LU_solve(mt, &res, mt0, NONSYM)) == 1) {
603 G_warning(_(
"Problem in LA procedure."));
614 #warning G_matrix_inverse() not compiled; requires BLAS and LAPACK libraries
629 void G_matrix_free(mat_struct * mt)
649 void G_matrix_print(mat_struct * mt)
652 char buf[64], numbuf[64];
654 for (i = 0; i < mt->rows; i++) {
657 for (j = 0; j < mt->cols; j++) {
659 sprintf(numbuf,
"%14.6f", G_matrix_get_element(mt, i, j));
661 if (j < mt->cols - 1)
668 fprintf(stderr,
"\n");
688 int G_matrix_set_element(mat_struct * mt,
int rowval,
int colval,
double val)
691 G_warning(_(
"Element array has not been allocated"));
695 if (rowval >= mt->rows || colval >= mt->cols || rowval < 0 || colval < 0) {
696 G_warning(_(
"Specified element is outside array bounds"));
700 mt->vals[rowval + colval * mt->ldim] = (doublereal) val;
721 double G_matrix_get_element(mat_struct * mt,
int rowval,
int colval)
728 return (val = (
double)mt->vals[rowval + colval * mt->ldim]);
744 vec_struct *G_matvect_get_column(mat_struct * mt,
int col)
749 if (col < 0 || col >= mt->cols) {
750 G_warning(_(
"Specified matrix column index is outside range"));
755 G_warning(_(
"Matrix is not initialised"));
759 if ((vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) ==
NULL) {
760 G_warning(_(
"Could not allocate space for vector structure"));
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));
785 vec_struct *G_matvect_get_row(mat_struct * mt,
int row)
790 if (row < 0 || row >= mt->cols) {
791 G_warning(_(
"Specified matrix row index is outside range"));
796 G_warning(_(
"Matrix is not initialised"));
800 if ((vc1 = G_vector_init(mt->cols, mt->ldim, RVEC)) ==
NULL) {
801 G_warning(_(
"Could not allocate space for vector structure"));
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));
828 int G_matvect_extract_vector(mat_struct * mt, vtype vt,
int indx)
830 if (vt == RVEC && indx >= mt->rows) {
831 G_warning(_(
"Specified row index is outside range"));
835 else if (vt == CVEC && indx >= mt->cols) {
836 G_warning(_(
"Specified column index is outside range"));
878 int G_matvect_retrieve_matrix(vec_struct * vc)
905 vec_struct *G_vector_init(
int cells,
int ldim, vtype vt)
907 vec_struct *tmp_arry;
909 if ((cells < 1) || (vt == RVEC && ldim < 1)
910 || (vt == CVEC && ldim < cells) || ldim < 0) {
911 G_warning(
"Vector dimensions out of range.");
915 tmp_arry = (vec_struct *) G_malloc(
sizeof(vec_struct));
919 tmp_arry->cols = cells;
920 tmp_arry->ldim = ldim;
921 tmp_arry->type = ROWVEC_;
924 else if (vt == CVEC) {
925 tmp_arry->rows = cells;
927 tmp_arry->ldim = ldim;
928 tmp_arry->type = COLVEC_;
931 tmp_arry->v_indx = 0;
933 tmp_arry->vals = (doublereal *) G_calloc(ldim * tmp_arry->cols,
935 tmp_arry->is_init = 1;
952 void G_vector_free(vec_struct * v)
975 vec_struct *G_vector_sub(vec_struct * v1, vec_struct * v2, vec_struct * out)
977 int idx1, idx2, idx0;
981 G_warning(_(
"Output vector is uninitialized"));
985 if (v1->type != v2->type) {
986 G_warning(_(
"Vectors are not of the same type"));
990 if (v1->type != out->type) {
991 G_warning(_(
"Output vector is of incorrect type"));
995 if (v1->type == MATRIX_) {
1000 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1001 (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1002 G_warning(_(
"Vectors have differing dimensions"));
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"));
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;
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));
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));
1050 int G_vector_set(vec_struct * A,
int cells,
int ldim, vtype vt,
int vindx)
1052 if ((cells < 1) || (vt == RVEC && ldim < 1)
1053 || (vt == CVEC && ldim < cells) || ldim < 0) {
1054 G_warning(_(
"Vector dimensions out of range"));
1058 if ((vt == RVEC && vindx >= A->cols) || (vt == CVEC && vindx >= A->rows)) {
1059 G_warning(_(
"Row/column out of range"));
1081 A->vals = (doublereal *) G_calloc(ldim * A->cols,
sizeof(doublereal));
1089 #if defined(HAVE_LIBBLAS)
1103 double G_vector_norm_euclid(vec_struct * vc)
1106 doublereal *startpt;
1111 if (vc->type == ROWVEC_) {
1112 Nval = (integer) vc->cols;
1113 incr = (integer) vc->ldim;
1117 startpt = vc->vals + vc->v_indx;
1120 Nval = (integer) vc->rows;
1125 startpt = vc->vals + vc->v_indx * vc->ldim;
1129 return (
double)f77_dnrm2(&Nval, startpt, &incr);
1133 #warning G_vector_norm_euclid() not compiled; requires BLAS library
1154 double G_vector_norm_maxval(vec_struct * vc,
int vflag)
1156 doublereal xval, *startpt, *curpt;
1163 if (vc->type == ROWVEC_) {
1164 ncells = (integer) vc->cols;
1165 incr = (integer) vc->ldim;
1169 startpt = vc->vals + vc->v_indx;
1172 ncells = (integer) vc->rows;
1177 startpt = vc->vals + vc->v_indx * vc->ldim;
1183 while (ncells > 0) {
1184 if (curpt != startpt) {
1203 cellval = (double)(*curpt);
1204 if (hypot(cellval, cellval) > (double)xval)
1214 return (
double)xval;
1229 double G_vector_norm1(vec_struct * vc)
1231 double result = 0.0;
1236 G_warning(_(
"Matrix is not initialised"));
1240 idx = (vc->v_indx > 0) ? vc->v_indx : 0;
1242 if (vc->type == ROWVEC_) {
1243 for (i = 0; i < vc->cols; i++)
1244 result += fabs(G_matrix_get_element(vc, idx, i));
1247 for (i = 0; i < vc->rows; i++)
1248 result += fabs(G_matrix_get_element(vc, i, idx));
1266 vec_struct *G_vector_copy(
const vec_struct * vc1,
int comp_flag)
1268 vec_struct *tmp_arry;
1270 doublereal *startpt1, *startpt2, *curpt1, *curpt2;
1273 if (!vc1->is_init) {
1274 G_warning(_(
"Vector structure is not initialised"));
1278 tmp_arry = (vec_struct *) G_malloc(
sizeof(vec_struct));
1280 if (comp_flag == DO_COMPACT) {
1281 if (vc1->type == ROWVEC_) {
1283 tmp_arry->cols = vc1->cols;
1285 tmp_arry->type = ROWVEC_;
1286 tmp_arry->v_indx = 0;
1288 else if (vc1->type == COLVEC_) {
1289 tmp_arry->rows = vc1->rows;
1291 tmp_arry->ldim = vc1->ldim;
1292 tmp_arry->type = COLVEC_;
1293 tmp_arry->v_indx = 0;
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;
1308 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
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;
1324 else if (tmp_arry->type == COLVEC_) {
1325 startpt1 = tmp_arry->vals;
1326 startpt2 = vc1->vals + vc1->v_indx * vc1->ldim;
1334 G_warning(
"Structure type is not vector.");
1338 else if (comp_flag == NO_COMPACT) {
1339 startpt1 = tmp_arry->vals;
1340 startpt2 = vc1->vals;
1345 cnt = vc1->ldim * vc1->cols;
1348 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1353 memcpy(curpt1, curpt2,
sizeof(doublereal));
1359 tmp_arry->is_init = 1;
1379 int G_matrix_read(FILE * fp, mat_struct * out)
1388 if (!
G_getl(buff,
sizeof(buff), fp))
1394 if (sscanf(buff,
"Matrix: %d by %d", &rows, &cols) != 2) {
1399 G_matrix_set(out, rows, cols, rows);
1401 for (i = 0; i < rows; i++) {
1402 if (fscanf(fp,
"row%d:", &row) != 1 || row != i) {
1406 for (j = 0; j <
cols; j++) {
1407 if (fscanf(fp,
"%lf:", &val) != 1) {
1412 G_matrix_set_element(out, i, j, val);
1433 int G_matrix_stdin(mat_struct * out)
1435 return G_matrix_read(stdin, out);
1451 int G_matrix_eigen_sort(vec_struct * d, mat_struct * m)
1457 G_matrix_set(&tmp, m->rows + 1, m->cols, m->ldim + 1);
1459 idx = (d->v_indx > 0) ? d->v_indx : 0;
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));
1469 G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, i, idx));
1473 qsort(tmp.vals, tmp.cols, tmp.ldim *
sizeof(doublereal), egcmp);
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));
1483 G_matrix_set_element(d, i, idx, G_matrix_get_element(&tmp, 0, i));
1492 static int egcmp(
const void *pa,
const void *pb)
1494 double a = *(doublereal *
const)pa;
1495 double b = *(doublereal *
const)pb;