27 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
36 #define LAPACK_COMPLEX_CUSTOM
37 #define lapack_complex_float _Fcomplex
38 #define lapack_complex_double _Dcomplex
42 #if defined(HAVE_CBLAS_ATLAS_H)
43 #include <cblas-atlas.h>
52 static int egcmp(
const void *pa,
const void *pb);
70 if (rows < 1 || cols < 1 || ldim < rows) {
71 G_warning(
_(
"Matrix dimensions out of range"));
76 tmp_arry->
rows = rows;
77 tmp_arry->
cols = cols;
78 tmp_arry->
ldim = ldim;
82 tmp_arry->
vals = (
double *)
G_calloc(ldim * cols,
sizeof(
double));
101 memset(A->
vals, 0, (A->
ldim * A->
cols) *
sizeof(
double));
122 if (rows < 1 || cols < 1 || ldim < 0) {
123 G_warning(
_(
"Matrix dimensions out of range"));
133 A->
vals = (
double *)
G_calloc(ldim * cols,
sizeof(
double));
154 G_warning(
_(
"Matrix is not initialised fully."));
159 G_warning(
_(
"Unable to allocate space for matrix copy"));
164 (
size_t)A->
cols * A->
ldim *
sizeof(
double));
221 if (matrix ==
NULL) {
222 G_warning(
_(
"Input matrix is uninitialized"));
235 for (i = 0; i < m; i++) {
236 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++) {
353 double unity = 1., zero = 0.;
354 int rows, cols, interdim, lda, ldb;
357 G_warning(
_(
"One or both input matrices uninitialised"));
367 G_warning(
_(
"Unable to allocate space for matrix product"));
373 rows = (int)mt1->
rows;
374 interdim = (
int)mt1->
cols;
375 cols = (int)mt2->
cols;
377 lda = (
int)mt1->
ldim;
378 ldb = (int)mt2->
ldim;
380 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, rows, cols, interdim,
381 unity, mt1->
vals, lda, mt2->
vals, ldb, zero, mt3->
vals, lda);
402 double *dbo, *dbt, *dbx, *dby;
406 if (mt->
cols % 2 == 0)
418 for (cnt = 0; cnt < mt->
cols; cnt++) {
422 for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
430 if (cnt < mt->cols - 1) {
470 G_warning(
_(
"Input: one or both data matrices uninitialised"));
475 G_warning(
_(
"Principal matrix is not properly dimensioned"));
479 if (bmat->
cols < 1) {
480 G_warning(
_(
"Input: you must have at least one array to solve"));
486 G_warning(
_(
"Could not allocate space for solution matrix"));
492 G_warning(
_(
"Could not allocate space for working matrix"));
500 G_warning(
_(
"Could not allocate space for working matrix"));
509 int num_eqns, nrhs, lda, ldb;
514 num_eqns = (int)mt1->
rows;
515 nrhs = (
int)wmat->
cols;
516 lda = (int)mt1->
ldim;
517 ldb = (
int)wmat->
ldim;
520 res_info = LAPACKE_dgesv(LAPACK_COL_MAJOR, num_eqns, nrhs, mtx->
vals,
521 lda, perm, wmat->
vals, ldb);
544 (
size_t)wmat->
cols * wmat->
ldim *
sizeof(
double));
553 _(
"Matrix (or submatrix is singular). Solution undetermined"));
556 else if (res_info < 0) {
563 G_warning(
_(
"Procedure not yet available for selected matrix type"));
590 G_warning(
_(
"Matrix is not square. Cannot determine inverse"));
595 G_warning(
_(
"Unable to allocate space for matrix"));
600 for (i = 0; i < mt0->
rows - 1; i++) {
603 for (j = i + 1; j < mt0->
cols; j++) {
658 char buf[2048], numbuf[64];
660 for (i = 0; i < mt->
rows; i++) {
663 for (j = 0; j < mt->
cols; j++) {
664 snprintf(numbuf,
sizeof(numbuf),
"%14.6f",
667 if (j < mt->cols - 1)
674 fprintf(stderr,
"\n");
696 G_warning(
_(
"Element array has not been allocated"));
700 if (rowval >= mt->
rows || colval >= mt->
cols || rowval < 0 || colval < 0) {
701 G_warning(
_(
"Specified element is outside array bounds"));
705 mt->
vals[rowval + colval * mt->
ldim] = (double)val;
731 return (val = (
double)mt->
vals[rowval + colval * mt->
ldim]);
750 if (col < 0 || col >= mt->
cols) {
751 G_warning(
_(
"Specified matrix column index is outside range"));
761 G_warning(
_(
"Could not allocate space for vector structure"));
765 for (i = 0; i < mt->
rows; i++)
789 if (row < 0 || row >= mt->
cols) {
790 G_warning(
_(
"Specified matrix row index is outside range"));
800 G_warning(
_(
"Could not allocate space for vector structure"));
804 for (i = 0; i < mt->
cols; i++)
827 if (vt ==
RVEC && indx >= mt->
rows) {
828 G_warning(
_(
"Specified row index is outside range"));
832 else if (vt ==
CVEC && indx >= mt->
cols) {
833 G_warning(
_(
"Specified column index is outside range"));
896 unsigned int i, m, n, j;
901 if (A->
cols !=
b->cols) {
902 G_warning(
_(
"Input matrix and vector have differing dimensions1"));
907 G_warning(
_(
"Output vector is uninitialized"));
919 for (i = 0; i < m; i++) {
923 for (j = 0; j < n; j++) {
951 if ((cells < 1) || (vt ==
RVEC && ldim < 1) ||
952 (vt ==
CVEC && ldim < cells) || ldim < 0) {
953 G_warning(
"Vector dimensions out of range.");
961 tmp_arry->
cols = cells;
962 tmp_arry->
ldim = ldim;
966 else if (vt ==
CVEC) {
967 tmp_arry->
rows = cells;
969 tmp_arry->
ldim = ldim;
1015 int idx1, idx2, idx0;
1019 G_warning(
_(
"Output vector is uninitialized"));
1024 G_warning(
_(
"Vectors are not of the same type"));
1029 G_warning(
_(
"Output vector is of incorrect type"));
1040 G_warning(
_(
"Vectors have differing dimensions"));
1046 G_warning(
_(
"Output vector has incorrect dimension"));
1055 for (i = 0; i < v1->
cols; i++)
1061 for (i = 0; i < v1->
rows; i++)
1089 if ((cells < 1) || (vt ==
RVEC && ldim < 1) ||
1090 (vt ==
CVEC && ldim < cells) || ldim < 0) {
1091 G_warning(
_(
"Vector dimensions out of range"));
1095 if ((vt ==
RVEC && vindx >= A->
cols) || (vt ==
CVEC && vindx >= A->
rows)) {
1144 Nval = (int)vc->
cols;
1145 incr = (
int)vc->
ldim;
1152 Nval = (int)vc->
rows;
1161 return cblas_dnrm2(Nval, startpt, incr);
1182 double xval, *startpt, *curpt;
1190 ncells = (int)vc->
cols;
1191 incr = (
int)vc->
ldim;
1198 ncells = (int)vc->
rows;
1209 while (ncells > 0) {
1210 if (curpt != startpt) {
1226 cellval = (double)(*curpt);
1227 if (hypot(cellval, cellval) > (double)xval)
1237 return (
double)xval;
1252 double result = 0.0;
1264 for (i = 0; i < vc->
cols; i++)
1268 for (i = 0; i < vc->
rows; i++)
1290 int idx1, idx2, idx0;
1294 G_warning(
_(
"Output vector is uninitialized"));
1299 G_warning(
_(
"Vectors are not of the same type"));
1304 G_warning(
_(
"Output vector is not the same type as others"));
1315 G_warning(
_(
"Vectors have differing dimensions"));
1321 G_warning(
_(
"Output vector has incorrect dimension"));
1330 for (i = 0; i < v1->
cols; i++)
1336 for (i = 0; i < v1->
rows; i++)
1359 double *startpt1, *startpt2, *curpt1, *curpt2;
1363 G_warning(
_(
"Vector structure is not initialised"));
1397 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1405 startpt1 = tmp_arry->
vals;
1414 startpt1 = tmp_arry->
vals;
1423 G_warning(
"Structure type is not vector.");
1428 startpt1 = tmp_arry->
vals;
1429 startpt2 = vc1->
vals;
1437 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1442 memcpy(curpt1, curpt2,
sizeof(
double));
1475 if (!
G_getl(buff,
sizeof(buff), fp))
1481 if (sscanf(buff,
"Matrix: %d by %d", &rows, &cols) != 2) {
1488 for (i = 0; i < rows; i++) {
1489 if (fscanf(fp,
"row%d:", &row) != 1 || row != i) {
1493 for (j = 0; j < cols; j++) {
1494 if (fscanf(fp,
"%lf:", &val) != 1) {
1525 for (i = 0; i < rows; i++)
1526 for (j = 0; j < cols; j++)
1530 int old_size = in->
rows * in->
cols;
1531 int new_size = rows * cols;
1533 if (new_size > old_size)
1534 for (p = old_size; p < new_size; p++)
1579 for (i = 0; i < m->
cols; i++) {
1580 for (j = 0; j < m->
rows; j++)
1589 qsort(tmp.
vals, tmp.
cols, tmp.
ldim *
sizeof(
double), egcmp);
1592 for (i = 0; i < m->
cols; i++) {
1593 for (j = 0; j < m->
rows; j++)
1606 static int egcmp(
const void *pa,
const void *pb)
1608 double a = *(
double *
const)pa;
1609 double b = *(
double *
const)pb;
void G_free(void *)
Free allocated memory.
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_warning(const char *,...) __attribute__((format(printf
void G_message(const char *,...) __attribute__((format(printf
int G_getl(char *, int, FILE *)
Gets a line of text from a file.
size_t G_strlcpy(char *, const char *, size_t)
Safe string copy function.
vec_struct * G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
Calculates the matrix-vector product.
double G_vector_norm_maxval(vec_struct *vc, int vflag)
Calculates maximum value.
vec_struct * G_matvect_get_column(mat_struct *mt, int col)
Retrieve a column of the matrix to a vector structure.
vec_struct * G_vector_sub(vec_struct *v1, vec_struct *v2, vec_struct *out)
Subtract two vectors.
int G_matvect_retrieve_matrix(vec_struct *vc)
Revert a vector to matrix.
double G_vector_norm_euclid(vec_struct *vc)
Calculates euclidean norm.
mat_struct * G_matrix_subtract(mat_struct *mt1, mat_struct *mt2)
Subtract two matrices.
mat_struct * G_matrix_add(mat_struct *mt1, mat_struct *mt2)
Adds two matrices.
double G_matrix_get_element(mat_struct *mt, int rowval, int colval)
Retrieve value of the (i,j)th element.
mat_struct * G__matrix_add(mat_struct *mt1, mat_struct *mt2, const double c1, const double c2)
mat_struct * G_matrix_inverse(mat_struct *mt)
Returns the matrix inverse.
int G_matrix_stdin(mat_struct *out)
Read a matrix from standard input.
void G_matrix_print(mat_struct *mt)
Print out a matrix.
mat_struct * G_matrix_copy(const mat_struct *A)
Copy a matrix.
mat_struct * G_matrix_resize(mat_struct *in, int rows, int cols)
Resizes a matrix.
int G_matvect_extract_vector(mat_struct *mt, vtype vt, int indx)
Convert matrix to vector.
vec_struct * G_vector_init(int cells, int ldim, vtype vt)
Initialize a vector structure.
mat_struct * G_matrix_product(mat_struct *mt1, mat_struct *mt2)
Returns product of two matrices.
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.
vec_struct * G_vector_product(vec_struct *v1, vec_struct *v2, vec_struct *out)
Calculates the vector product.
mat_struct * G_matrix_init(int rows, int cols, int ldim)
Initialize a matrix structure.
int suppress_empty_translation_unit_compiler_warning
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.
int G_matrix_eigen_sort(vec_struct *d, mat_struct *m)
Sort eigenvectors according to eigenvalues.
mat_struct * G_matrix_scale(mat_struct *mt1, const double c)
Scale a matrix by a scalar value.
void G_vector_free(vec_struct *v)
Free an allocated vector structure.
vec_struct * G_matvect_get_row(mat_struct *mt, int row)
Retrieve a row of the matrix to a vector structure.
void G_matrix_free(mat_struct *mt)
Free up allocated matrix.
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.
int G_vector_set(vec_struct *A, int cells, int ldim, vtype vt, int vindx)
int G_matrix_set_element(mat_struct *mt, int rowval, int colval, double val)
mat_struct * G_matrix_transpose(mat_struct *mt)
Transpose a matrix.
mat_struct * G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
Calculates the scalar-matrix multiplication.
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.
Wrapper headers for BLAS/LAPACK.