30#include <grass/config.h>
32#if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS)
35#include <grass/glocale.h>
38static int egcmp(
const void *pa,
const void *pb);
57 if (rows < 1 || cols < 1 || ldim < rows) {
58 G_warning(_(
"Matrix dimensions out of range"));
62 tmp_arry = (mat_struct *)G_malloc(
sizeof(mat_struct));
63 tmp_arry->rows = rows;
64 tmp_arry->cols = cols;
65 tmp_arry->ldim = ldim;
66 tmp_arry->type = MATRIX_;
67 tmp_arry->v_indx = -1;
69 tmp_arry->vals = (doublereal *)G_calloc(ldim * cols,
sizeof(doublereal));
70 tmp_arry->is_init = 1;
89 memset(A->vals, 0, (A->ldim * A->cols) *
sizeof(doublereal));
111 if (rows < 1 || cols < 1 || ldim < 0) {
112 G_warning(_(
"Matrix dimensions out of range"));
122 A->vals = (doublereal *)G_calloc(ldim * cols,
sizeof(doublereal));
144 G_warning(_(
"Matrix is not initialised fully."));
149 G_warning(_(
"Unable to allocate space for matrix copy"));
153 memcpy(&B->vals[0], &A->vals[0], A->cols * A->ldim *
sizeof(doublereal));
212 if (matrix ==
NULL) {
213 G_warning(_(
"Input matrix is uninitialized"));
218 out =
G_matrix_init(matrix->rows, matrix->cols, matrix->rows);
220 if (out->rows != matrix->rows || out->cols != matrix->cols)
226 for (i = 0; i < m; i++) {
227 for (j = 0; j < n; j++) {
278 G_warning(_(
"First scalar multiplier must be non-zero"));
284 G_warning(_(
"One or both input matrices uninitialised"));
291 if (!((mt1->is_init) && (mt2->is_init))) {
292 G_warning(_(
"One or both input matrices uninitialised"));
296 if (mt1->rows != mt2->rows || mt1->cols != mt2->cols) {
297 G_warning(_(
"Matrix order does not match"));
303 G_warning(_(
"Unable to allocate space for matrix sum"));
309 for (i = 0; i < mt3->rows; i++) {
310 for (j = 0; j < mt3->cols; j++) {
311 mt3->vals[i + mt3->ldim * j] =
312 c1 * mt1->vals[i + mt1->ldim * j];
319 for (i = 0; i < mt3->rows; i++) {
320 for (j = 0; j < mt3->cols; j++) {
321 mt3->vals[i + mt3->ldim * j] =
322 c1 * mt1->vals[i + mt1->ldim * j] +
323 c2 * mt2->vals[i + mt2->ldim * j];
331#if defined(HAVE_LIBBLAS)
349 doublereal unity = 1, zero = 0;
350 integer rows, cols, interdim, lda, ldb;
351 integer1 no_trans =
'n';
353 if (!((mt1->is_init) || (mt2->is_init))) {
354 G_warning(_(
"One or both input matrices uninitialised"));
358 if (mt1->cols != mt2->rows) {
359 G_warning(_(
"Matrix order does not match"));
364 G_warning(_(
"Unable to allocate space for matrix product"));
370 rows = (integer)mt1->rows;
371 interdim = (integer)mt1->cols;
372 cols = (integer)mt2->cols;
374 lda = (integer)mt1->ldim;
375 ldb = (integer)mt2->ldim;
377 f77_dgemm(&no_trans, &no_trans, &rows, &cols, &interdim, &unity, mt1->vals,
378 &lda, mt2->vals, &ldb, &zero, mt3->vals, &lda);
384#warning G_matrix_product() not compiled; requires BLAS library
404 doublereal *dbo, *dbt, *dbx, *dby;
408 if (mt->cols % 2 == 0)
420 for (cnt = 0; cnt < mt->cols; cnt++) {
424 for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
432 if (cnt < mt->cols - 1) {
441#if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
469 const mat_struct *bmat, mat_type mtype)
471 mat_struct *wmat, *xmat, *mtx;
473 if (mt1->is_init == 0 || bmat->is_init == 0) {
474 G_warning(_(
"Input: one or both data matrices uninitialised"));
478 if (mt1->rows != mt1->cols || mt1->rows < 1) {
479 G_warning(_(
"Principal matrix is not properly dimensioned"));
483 if (bmat->cols < 1) {
484 G_warning(_(
"Input: you must have at least one array to solve"));
490 G_warning(_(
"Could not allocate space for solution matrix"));
496 G_warning(_(
"Could not allocate space for working matrix"));
504 G_warning(_(
"Could not allocate space for working matrix"));
512 integer *perm, res_info;
513 integer num_eqns, nrhs, lda, ldb;
515 perm = (integer *)G_malloc(wmat->rows *
sizeof(integer));
518 num_eqns = (integer)mt1->rows;
519 nrhs = (integer)wmat->cols;
520 lda = (integer)mt1->ldim;
521 ldb = (integer)wmat->ldim;
524 f77_dgesv(&num_eqns, &nrhs, mtx->vals, &lda, perm, wmat->vals, &ldb,
547 memcpy(xmat->vals, wmat->vals,
548 wmat->cols * wmat->ldim *
sizeof(doublereal));
557 _(
"Matrix (or submatrix is singular). Solution undetermined"));
560 else if (res_info < 0) {
567 G_warning(_(
"Procedure not yet available for selected matrix type"));
578#warning G_matrix_LU_solve() not compiled; requires BLAS and LAPACK libraries
581#if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
597 mat_struct *mt0, *res;
600 if (mt->rows != mt->cols) {
601 G_warning(_(
"Matrix is not square. Cannot determine inverse"));
606 G_warning(_(
"Unable to allocate space for matrix"));
611 for (i = 0; i < mt0->rows - 1; i++) {
612 mt0->vals[i + i * mt0->ldim] = 1.0;
614 for (j = i + 1; j < mt0->cols; j++) {
615 mt0->vals[i + j * mt0->ldim] = mt0->vals[j + i * mt0->ldim] = 0.0;
619 mt0->vals[mt0->rows - 1 + (mt0->rows - 1) * mt0->ldim] = 1.0;
628 G_warning(_(
"Problem in LA procedure."));
639#warning G_matrix_inverse() not compiled; requires BLAS and LAPACK libraries
675 char buf[64], numbuf[64];
677 for (i = 0; i < mt->rows; i++) {
680 for (j = 0; j < mt->cols; j++) {
684 if (j < mt->cols - 1)
691 fprintf(stderr,
"\n");
714 G_warning(_(
"Element array has not been allocated"));
718 if (rowval >= mt->rows || colval >= mt->cols || rowval < 0 || colval < 0) {
719 G_warning(_(
"Specified element is outside array bounds"));
723 mt->vals[rowval + colval * mt->ldim] = (doublereal)val;
750 return (val = (
double)mt->vals[rowval + colval * mt->ldim]);
770 if (col < 0 || col >= mt->cols) {
771 G_warning(_(
"Specified matrix column index is outside range"));
776 G_warning(_(
"Matrix is not initialised"));
781 G_warning(_(
"Could not allocate space for vector structure"));
785 for (i = 0; i < mt->rows; i++)
810 if (row < 0 || row >= mt->cols) {
811 G_warning(_(
"Specified matrix row index is outside range"));
816 G_warning(_(
"Matrix is not initialised"));
821 G_warning(_(
"Could not allocate space for vector structure"));
825 for (i = 0; i < mt->cols; i++)
849 if (vt == RVEC && indx >= mt->rows) {
850 G_warning(_(
"Specified row index is outside range"));
854 else if (vt == CVEC && indx >= mt->cols) {
855 G_warning(_(
"Specified column index is outside range"));
920 unsigned int i, m, n, j;
921 register doublereal sum;
925 if (A->cols !=
b->cols) {
926 G_warning(_(
"Input matrix and vector have differing dimensions1"));
931 G_warning(_(
"Output vector is uninitialized"));
943 for (i = 0; i < m; i++) {
947 for (j = 0; j < n; j++) {
974 vec_struct *tmp_arry;
976 if ((cells < 1) || (vt == RVEC && ldim < 1) ||
977 (vt == CVEC && ldim < cells) || ldim < 0) {
978 G_warning(
"Vector dimensions out of range.");
982 tmp_arry = (vec_struct *)G_malloc(
sizeof(vec_struct));
986 tmp_arry->cols = cells;
987 tmp_arry->ldim = ldim;
988 tmp_arry->type = ROWVEC_;
991 else if (vt == CVEC) {
992 tmp_arry->rows = cells;
994 tmp_arry->ldim = ldim;
995 tmp_arry->type = COLVEC_;
998 tmp_arry->v_indx = 0;
1001 (doublereal *)G_calloc(ldim * tmp_arry->cols,
sizeof(doublereal));
1002 tmp_arry->is_init = 1;
1043 int idx1, idx2, idx0;
1046 if (!out->is_init) {
1047 G_warning(_(
"Output vector is uninitialized"));
1051 if (v1->type != v2->type) {
1052 G_warning(_(
"Vectors are not of the same type"));
1056 if (v1->type != out->type) {
1057 G_warning(_(
"Output vector is of incorrect type"));
1061 if (v1->type == MATRIX_) {
1066 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1067 (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1068 G_warning(_(
"Vectors have differing dimensions"));
1072 if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1073 (v1->type == COLVEC_ && v1->rows != out->rows)) {
1074 G_warning(_(
"Output vector has incorrect dimension"));
1078 idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
1079 idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
1080 idx0 = (out->v_indx > 0) ? out->v_indx : 0;
1082 if (v1->type == ROWVEC_) {
1083 for (i = 0; i < v1->cols; i++)
1089 for (i = 0; i < v1->rows; i++)
1118 if ((cells < 1) || (vt == RVEC && ldim < 1) ||
1119 (vt == CVEC && ldim < cells) || ldim < 0) {
1120 G_warning(_(
"Vector dimensions out of range"));
1124 if ((vt == RVEC && vindx >= A->cols) || (vt == CVEC && vindx >= A->rows)) {
1125 G_warning(_(
"Row/column out of range"));
1147 A->vals = (doublereal *)G_calloc(ldim * A->cols,
sizeof(doublereal));
1153#if defined(HAVE_LIBBLAS)
1170 doublereal *startpt;
1175 if (vc->type == ROWVEC_) {
1176 Nval = (integer)vc->cols;
1177 incr = (integer)vc->ldim;
1181 startpt = vc->vals + vc->v_indx;
1184 Nval = (integer)vc->rows;
1189 startpt = vc->vals + vc->v_indx * vc->ldim;
1193 return (
double)f77_dnrm2(&Nval, startpt, &incr);
1197#warning G_vector_norm_euclid() not compiled; requires BLAS library
1219 doublereal xval, *startpt, *curpt;
1226 if (vc->type == ROWVEC_) {
1227 ncells = (integer)vc->cols;
1228 incr = (integer)vc->ldim;
1232 startpt = vc->vals + vc->v_indx;
1235 ncells = (integer)vc->rows;
1240 startpt = vc->vals + vc->v_indx * vc->ldim;
1246 while (ncells > 0) {
1247 if (curpt != startpt) {
1263 cellval = (double)(*curpt);
1264 if (hypot(cellval, cellval) > (double)xval)
1274 return (
double)xval;
1290 double result = 0.0;
1295 G_warning(_(
"Matrix is not initialised"));
1299 idx = (vc->v_indx > 0) ? vc->v_indx : 0;
1301 if (vc->type == ROWVEC_) {
1302 for (i = 0; i < vc->cols; i++)
1306 for (i = 0; i < vc->rows; i++)
1328 int idx1, idx2, idx0;
1331 if (!out->is_init) {
1332 G_warning(_(
"Output vector is uninitialized"));
1336 if (v1->type != v2->type) {
1337 G_warning(_(
"Vectors are not of the same type"));
1341 if (v1->type != out->type) {
1342 G_warning(_(
"Output vector is not the same type as others"));
1346 if (v1->type == MATRIX_) {
1351 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1352 (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1353 G_warning(_(
"Vectors have differing dimensions"));
1357 if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1358 (v1->type == COLVEC_ && v1->rows != out->rows)) {
1359 G_warning(_(
"Output vector has incorrect dimension"));
1363#if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS)
1364 f77_dhad(v1->cols, 1.0, v1->vals, 1, v2->vals, 1, 0.0, out->vals, 1.0);
1366 idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
1367 idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
1368 idx0 = (out->v_indx > 0) ? out->v_indx : 0;
1370 if (v1->type == ROWVEC_) {
1371 for (i = 0; i < v1->cols; i++)
1377 for (i = 0; i < v1->rows; i++)
1400 vec_struct *tmp_arry;
1402 doublereal *startpt1, *startpt2, *curpt1, *curpt2;
1405 if (!vc1->is_init) {
1406 G_warning(_(
"Vector structure is not initialised"));
1410 tmp_arry = (vec_struct *)G_malloc(
sizeof(vec_struct));
1412 if (comp_flag == DO_COMPACT) {
1413 if (vc1->type == ROWVEC_) {
1415 tmp_arry->cols = vc1->cols;
1417 tmp_arry->type = ROWVEC_;
1418 tmp_arry->v_indx = 0;
1420 else if (vc1->type == COLVEC_) {
1421 tmp_arry->rows = vc1->rows;
1423 tmp_arry->ldim = vc1->ldim;
1424 tmp_arry->type = COLVEC_;
1425 tmp_arry->v_indx = 0;
1432 else if (comp_flag == NO_COMPACT) {
1433 tmp_arry->v_indx = vc1->v_indx;
1434 tmp_arry->rows = vc1->rows;
1435 tmp_arry->cols = vc1->cols;
1436 tmp_arry->ldim = vc1->ldim;
1437 tmp_arry->type = vc1->type;
1440 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1444 tmp_arry->vals = (doublereal *)G_calloc(tmp_arry->ldim * tmp_arry->cols,
1445 sizeof(doublereal));
1446 if (comp_flag == DO_COMPACT) {
1447 if (tmp_arry->type == ROWVEC_) {
1448 startpt1 = tmp_arry->vals;
1449 startpt2 = vc1->vals + vc1->v_indx;
1456 else if (tmp_arry->type == COLVEC_) {
1457 startpt1 = tmp_arry->vals;
1458 startpt2 = vc1->vals + vc1->v_indx * vc1->ldim;
1466 G_warning(
"Structure type is not vector.");
1470 else if (comp_flag == NO_COMPACT) {
1471 startpt1 = tmp_arry->vals;
1472 startpt2 = vc1->vals;
1477 cnt = vc1->ldim * vc1->cols;
1480 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1485 memcpy(curpt1, curpt2,
sizeof(doublereal));
1491 tmp_arry->is_init = 1;
1519 if (!
G_getl(buff,
sizeof(buff), fp))
1525 if (sscanf(buff,
"Matrix: %d by %d", &rows, &cols) != 2) {
1532 for (i = 0; i < rows; i++) {
1533 if (fscanf(fp,
"row%d:", &row) != 1 || row != i) {
1537 for (j = 0; j < cols; j++) {
1538 if (fscanf(fp,
"%lf:", &val) != 1) {
1570 for (i = 0; i < rows; i++)
1571 for (j = 0; j < cols; j++)
1575 int old_size = in->rows * in->cols;
1576 int new_size = rows * cols;
1578 if (new_size > old_size)
1579 for (p = old_size; p < new_size; p++)
1623 idx = (d->v_indx > 0) ? d->v_indx : 0;
1626 for (i = 0; i < m->cols; i++) {
1627 for (j = 0; j < m->rows; j++)
1629 if (d->type == ROWVEC_)
1636 qsort(tmp.vals, tmp.cols, tmp.ldim *
sizeof(doublereal), egcmp);
1639 for (i = 0; i < m->cols; i++) {
1640 for (j = 0; j < m->rows; j++)
1642 if (d->type == ROWVEC_)
1653static int egcmp(
const void *pa,
const void *pb)
1655 double a = *(doublereal *
const)pa;
1656 double b = *(doublereal *
const)pb;
void G_free(void *buf)
Free allocated memory.
int G_getl(char *buf, int n, FILE *fd)
Gets a line of text from a file.
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
void G_message(const char *msg,...)
Print a message to stderr.
void G_warning(const char *msg,...)
Print a warning message to stderr.
vec_struct * G_matvect_get_column(mat_struct *mt, int col)
Retrieve a column of the matrix to a vector structure.
mat_struct * G_matrix_copy(const mat_struct *A)
Copy a matrix.
vec_struct * G_vector_sub(vec_struct *v1, vec_struct *v2, vec_struct *out)
Subtract two vectors.
double G_vector_norm_maxval(vec_struct *vc, int vflag)
Calculates maximum value.
mat_struct * G__matrix_add(mat_struct *mt1, mat_struct *mt2, const double c1, const double c2)
vec_struct * G_vector_product(vec_struct *v1, vec_struct *v2, vec_struct *out)
Calculates the vector product.
mat_struct * G_matrix_subtract(mat_struct *mt1, mat_struct *mt2)
Subtract two matricies.
int G_matvect_retrieve_matrix(vec_struct *vc)
Revert a vector to matrix.
double G_vector_norm_euclid(vec_struct *vc)
Calculates euclidean norm.
vec_struct * G_vector_init(int cells, int ldim, vtype vt)
Initialize a vector structure.
mat_struct * G_matrix_inverse(mat_struct *mt)
Returns the matrix inverse.
double G_matrix_get_element(mat_struct *mt, int rowval, int colval)
Retrieve value of the (i,j)th element.
mat_struct * G_matrix_scale(mat_struct *mt1, const double c)
Scale a matrix by a scalar value.
int G_matrix_stdin(mat_struct *out)
mat_struct * G_matrix_resize(mat_struct *in, int rows, int cols)
Resizes a matrix.
void G_matrix_print(mat_struct *mt)
Print out a matrix.
mat_struct * G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
Calculates the scalar-matrix multiplication.
vec_struct * G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
Calculates the matrix-vector product.
mat_struct * G_matrix_add(mat_struct *mt1, mat_struct *mt2)
Adds two matricies.
int G_matvect_extract_vector(mat_struct *mt, vtype vt, int indx)
Convert matrix to vector.
int G_matrix_LU_solve(const mat_struct *mt1, mat_struct **xmat0, const mat_struct *bmat, mat_type mtype)
Solve a general system A.X = B.
mat_struct * G_matrix_init(int rows, int cols, int ldim)
Initialize a matrix structure.
double G_vector_norm1(vec_struct *vc)
Calculates the 1-norm of a vector.
int G_matrix_read(FILE *fp, mat_struct *out)
Read a matrix from a file stream.
vec_struct * G_matvect_get_row(mat_struct *mt, int row)
Retrieve a row of the matrix to a vector structure.
int G_matrix_eigen_sort(vec_struct *d, mat_struct *m)
Sort eigenvectors according to eigenvalues.
vec_struct * G_vector_copy(const vec_struct *vc1, int comp_flag)
Returns a vector copied from vc1. Underlying structure is preserved unless DO_COMPACT flag.
void G_vector_free(vec_struct *v)
Free an allocated vector structure.
void G_matrix_free(mat_struct *mt)
Free up allocated matrix.
mat_struct * G_matrix_product(mat_struct *mt1, mat_struct *mt2)
Returns product of two matricies.
int G_vector_set(vec_struct *A, int cells, int ldim, vtype vt, int vindx)
mat_struct * G_matrix_transpose(mat_struct *mt)
Transpose a matrix.
int G_matrix_set_element(mat_struct *mt, int rowval, int colval, double val)
int G_matrix_set(mat_struct *A, int rows, int cols, int ldim)
Set parameters for an initialized matrix.
int G_matrix_zero(mat_struct *A)
Clears (or resets) the matrix values to 0.