32 #if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS)
38 static int egcmp(
const void *pa,
const void *pb);
57 if (rows < 1 || cols < 1 || ldim < rows) {
58 G_warning(
_(
"Matrix dimensions out of range"));
63 tmp_arry->
rows = rows;
64 tmp_arry->
cols = cols;
65 tmp_arry->
ldim = ldim;
111 if (rows < 1 || cols < 1 || ldim < 0) {
112 G_warning(
_(
"Matrix dimensions out of range"));
144 G_warning(
_(
"Matrix is not initialised fully."));
149 G_warning(
_(
"Unable to allocate space for matrix copy"));
212 if (matrix ==
NULL) {
213 G_warning(
_(
"Input matrix is uninitialized"));
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"));
292 G_warning(
_(
"One or both input matrices uninitialised"));
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++) {
319 for (i = 0; i < mt3->
rows; i++) {
320 for (j = 0; j < mt3->
cols; j++) {
331 #if defined(HAVE_LIBBLAS)
350 integer rows, cols, interdim, lda, ldb;
354 G_warning(
_(
"One or both input matrices uninitialised"));
364 G_warning(
_(
"Unable to allocate space for matrix product"));
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
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)
474 G_warning(
_(
"Input: one or both data matrices uninitialised"));
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"));
513 integer num_eqns, nrhs, lda, ldb;
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)
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++) {
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"));
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"));
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"));
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;
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++) {
976 if ((cells < 1) || (vt ==
RVEC && ldim < 1) ||
977 (vt ==
CVEC && ldim < cells) || ldim < 0) {
978 G_warning(
"Vector dimensions out of range.");
986 tmp_arry->
cols = cells;
987 tmp_arry->
ldim = ldim;
991 else if (vt ==
CVEC) {
992 tmp_arry->
rows = cells;
994 tmp_arry->
ldim = ldim;
1043 int idx1, idx2, idx0;
1047 G_warning(
_(
"Output vector is uninitialized"));
1052 G_warning(
_(
"Vectors are not of the same type"));
1057 G_warning(
_(
"Output vector is of incorrect type"));
1068 G_warning(
_(
"Vectors have differing dimensions"));
1074 G_warning(
_(
"Output vector has incorrect dimension"));
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)) {
1153 #if defined(HAVE_LIBBLAS)
1193 return (
double)
f77_dnrm2(&Nval, startpt, &incr);
1197 #warning G_vector_norm_euclid() not compiled; requires BLAS library
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;
1302 for (i = 0; i < vc->
cols; i++)
1306 for (i = 0; i < vc->
rows; i++)
1328 int idx1, idx2, idx0;
1332 G_warning(
_(
"Output vector is uninitialized"));
1337 G_warning(
_(
"Vectors are not of the same type"));
1342 G_warning(
_(
"Output vector is not the same type as others"));
1353 G_warning(
_(
"Vectors have differing dimensions"));
1359 G_warning(
_(
"Output vector has incorrect dimension"));
1363 #if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS)
1371 for (i = 0; i < v1->
cols; i++)
1377 for (i = 0; i < v1->
rows; i++)
1402 doublereal *startpt1, *startpt2, *curpt1, *curpt2;
1406 G_warning(
_(
"Vector structure is not initialised"));
1440 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
1448 startpt1 = tmp_arry->
vals;
1457 startpt1 = tmp_arry->
vals;
1466 G_warning(
"Structure type is not vector.");
1471 startpt1 = tmp_arry->
vals;
1472 startpt2 = vc1->
vals;
1480 G_warning(
"Copy method must be specified: [DO,NO]_COMPACT.\n");
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++)
1626 for (i = 0; i < m->
cols; i++) {
1627 for (j = 0; j < m->
rows; j++)
1639 for (i = 0; i < m->
cols; i++) {
1640 for (j = 0; j < m->
rows; j++)
1653 static int egcmp(
const void *pa,
const void *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.
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.
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.