33 #if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS) 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"));
66 tmp_arry->
rows = rows;
67 tmp_arry->
cols = cols;
68 tmp_arry->
ldim = ldim;
116 if (rows < 1 || cols < 1 || ldim < 0) {
117 G_warning(
_(
"Matrix dimensions out of range"));
150 G_warning(
_(
"Matrix is not initialised fully."));
155 G_warning(
_(
"Unable to allocate space for matrix copy"));
219 if (matrix ==
NULL) {
220 G_warning (
_(
"Input matrix is uninitialized"));
233 for (i = 0; i < m; i++) {
234 for (j = 0; j < n; j++) {
285 G_warning(
_(
"First scalar multiplier must be non-zero"));
291 G_warning(
_(
"One or both input matrices uninitialised"));
299 G_warning(
_(
"One or both input matrices uninitialised"));
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++) {
326 for (i = 0; i < mt3->
rows; i++) {
327 for (j = 0; j < mt3->
cols; j++) {
341 #if defined(HAVE_LIBBLAS) 360 integer rows, cols, interdim, lda, ldb;
364 G_warning(
_(
"One or both input matrices uninitialised"));
374 G_warning(
_(
"Unable to allocate space for matrix product"));
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 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) 486 G_warning(
_(
"Input: one or both data matrices uninitialised"));
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"));
526 integer num_eqns, nrhs, lda, ldb;
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) 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++) {
627 for (j = i + 1; j < mt0->
cols; j++) {
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"));
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"));
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"));
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;
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++) {
999 if ((cells < 1) || (vt ==
RVEC && ldim < 1)
1000 || (vt ==
CVEC && ldim < cells) || ldim < 0) {
1001 G_warning(
"Vector dimensions out of range.");
1009 tmp_arry->
cols = cells;
1010 tmp_arry->
ldim = ldim;
1014 else if (vt ==
CVEC) {
1015 tmp_arry->
rows = cells;
1017 tmp_arry->
ldim = ldim;
1067 int idx1, idx2, idx0;
1071 G_warning(
_(
"Output vector is uninitialized"));
1076 G_warning(
_(
"Vectors are not of the same type"));
1081 G_warning(
_(
"Output vector is of incorrect type"));
1092 G_warning(
_(
"Vectors have differing dimensions"));
1098 G_warning(
_(
"Output vector has incorrect dimension"));
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)) {
1178 #if defined(HAVE_LIBBLAS) 1218 return (
double)
f77_dnrm2(&Nval, startpt, &incr);
1222 #warning G_vector_norm_euclid() not compiled; requires BLAS library 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;
1332 for (i = 0; i < vc->
cols; i++)
1336 for (i = 0; i < vc->
rows; i++)
1357 int idx1, idx2, idx0;
1361 G_warning (
_(
"Output vector is uninitialized"));
1366 G_warning (
_(
"Vectors are not of the same type"));
1371 G_warning (
_(
"Output vector is not the same type as others"));
1383 G_warning (
_(
"Vectors have differing dimensions"));
1390 G_warning (
_(
"Output vector has incorrect dimension"));
1394 #if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS) 1402 for (i = 0; i < v1->
cols; i++)
1407 for (i = 0; i < v1->
rows; i++)
1433 doublereal *startpt1, *startpt2, *curpt1, *curpt2;
1437 G_warning(
_(
"Vector structure is not initialised"));
1471 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1479 startpt1 = tmp_arry->
vals;
1488 startpt1 = tmp_arry->
vals;
1497 G_warning(
"Structure type is not vector.");
1502 startpt1 = tmp_arry->
vals;
1503 startpt2 = vc1->
vals;
1511 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
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++)
1659 for (i = 0; i < m->
cols; i++) {
1660 for (j = 0; j < m->
rows; j++)
1673 for (i = 0; i < m->
cols; i++) {
1674 for (j = 0; j < m->
rows; j++)
1689 static int egcmp(
const void *pa,
const void *pb)
mat_struct * G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
Calculates the scalar-matrix multiplication.
int G_getl(char *, int, FILE *)
Gets a line of text from a file.
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
mat_struct * G_matrix_scale(mat_struct *mt1, const double c)
Scale a matrix by a scalar value.
vec_struct * G_vector_init(int cells, int ldim, vtype vt)
Initialize a vector structure.
void G_matrix_free(mat_struct *mt)
Free up allocated matrix.
int G_matvect_extract_vector(mat_struct *mt, vtype vt, int indx)
Convert matrix to vector.
mat_struct * G_matrix_transpose(mat_struct *mt)
Transpose a matrix.
int G_matrix_set(mat_struct *A, int rows, int cols, int ldim)
Set parameters for an initialized matrix.
void G_free(void *)
Free allocated memory.
vec_struct * G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
Calculates the matrix-vector product.
mat_struct * G_matrix_init(int rows, int cols, int ldim)
Initialize a matrix structure.
vec_struct * G_vector_product(vec_struct *v1, vec_struct *v2, vec_struct *out)
Calculates the vector product.
void G_matrix_print(mat_struct *mt)
Print out a matrix.
void G_message(const char *,...) __attribute__((format(printf
mat_struct * G_matrix_add(mat_struct *mt1, mat_struct *mt2)
Adds two matricies.
double G_vector_norm_maxval(vec_struct *vc, int vflag)
Calculates maximum value.
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...
double G_matrix_get_element(mat_struct *mt, int rowval, int colval)
Retrieve value of the (i,j)th element.
int G_vector_set(vec_struct *A, int cells, int ldim, vtype vt, int vindx)
Set parameters for vector structure.
double G_vector_norm_euclid(vec_struct *vc)
Calculates euclidean norm.
mat_struct * G__matrix_add(mat_struct *mt1, mat_struct *mt2, const double c1, const double c2)
General add/subtract/scalar multiply routine.
int G_matrix_eigen_sort(vec_struct *d, mat_struct *m)
Sort eigenvectors according to eigenvalues.
vec_struct * G_matvect_get_row(mat_struct *mt, int row)
Retrieve a row 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.
void G_warning(const char *,...) __attribute__((format(printf
int G_matrix_set_element(mat_struct *mt, int rowval, int colval, double val)
Set the value of the (i, j)th element.
void G_vector_free(vec_struct *v)
Free an allocated vector structure.
int G_matrix_read(FILE *fp, mat_struct *out)
Read a matrix from a file stream.
mat_struct * G_matrix_subtract(mat_struct *mt1, mat_struct *mt2)
Subtract two matricies.
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_resize(mat_struct *in, int rows, int cols)
Resizes a matrix.
int G_matrix_zero(mat_struct *A)
Clears (or resets) the matrix values to 0.
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_product(mat_struct *mt1, mat_struct *mt2)
Returns product of two matricies.
double G_vector_norm1(vec_struct *vc)
Calculates the 1-norm of a vector.
int G_matrix_stdin(mat_struct *out)
int G_matvect_retrieve_matrix(vec_struct *vc)
Revert a vector to matrix.
mat_struct * G_matrix_inverse(mat_struct *mt)
Returns the matrix inverse.