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);
71 if (rows < 1 || cols < 1 || ldim < rows) {
72 G_warning(
_(
"Matrix dimensions out of range"));
77 tmp_arry->
rows = rows;
78 tmp_arry->
cols = cols;
79 tmp_arry->
ldim = ldim;
83 tmp_arry->
vals = (
double *)
G_calloc(ldim * cols,
sizeof(
double));
103 memset(A->
vals, 0, (A->
ldim * A->
cols) *
sizeof(
double));
125 if (rows < 1 || cols < 1 || ldim < 0) {
126 G_warning(
_(
"Matrix dimensions out of range"));
136 A->
vals = (
double *)
G_calloc(ldim * cols,
sizeof(
double));
158 G_warning(
_(
"Matrix is not initialised fully."));
163 G_warning(
_(
"Unable to allocate space for matrix copy"));
168 (
size_t)A->
cols * A->
ldim *
sizeof(
double));
227 if (matrix ==
NULL) {
228 G_warning(
_(
"Input matrix is uninitialized"));
241 for (i = 0; i < m; i++) {
242 for (j = 0; j < n; j++) {
293 G_warning(
_(
"First scalar multiplier must be non-zero"));
299 G_warning(
_(
"One or both input matrices uninitialised"));
307 G_warning(
_(
"One or both input matrices uninitialised"));
318 G_warning(
_(
"Unable to allocate space for matrix sum"));
324 for (i = 0; i < mt3->
rows; i++) {
325 for (j = 0; j < mt3->
cols; j++) {
334 for (i = 0; i < mt3->
rows; i++) {
335 for (j = 0; j < mt3->
cols; j++) {
362 double unity = 1., zero = 0.;
363 int rows, cols, interdim, lda, ldb;
366 G_warning(
_(
"One or both input matrices uninitialised"));
376 G_warning(
_(
"Unable to allocate space for matrix product"));
382 rows = (int)mt1->
rows;
383 interdim = (
int)mt1->
cols;
384 cols = (int)mt2->
cols;
386 lda = (
int)mt1->
ldim;
387 ldb = (int)mt2->
ldim;
389 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, rows, cols, interdim,
390 unity, mt1->
vals, lda, mt2->
vals, ldb, zero, mt3->
vals, lda);
412 double *dbo, *dbt, *dbx, *dby;
416 if (mt->
cols % 2 == 0)
428 for (cnt = 0; cnt < mt->
cols; cnt++) {
432 for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
440 if (cnt < mt->cols - 1) {
480 G_warning(
_(
"Input: one or both data matrices uninitialised"));
485 G_warning(
_(
"Principal matrix is not properly dimensioned"));
489 if (bmat->
cols < 1) {
490 G_warning(
_(
"Input: you must have at least one array to solve"));
496 G_warning(
_(
"Could not allocate space for solution matrix"));
502 G_warning(
_(
"Could not allocate space for working matrix"));
510 G_warning(
_(
"Could not allocate space for working matrix"));
519 int num_eqns, nrhs, lda, ldb;
524 num_eqns = (int)mt1->
rows;
525 nrhs = (
int)wmat->
cols;
526 lda = (int)mt1->
ldim;
527 ldb = (
int)wmat->
ldim;
530 res_info = LAPACKE_dgesv(LAPACK_COL_MAJOR, num_eqns, nrhs, mtx->
vals,
531 lda, perm, wmat->
vals, ldb);
554 (
size_t)wmat->
cols * wmat->
ldim *
sizeof(
double));
563 _(
"Matrix (or submatrix is singular). Solution undetermined"));
566 else if (res_info < 0) {
573 G_warning(
_(
"Procedure not yet available for selected matrix type"));
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++) {
614 for (j = i + 1; j < mt0->
cols; j++) {
671 char buf[2048], numbuf[64];
673 for (i = 0; i < mt->
rows; i++) {
676 for (j = 0; j < mt->
cols; j++) {
677 snprintf(numbuf,
sizeof(numbuf),
"%14.6f",
680 if (j < mt->cols - 1)
687 fprintf(stderr,
"\n");
710 G_warning(
_(
"Element array has not been allocated"));
714 if (rowval >= mt->
rows || colval >= mt->
cols || rowval < 0 || colval < 0) {
715 G_warning(
_(
"Specified element is outside array bounds"));
719 mt->
vals[rowval + colval * mt->
ldim] = (double)val;
746 return (val = (
double)mt->
vals[rowval + colval * mt->
ldim]);
766 if (col < 0 || col >= mt->
cols) {
767 G_warning(
_(
"Specified matrix column index is outside range"));
777 G_warning(
_(
"Could not allocate space for vector structure"));
781 for (i = 0; i < mt->
rows; i++)
806 if (row < 0 || row >= mt->
cols) {
807 G_warning(
_(
"Specified matrix row index is outside range"));
817 G_warning(
_(
"Could not allocate space for vector structure"));
821 for (i = 0; i < mt->
cols; i++)
845 if (vt ==
RVEC && indx >= mt->
rows) {
846 G_warning(
_(
"Specified row index is outside range"));
850 else if (vt ==
CVEC && indx >= mt->
cols) {
851 G_warning(
_(
"Specified column index is outside range"));
916 unsigned int i, m, n, j;
921 if (A->
cols !=
b->cols) {
922 G_warning(
_(
"Input matrix and vector have differing dimensions1"));
927 G_warning(
_(
"Output vector is uninitialized"));
939 for (i = 0; i < m; i++) {
943 for (j = 0; j < n; j++) {
972 if ((cells < 1) || (vt ==
RVEC && ldim < 1) ||
973 (vt ==
CVEC && ldim < cells) || ldim < 0) {
974 G_warning(
"Vector dimensions out of range.");
982 tmp_arry->
cols = cells;
983 tmp_arry->
ldim = ldim;
987 else if (vt ==
CVEC) {
988 tmp_arry->
rows = cells;
990 tmp_arry->
ldim = ldim;
1038 int idx1, idx2, idx0;
1042 G_warning(
_(
"Output vector is uninitialized"));
1047 G_warning(
_(
"Vectors are not of the same type"));
1052 G_warning(
_(
"Output vector is of incorrect type"));
1063 G_warning(
_(
"Vectors have differing dimensions"));
1069 G_warning(
_(
"Output vector has incorrect dimension"));
1078 for (i = 0; i < v1->
cols; i++)
1084 for (i = 0; i < v1->
rows; i++)
1113 if ((cells < 1) || (vt ==
RVEC && ldim < 1) ||
1114 (vt ==
CVEC && ldim < cells) || ldim < 0) {
1115 G_warning(
_(
"Vector dimensions out of range"));
1119 if ((vt ==
RVEC && vindx >= A->
cols) || (vt ==
CVEC && vindx >= A->
rows)) {
1169 Nval = (int)vc->
cols;
1170 incr = (
int)vc->
ldim;
1177 Nval = (int)vc->
rows;
1186 return cblas_dnrm2(Nval, startpt, incr);
1208 double xval, *startpt, *curpt;
1216 ncells = (int)vc->
cols;
1217 incr = (
int)vc->
ldim;
1224 ncells = (int)vc->
rows;
1235 while (ncells > 0) {
1236 if (curpt != startpt) {
1252 cellval = (double)(*curpt);
1253 if (hypot(cellval, cellval) > (double)xval)
1263 return (
double)xval;
1279 double result = 0.0;
1291 for (i = 0; i < vc->
cols; i++)
1295 for (i = 0; i < vc->
rows; i++)
1317 int idx1, idx2, idx0;
1321 G_warning(
_(
"Output vector is uninitialized"));
1326 G_warning(
_(
"Vectors are not of the same type"));
1331 G_warning(
_(
"Output vector is not the same type as others"));
1342 G_warning(
_(
"Vectors have differing dimensions"));
1348 G_warning(
_(
"Output vector has incorrect dimension"));
1357 for (i = 0; i < v1->
cols; i++)
1363 for (i = 0; i < v1->
rows; i++)
1387 double *startpt1, *startpt2, *curpt1, *curpt2;
1391 G_warning(
_(
"Vector structure is not initialised"));
1425 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1433 startpt1 = tmp_arry->
vals;
1442 startpt1 = tmp_arry->
vals;
1451 G_warning(
"Structure type is not vector.");
1456 startpt1 = tmp_arry->
vals;
1457 startpt2 = vc1->
vals;
1465 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1470 memcpy(curpt1, curpt2,
sizeof(
double));
1504 if (!
G_getl(buff,
sizeof(buff), fp))
1510 if (sscanf(buff,
"Matrix: %d by %d", &rows, &cols) != 2) {
1517 for (i = 0; i < rows; i++) {
1518 if (fscanf(fp,
"row%d:", &row) != 1 || row != i) {
1522 for (j = 0; j < cols; j++) {
1523 if (fscanf(fp,
"%lf:", &val) != 1) {
1555 for (i = 0; i < rows; i++)
1556 for (j = 0; j < cols; j++)
1560 int old_size = in->
rows * in->
cols;
1561 int new_size = rows * cols;
1563 if (new_size > old_size)
1564 for (p = old_size; p < new_size; p++)
1611 for (i = 0; i < m->
cols; i++) {
1612 for (j = 0; j < m->
rows; j++)
1621 qsort(tmp.
vals, tmp.
cols, tmp.
ldim *
sizeof(
double), egcmp);
1624 for (i = 0; i < m->
cols; i++) {
1625 for (j = 0; j < m->
rows; j++)
1638 static int egcmp(
const void *pa,
const void *pb)
1640 double a = *(
double *
const)pa;
1641 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 matricies.
mat_struct * G_matrix_add(mat_struct *mt1, mat_struct *mt2)
Adds two matricies.
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)
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 matricies.
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.