GRASS GIS 8 Programmer's Manual 8.2.1RC1(2022)-exported
la.c
Go to the documentation of this file.
1
2/******************************************************************************
3 * la.c
4 * wrapper modules for linear algebra problems
5 * linking to BLAS / LAPACK (and others?)
6
7 * @Copyright David D.Gray <ddgray@armadce.demon.co.uk>
8 * 26th. Sep. 2000
9 * Last updated:
10 * 2006-11-23
11 * 2015-01-20
12
13 * This file is part of GRASS GIS. It is free software. You can
14 * redistribute it and/or modify it under the terms of
15 * the GNU General Public License as published by the Free Software
16 * Foundation; either version 2 of the License, or (at your option)
17 * any later version.
18
19 * This program is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU General Public License for more details.
23
24 ******************************************************************************/
25
26#include <stdio.h> /* needed here for ifdef/else */
27#include <stdlib.h>
28#include <string.h>
29#include <math.h>
30
31#include <grass/config.h>
32
33#if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS)
34
35#include <grass/gis.h>
36#include <grass/glocale.h>
37#include <grass/la.h>
38
39
40static int egcmp(const void *pa, const void *pb);
41
42
43/*!
44 * \fn mat_struct *G_matrix_init(int rows, int cols, int ldim)
45 *
46 * \brief Initialize a matrix structure
47 *
48 * Initialize a matrix structure
49 *
50 * \param rows
51 * \param cols
52 * \param ldim
53 * \return mat_struct
54 */
55
56mat_struct *G_matrix_init(int rows, int cols, int ldim)
57{
58 mat_struct *tmp_arry;
59
60 if (rows < 1 || cols < 1 || ldim < rows) {
61 G_warning(_("Matrix dimensions out of range"));
62 return NULL;
63 }
64
65 tmp_arry = (mat_struct *) G_malloc(sizeof(mat_struct));
66 tmp_arry->rows = rows;
67 tmp_arry->cols = cols;
68 tmp_arry->ldim = ldim;
69 tmp_arry->type = MATRIX_;
70 tmp_arry->v_indx = -1;
71
72 tmp_arry->vals = (doublereal *) G_calloc(ldim * cols, sizeof(doublereal));
73 tmp_arry->is_init = 1;
74
75 return tmp_arry;
76}
77
78
79/*!
80 * \fn int G_matrix_zero (mat_struct *A)
81 *
82 * \brief Clears (or resets) the matrix values to 0
83 *
84 * \param A
85 * \return 0 on error; 1 on success
86 */
87
88int G_matrix_zero(mat_struct * A)
89{
90 if (!A->vals)
91 return 0;
92
93 memset(A->vals, 0, (A->ldim * A->cols) * sizeof(doublereal));
94
95 return 1;
96}
97
98
99/*!
100 * \fn int G_matrix_set(mat_struct *A, int rows, int cols, int ldim)
101 *
102 * \brief Set parameters for an initialized matrix
103 *
104 * Set parameters for matrix <b>A</b> that is allocated,
105 * but not yet fully initialized. Is an alternative to G_matrix_init().
106 *
107 * \param A
108 * \param rows
109 * \param cols
110 * \param ldim
111 * \return int
112 */
113
114int G_matrix_set(mat_struct * A, int rows, int cols, int ldim)
115{
116 if (rows < 1 || cols < 1 || ldim < 0) {
117 G_warning(_("Matrix dimensions out of range"));
118 return -1;
119 }
120
121 A->rows = rows;
122 A->cols = cols;
123 A->ldim = ldim;
124 A->type = MATRIX_;
125 A->v_indx = -1;
126
127 A->vals = (doublereal *) G_calloc(ldim * cols, sizeof(doublereal));
128 A->is_init = 1;
129
130 return 0;
131}
132
133
134/*!
135 * \fn mat_struct *G_matrix_copy (const mat_struct *A)
136 *
137 * \brief Copy a matrix
138 *
139 * Copy matrix <b>A</b> by exactly duplicating its contents.
140 *
141 * \param A
142 * \return mat_struct
143 */
144
145mat_struct *G_matrix_copy(const mat_struct * A)
146{
147 mat_struct *B;
148
149 if (!A->is_init) {
150 G_warning(_("Matrix is not initialised fully."));
151 return NULL;
152 }
153
154 if ((B = G_matrix_init(A->rows, A->cols, A->ldim)) == NULL) {
155 G_warning(_("Unable to allocate space for matrix copy"));
156 return NULL;
157 }
158
159 memcpy(&B->vals[0], &A->vals[0], A->cols * A->ldim * sizeof(doublereal));
160
161 return B;
162}
163
164
165/*!
166 * \fn mat_struct *G_matrix_add (mat_struct *mt1, mat_struct *mt2)
167 *
168 * \brief Adds two matricies
169 *
170 * Adds two matricies <b>mt1</b> and <b>mt2</b> and returns a
171 * resulting matrix. The return structure is automatically initialized.
172 *
173 * \param mt1
174 * \param mt2
175 * \return mat_struct
176 */
177
178mat_struct *G_matrix_add(mat_struct * mt1, mat_struct * mt2)
179{
180 return G__matrix_add(mt1, mt2, 1, 1);
181}
182
183
184/*!
185 * \fn mat_struct *G_matrix_subtract (mat_struct *mt1, mat_struct *mt2)
186 *
187 * \brief Subtract two matricies
188 *
189 * Subtracts two matricies <b>mt1</b> and <b>mt2</b> and returns
190 * a resulting matrix. The return matrix is automatically initialized.
191 *
192 * \param mt1
193 * \param mt2
194 * \return mat_struct
195 */
196
197mat_struct *G_matrix_subtract(mat_struct * mt1, mat_struct * mt2)
198{
199 return G__matrix_add(mt1, mt2, 1, -1);
200}
201
202/*!
203 * \fn mat_struct *G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
204 *
205 * \brief Calculates the scalar-matrix multiplication
206 *
207 * Calculates the scalar-matrix multiplication
208 *
209 * \param scalar
210 * \param A
211 * \return mat_struct
212 */
213
214mat_struct *G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
215{
216 int m, n, i, j;
217 int index = 0;
218
219 if (matrix == NULL) {
220 G_warning (_("Input matrix is uninitialized"));
221 return NULL;
222 }
223
224 if (out == NULL)
225 out = G_matrix_init(matrix->rows, matrix->cols, matrix->rows);
226
227 if (out->rows != matrix->rows || out->cols != matrix->cols)
228 out = G_matrix_resize(out, matrix->rows, matrix->cols);
229
230 m = matrix->rows;
231 n = matrix->cols;
232
233 for (i = 0; i < m; i++) {
234 for (j = 0; j < n; j++) {
235 doublereal value = scalar * G_matrix_get_element(matrix, i, j);
236 G_matrix_set_element (out, i,j, value);
237 }
238 }
239
240 return (out);
241}
242
243
244/*!
245 * \fn mat_struct *G_matrix_scale (mat_struct *mt1, const double c)
246 *
247 * \brief Scale a matrix by a scalar value
248 *
249 * Scales matrix <b>mt1</b> by scalar value <b>c</b>. The
250 * resulting matrix is automatically initialized.
251 *
252 * \param mt1
253 * \param c
254 * \return mat_struct
255 */
256
257mat_struct *G_matrix_scale(mat_struct * mt1, const double c)
258{
259 return G__matrix_add(mt1, NULL, c, 0);
260}
261
262
263/*!
264 * \fn mat_struct *G__matrix_add (mat_struct *mt1, mat_struct *mt2, const double c1, const double c2)
265 *
266 * \brief General add/subtract/scalar multiply routine
267 *
268 * General add/subtract/scalar multiply routine. <b>c2</b> may be
269 * zero, but <b>c1</b> must be non-zero.
270 *
271 * \param mt1
272 * \param mt2
273 * \param c1
274 * \param c2
275 * \return mat_struct
276 */
277
278mat_struct *G__matrix_add(mat_struct * mt1, mat_struct * mt2, const double c1,
279 const double c2)
280{
281 mat_struct *mt3;
282 int i, j; /* loop variables */
283
284 if (c1 == 0) {
285 G_warning(_("First scalar multiplier must be non-zero"));
286 return NULL;
287 }
288
289 if (c2 == 0) {
290 if (!mt1->is_init) {
291 G_warning(_("One or both input matrices uninitialised"));
292 return NULL;
293 }
294 }
295
296 else {
297
298 if (!((mt1->is_init) && (mt2->is_init))) {
299 G_warning(_("One or both input matrices uninitialised"));
300 return NULL;
301 }
302
303 if (mt1->rows != mt2->rows || mt1->cols != mt2->cols) {
304 G_warning(_("Matrix order does not match"));
305 return NULL;
306 }
307 }
308
309 if ((mt3 = G_matrix_init(mt1->rows, mt1->cols, mt1->ldim)) == NULL) {
310 G_warning(_("Unable to allocate space for matrix sum"));
311 return NULL;
312 }
313
314 if (c2 == 0) {
315
316 for (i = 0; i < mt3->rows; i++) {
317 for (j = 0; j < mt3->cols; j++) {
318 mt3->vals[i + mt3->ldim * j] =
319 c1 * mt1->vals[i + mt1->ldim * j];
320 }
321 }
322 }
323
324 else {
325
326 for (i = 0; i < mt3->rows; i++) {
327 for (j = 0; j < mt3->cols; j++) {
328 mt3->vals[i + mt3->ldim * j] =
329 c1 * mt1->vals[i + mt1->ldim * j] + c2 * mt2->vals[i +
330 mt2->
331 ldim *
332 j];
333 }
334 }
335 }
336
337 return mt3;
338}
339
340
341#if defined(HAVE_LIBBLAS)
342
343/*!
344 * \fn mat_struct *G_matrix_product (mat_struct *mt1, mat_struct *mt2)
345 *
346 * \brief Returns product of two matricies
347 *
348 * Returns a matrix with the product of matrix <b>mt1</b> and
349 * <b>mt2</b>. The return matrix is automatically initialized.
350 *
351 * \param mt1
352 * \param mt2
353 * \return mat_struct
354 */
355
356mat_struct *G_matrix_product(mat_struct * mt1, mat_struct * mt2)
357{
358 mat_struct *mt3;
359 doublereal unity = 1, zero = 0;
360 integer rows, cols, interdim, lda, ldb;
361 integer1 no_trans = 'n';
362
363 if (!((mt1->is_init) || (mt2->is_init))) {
364 G_warning(_("One or both input matrices uninitialised"));
365 return NULL;
366 }
367
368 if (mt1->cols != mt2->rows) {
369 G_warning(_("Matrix order does not match"));
370 return NULL;
371 }
372
373 if ((mt3 = G_matrix_init(mt1->rows, mt2->cols, mt1->ldim)) == NULL) {
374 G_warning(_("Unable to allocate space for matrix product"));
375 return NULL;
376 }
377
378 /* Call the driver */
379
380 rows = (integer) mt1->rows;
381 interdim = (integer) mt1->cols;
382 cols = (integer) mt2->cols;
383
384 lda = (integer) mt1->ldim;
385 ldb = (integer) mt2->ldim;
386
387 f77_dgemm(&no_trans, &no_trans, &rows, &cols, &interdim, &unity,
388 mt1->vals, &lda, mt2->vals, &ldb, &zero, mt3->vals, &lda);
389
390 return mt3;
391}
392
393#else /* defined(HAVE_LIBBLAS) */
394#warning G_matrix_product() not compiled; requires BLAS library
395#endif /* defined(HAVE_LIBBLAS) */
396
397/*!
398 * \fn mat_struct *G_matrix_transpose (mat_struct *mt)
399 *
400 * \brief Transpose a matrix
401 *
402 * Transpose matrix <b>m1</b> by creating a new one and
403 * populating with transposed elements. The return matrix is
404 * automatically initialized.
405 *
406 * \param mt
407 * \return mat_struct
408 */
409
410mat_struct *G_matrix_transpose(mat_struct * mt)
411{
412 mat_struct *mt1;
413 int ldim, ldo;
414 doublereal *dbo, *dbt, *dbx, *dby;
415 int cnt, cnt2;
416
417 /* Word align the workspace blocks */
418 if (mt->cols % 2 == 0)
419 ldim = mt->cols;
420 else
421 ldim = mt->cols + 1;
422
423 mt1 = G_matrix_init(mt->cols, mt->rows, ldim);
424
425 /* Set initial values for reading arrays */
426 dbo = &mt->vals[0];
427 dbt = &mt1->vals[0];
428 ldo = mt->ldim;
429
430 for (cnt = 0; cnt < mt->cols; cnt++) {
431 dbx = dbo;
432 dby = dbt;
433
434 for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) {
435 *dby = *dbx;
436 dby += ldim;
437 dbx++;
438 }
439
440 *dby = *dbx;
441
442 if (cnt < mt->cols - 1) {
443 dbo += ldo;
444 dbt++;
445 }
446 }
447
448 return mt1;
449}
450
451
452#if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
453
454/*!
455 * \fn int G_matrix_LU_solve (const mat_struct *mt1, mat_struct **xmat0,
456 * const mat_struct *bmat, mat_type mtype)
457 *
458 * \brief Solve a general system A.X = B
459 *
460 * Solve a general system A.X = B, where A is a NxN matrix, X and B are
461 * NxC matrices, and we are to solve for C arrays in X given B. Uses LU
462 * decomposition.<br>
463 * Links to LAPACK function dgesv_() and similar to perform the core routine.
464 * (By default solves for a general non-symmetric matrix.)<br>
465 * mtype is a flag to indicate what kind of matrix (real/complex, Hermitian,
466 * symmetric, general etc.) is used (NONSYM, SYM, HERMITIAN).<br>
467 * <b>Warning:</b> NOT YET COMPLETE: only some solutions' options
468 * available. Now, only general real matrix is supported.
469 *
470 * \param mt1
471 * \param xmat0
472 * \param bmat
473 * \param mtype
474 * \return int
475 */
476
477/*** NOT YET COMPLETE: only some solutions' options available ***/
478
479int
480G_matrix_LU_solve(const mat_struct * mt1, mat_struct ** xmat0,
481 const mat_struct * bmat, mat_type mtype)
482{
483 mat_struct *wmat, *xmat, *mtx;
484
485 if (mt1->is_init == 0 || bmat->is_init == 0) {
486 G_warning(_("Input: one or both data matrices uninitialised"));
487 return -1;
488 }
489
490 if (mt1->rows != mt1->cols || mt1->rows < 1) {
491 G_warning(_("Principal matrix is not properly dimensioned"));
492 return -1;
493 }
494
495 if (bmat->cols < 1) {
496 G_warning(_("Input: you must have at least one array to solve"));
497 return -1;
498 }
499
500 /* Now create solution matrix by copying the original coefficient matrix */
501 if ((xmat = G_matrix_copy(bmat)) == NULL) {
502 G_warning(_("Could not allocate space for solution matrix"));
503 return -1;
504 }
505
506 /* Create working matrix for the coefficient array */
507 if ((mtx = G_matrix_copy(mt1)) == NULL) {
508 G_warning(_("Could not allocate space for working matrix"));
509 return -1;
510 }
511
512 /* Copy the contents of the data matrix, to preserve the
513 original information
514 */
515 if ((wmat = G_matrix_copy(bmat)) == NULL) {
516 G_warning(_("Could not allocate space for working matrix"));
517 return -1;
518 }
519
520 /* Now call appropriate LA driver to solve equations */
521 switch (mtype) {
522
523 case NONSYM:
524 {
525 integer *perm, res_info;
526 integer num_eqns, nrhs, lda, ldb;
527
528 perm = (integer *) G_malloc(wmat->rows * sizeof(integer));
529
530 /* Set fields to pass to fortran routine */
531 num_eqns = (integer) mt1->rows;
532 nrhs = (integer) wmat->cols;
533 lda = (integer) mt1->ldim;
534 ldb = (integer) wmat->ldim;
535
536 /* Call LA driver */
537 f77_dgesv(&num_eqns, &nrhs, mtx->vals, &lda, perm, wmat->vals,
538 &ldb, &res_info);
539
540 /* Copy the results from the modified data matrix, taking account
541 of pivot permutations ???
542 */
543
544 /*
545 for(indx1 = 0; indx1 < num_eqns; indx1++) {
546 iperm = perm[indx1];
547 ptin = &wmat->vals[0] + indx1;
548 ptout = &xmat->vals[0] + iperm;
549
550 for(indx2 = 0; indx2 < nrhs - 1; indx2++) {
551 *ptout = *ptin;
552 ptin += wmat->ldim;
553 ptout += xmat->ldim;
554 }
555
556 *ptout = *ptin;
557 }
558 */
559
560 memcpy(xmat->vals, wmat->vals,
561 wmat->cols * wmat->ldim * sizeof(doublereal));
562
563 /* Free temp arrays */
564 G_free(perm);
565 G_matrix_free(wmat);
566 G_matrix_free(mtx);
567
568 if (res_info > 0) {
569 G_warning(_("Matrix (or submatrix is singular). Solution undetermined"));
570 return 1;
571 }
572 else if (res_info < 0) {
573 G_warning(_("Problem in LA routine."));
574 return -1;
575 }
576 break;
577 }
578 default:
579 {
580 G_warning(_("Procedure not yet available for selected matrix type"));
581 return -1;
582 }
583 } /* end switch */
584
585 *xmat0 = xmat;
586
587 return 0;
588}
589
590#else /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
591#warning G_matrix_LU_solve() not compiled; requires BLAS and LAPACK libraries
592#endif /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
593
594#if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
595
596/*!
597 * \fn mat_struct *G_matrix_inverse (mat_struct *mt)
598 *
599 * \brief Returns the matrix inverse
600 *
601 * Calls G_matrix_LU_solve() to obtain matrix inverse using LU
602 * decomposition. Returns NULL on failure.
603 *
604 * \param mt
605 * \return mat_struct
606 */
607
608mat_struct *G_matrix_inverse(mat_struct * mt)
609{
610 mat_struct *mt0, *res;
611 int i, j, k; /* loop */
612
613 if (mt->rows != mt->cols) {
614 G_warning(_("Matrix is not square. Cannot determine inverse"));
615 return NULL;
616 }
617
618 if ((mt0 = G_matrix_init(mt->rows, mt->rows, mt->ldim)) == NULL) {
619 G_warning(_("Unable to allocate space for matrix"));
620 return NULL;
621 }
622
623 /* Set `B' matrix to unit matrix */
624 for (i = 0; i < mt0->rows - 1; i++) {
625 mt0->vals[i + i * mt0->ldim] = 1.0;
626
627 for (j = i + 1; j < mt0->cols; j++) {
628 mt0->vals[i + j * mt0->ldim] = mt0->vals[j + i * mt0->ldim] = 0.0;
629 }
630 }
631
632 mt0->vals[mt0->rows - 1 + (mt0->rows - 1) * mt0->ldim] = 1.0;
633
634 /* Solve system */
635 if ((k = G_matrix_LU_solve(mt, &res, mt0, NONSYM)) == 1) {
636 G_warning(_("Matrix is singular"));
637 G_matrix_free(mt0);
638 return NULL;
639 }
640 else if (k < 0) {
641 G_warning(_("Problem in LA procedure."));
642 G_matrix_free(mt0);
643 return NULL;
644 }
645 else {
646 G_matrix_free(mt0);
647 return res;
648 }
649}
650
651#else /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
652#warning G_matrix_inverse() not compiled; requires BLAS and LAPACK libraries
653#endif /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */
654
655
656/*!
657 * \fn void G_matrix_free (mat_struct *mt)
658 *
659 * \brief Free up allocated matrix
660 *
661 * Free up allocated matrix.
662 *
663 * \param mt
664 * \return void
665 */
666
667void G_matrix_free(mat_struct * mt)
668{
669 if (mt->is_init)
670 G_free(mt->vals);
671
672 G_free(mt);
673}
674
675
676/*!
677 * \fn void G_matrix_print (mat_struct *mt)
678 *
679 * \brief Print out a matrix
680 *
681 * Print out a representation of the matrix to standard output.
682 *
683 * \param mt
684 * \return void
685 */
686
687void G_matrix_print(mat_struct * mt)
688{
689 int i, j;
690 char buf[64], numbuf[64];
691
692 for (i = 0; i < mt->rows; i++) {
693 strcpy(buf, "");
694
695 for (j = 0; j < mt->cols; j++) {
696
697 sprintf(numbuf, "%14.6f", G_matrix_get_element(mt, i, j));
698 strcat(buf, numbuf);
699 if (j < mt->cols - 1)
700 strcat(buf, ", ");
701 }
702
703 G_message("%s", buf);
704 }
705
706 fprintf(stderr, "\n");
707}
708
709
710/*!
711 * \fn int G_matrix_set_element (mat_struct *mt, int rowval, int colval, double val)
712 *
713 * \brief Set the value of the (i, j)th element
714 *
715 * Set the value of the (i, j)th element to a double value. Index values
716 * are C-like ie. zero-based. The row number is given first as is
717 * conventional. Returns -1 if the accessed cell is outside the bounds.
718 *
719 * \param mt
720 * \param rowval
721 * \param colval
722 * \param val
723 * \return int
724 */
725
726int G_matrix_set_element(mat_struct * mt, int rowval, int colval, double val)
727{
728 if (!mt->is_init) {
729 G_warning(_("Element array has not been allocated"));
730 return -1;
731 }
732
733 if (rowval >= mt->rows || colval >= mt->cols || rowval < 0 || colval < 0) {
734 G_warning(_("Specified element is outside array bounds"));
735 return -1;
736 }
737
738 mt->vals[rowval + colval * mt->ldim] = (doublereal) val;
739
740 return 0;
741}
742
743
744/*!
745 * \fn double G_matrix_get_element (mat_struct *mt, int rowval, int colval)
746 *
747 * \brief Retrieve value of the (i,j)th element
748 *
749 * Retrieve the value of the (i, j)th element to a double value. Index
750 * values are C-like ie. zero-based.
751 * <b>Note:</b> Does currently not set an error flag for bounds checking.
752 *
753 * \param mt
754 * \param rowval
755 * \param colval
756 * \return double
757 */
758
759double G_matrix_get_element(mat_struct * mt, int rowval, int colval)
760{
761 double val;
762
763 /* Should do some checks, but this would require an error control
764 system: later? */
765
766 return (val = (double)mt->vals[rowval + colval * mt->ldim]);
767}
768
769
770/*!
771 * \fn vec_struct *G_matvect_get_column (mat_struct *mt, int col)
772 *
773 * \brief Retrieve a column of the matrix to a vector structure
774 *
775 * Retrieve a column of matrix <b>mt</b> to a returning vector structure
776 *
777 * \param mt
778 * \param col
779 * \return vec_struct
780 */
781
782vec_struct *G_matvect_get_column(mat_struct * mt, int col)
783{
784 int i; /* loop */
785 vec_struct *vc1;
786
787 if (col < 0 || col >= mt->cols) {
788 G_warning(_("Specified matrix column index is outside range"));
789 return NULL;
790 }
791
792 if (!mt->is_init) {
793 G_warning(_("Matrix is not initialised"));
794 return NULL;
795 }
796
797 if ((vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) == NULL) {
798 G_warning(_("Could not allocate space for vector structure"));
799 return NULL;
800 }
801
802 for (i = 0; i < mt->rows; i++)
803 G_matrix_set_element((mat_struct *) vc1, i, 0,
804 G_matrix_get_element(mt, i, col));
805
806 return vc1;
807}
808
809
810/*!
811 * \fn vec_struct *G_matvect_get_row (mat_struct *mt, int row)
812 *
813 * \brief Retrieve a row of the matrix to a vector structure
814 *
815 * Retrieves a row from matrix <b>mt</b> and returns it in a vector
816 * structure.
817 *
818 * \param mt
819 * \param row
820 * \return vec_struct
821 */
822
823vec_struct *G_matvect_get_row(mat_struct * mt, int row)
824{
825 int i; /* loop */
826 vec_struct *vc1;
827
828 if (row < 0 || row >= mt->cols) {
829 G_warning(_("Specified matrix row index is outside range"));
830 return NULL;
831 }
832
833 if (!mt->is_init) {
834 G_warning(_("Matrix is not initialised"));
835 return NULL;
836 }
837
838 if ((vc1 = G_vector_init(mt->cols, mt->ldim, RVEC)) == NULL) {
839 G_warning(_("Could not allocate space for vector structure"));
840 return NULL;
841 }
842
843 for (i = 0; i < mt->cols; i++)
844 G_matrix_set_element((mat_struct *) vc1, 0, i,
845 G_matrix_get_element(mt, row, i));
846
847 return vc1;
848}
849
850
851/*!
852 * \fn int G_matvect_extract_vector (mat_struct *mt, vtype vt, int indx)
853 *
854 * \brief Convert matrix to vector
855 *
856 * Convert the matrix <b>mt</b> to a vector structure. The vtype,
857 * <b>vt</b>, is RVEC or CVEC which specifies a row vector or column
858 * vector. The index, <b>indx</b>, indicates the row/column number (zero based).
859 *
860 * \param mt
861 * \param vt
862 * \param indx
863 * \return int
864 */
865
866int G_matvect_extract_vector(mat_struct * mt, vtype vt, int indx)
867{
868 if (vt == RVEC && indx >= mt->rows) {
869 G_warning(_("Specified row index is outside range"));
870 return -1;
871 }
872
873 else if (vt == CVEC && indx >= mt->cols) {
874 G_warning(_("Specified column index is outside range"));
875 return -1;
876 }
877
878 switch (vt) {
879
880 case RVEC:
881 {
882 mt->type = ROWVEC_;
883 mt->v_indx = indx;
884 }
885
886 case CVEC:
887 {
888 mt->type = COLVEC_;
889 mt->v_indx = indx;
890 }
891
892 default:
893 {
894 G_warning(_("Unknown vector type."));
895 return -1;
896 }
897
898 }
899
900 return 0;
901
902}
903
904
905/*!
906 * \fn int G_matvect_retrieve_matrix (vec_struct *vc)
907 *
908 * \brief Revert a vector to matrix
909 *
910 * Revert vector <b>vc</b> to a matrix.
911 *
912 * \param vc
913 * \return int
914 */
915
916int G_matvect_retrieve_matrix(vec_struct * vc)
917{
918 /* We have to take the integrity of the vector structure
919 largely on trust
920 */
921
922 vc->type = MATRIX_;
923 vc->v_indx = -1;
924
925 return 0;
926}
927
928
929/*!
930 * \fn vec_struct *G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
931 *
932 * \brief Calculates the matrix-vector product
933 *
934 * Calculates the product of a matrix and a vector
935 *
936 * \param A
937 * \param b
938 * \return vec_struct
939 */
940
941vec_struct *G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
942{
943 unsigned int i, m, n, j;
944 register doublereal sum;
945
946/* G_message("A=%d,%d,%d", A->cols, A->rows, A->ldim); */
947/* G_message("B=%d,%d,%d", b->cols, b->rows, b->ldim); */
948 if (A->cols != b->cols) {
949 G_warning (_("Input matrix and vector have differing dimensions1"));
950
951 return NULL;
952
953 }
954 if (!out) {
955 G_warning (_("Output vector is uninitialized"));
956 return NULL;
957 }
958/* if (out->ldim != A->rows) {*/
959/* G_warning (_("Output vector has incorrect dimension"));*/
960/* exit(1);*/
961/* return NULL;*/
962/* }*/
963
964 m = A->rows;
965 n = A->cols;
966
967 for (i = 0; i < m; i++) {
968 sum = 0.0;
969 int width = A->rows;
970 for (j = 0; j < n; j++) {
971
972 sum +=G_matrix_get_element(A, i, j) * G_matrix_get_element(b, 0, j);
973 /*sum += A->vals[i + j * width] * b->vals[j];*/
974 out->vals[i] = sum;
975 }
976 }
977 return (out);
978}
979
980
981/*!
982 * \fn vec_struct *G_vector_init (int cells, int ldim, vtype vt)
983 *
984 * \brief Initialize a vector structure
985 *
986 * Returns an initialized vector structure with <b>cell</b> cells,
987 * of dimension <b>ldim</b>, and of type <b>vt</b>.
988 *
989 * \param cells
990 * \param ldim
991 * \param vt
992 * \return vec_struct
993 */
994
995vec_struct *G_vector_init(int cells, int ldim, vtype vt)
996{
997 vec_struct *tmp_arry;
998
999 if ((cells < 1) || (vt == RVEC && ldim < 1)
1000 || (vt == CVEC && ldim < cells) || ldim < 0) {
1001 G_warning("Vector dimensions out of range.");
1002 return NULL;
1003 }
1004
1005 tmp_arry = (vec_struct *) G_malloc(sizeof(vec_struct));
1006
1007 if (vt == RVEC) {
1008 tmp_arry->rows = 1;
1009 tmp_arry->cols = cells;
1010 tmp_arry->ldim = ldim;
1011 tmp_arry->type = ROWVEC_;
1012 }
1013
1014 else if (vt == CVEC) {
1015 tmp_arry->rows = cells;
1016 tmp_arry->cols = 1;
1017 tmp_arry->ldim = ldim;
1018 tmp_arry->type = COLVEC_;
1019 }
1020
1021 tmp_arry->v_indx = 0;
1022
1023 tmp_arry->vals = (doublereal *) G_calloc(ldim * tmp_arry->cols,
1024 sizeof(doublereal));
1025 tmp_arry->is_init = 1;
1026
1027 return tmp_arry;
1028}
1029
1030
1031/*!
1032 * \fn void G_vector_free (vec_struct *v)
1033 *
1034 * \brief Free an allocated vector structure
1035 *
1036 * Free an allocated vector structure.
1037 *
1038 * \param v
1039 * \return void
1040 */
1041
1042void G_vector_free(vec_struct * v)
1043{
1044 if (v->is_init)
1045 G_free(v->vals);
1046
1047 G_free(v);
1048}
1049
1050
1051/*!
1052 * \fn vec_struct *G_vector_sub (vec_struct *v1, vec_struct *v2, vec_struct *out)
1053 *
1054 * \brief Subtract two vectors
1055 *
1056 * Subtracts two vectors, <b>v1</b> and <b>v2</b>, and returns and
1057 * populates vector <b>out</b>.
1058 *
1059 * \param v1
1060 * \param v2
1061 * \param out
1062 * \return vec_struct
1063 */
1064
1065vec_struct *G_vector_sub(vec_struct * v1, vec_struct * v2, vec_struct * out)
1066{
1067 int idx1, idx2, idx0;
1068 int i;
1069
1070 if (!out->is_init) {
1071 G_warning(_("Output vector is uninitialized"));
1072 return NULL;
1073 }
1074
1075 if (v1->type != v2->type) {
1076 G_warning(_("Vectors are not of the same type"));
1077 return NULL;
1078 }
1079
1080 if (v1->type != out->type) {
1081 G_warning(_("Output vector is of incorrect type"));
1082 return NULL;
1083 }
1084
1085 if (v1->type == MATRIX_) {
1086 G_warning(_("Matrices not allowed"));
1087 return NULL;
1088 }
1089
1090 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1091 (v1->type == COLVEC_ && v1->rows != v2->rows)) {
1092 G_warning(_("Vectors have differing dimensions"));
1093 return NULL;
1094 }
1095
1096 if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1097 (v1->type == COLVEC_ && v1->rows != out->rows)) {
1098 G_warning(_("Output vector has incorrect dimension"));
1099 return NULL;
1100 }
1101
1102 idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
1103 idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
1104 idx0 = (out->v_indx > 0) ? out->v_indx : 0;
1105
1106 if (v1->type == ROWVEC_) {
1107 for (i = 0; i < v1->cols; i++)
1108 G_matrix_set_element(out, idx0, i,
1109 G_matrix_get_element(v1, idx1, i) -
1110 G_matrix_get_element(v2, idx2, i));
1111 }
1112 else {
1113 for (i = 0; i < v1->rows; i++)
1114 G_matrix_set_element(out, i, idx0,
1115 G_matrix_get_element(v1, i, idx1) -
1116 G_matrix_get_element(v2, i, idx2));
1117 }
1118
1119 return out;
1120}
1121
1122/*!
1123 * \fn int G_vector_set (vec_struct *A, int cells, int ldim, vtype vt, int vindx)
1124 *
1125 * \brief Set parameters for vector structure
1126 *
1127 * Set parameters for a vector structure that is
1128 * allocated but not yet initialised fully. The vtype is RVEC or
1129 * CVEC which specifies a row vector or column vector.
1130 *
1131 * \param A
1132 * \param cells
1133 * \param ldim
1134 * \param vt
1135 * \param vindx
1136 * \return int
1137 */
1138
1139int G_vector_set(vec_struct * A, int cells, int ldim, vtype vt, int vindx)
1140{
1141 if ((cells < 1) || (vt == RVEC && ldim < 1)
1142 || (vt == CVEC && ldim < cells) || ldim < 0) {
1143 G_warning(_("Vector dimensions out of range"));
1144 return -1;
1145 }
1146
1147 if ((vt == RVEC && vindx >= A->cols) || (vt == CVEC && vindx >= A->rows)) {
1148 G_warning(_("Row/column out of range"));
1149 return -1;
1150 }
1151
1152 if (vt == RVEC) {
1153 A->rows = 1;
1154 A->cols = cells;
1155 A->ldim = ldim;
1156 A->type = ROWVEC_;
1157 }
1158 else {
1159 A->rows = cells;
1160 A->cols = 1;
1161 A->ldim = ldim;
1162 A->type = COLVEC_;
1163 }
1164
1165 if (vindx < 0)
1166 A->v_indx = 0;
1167 else
1168 A->v_indx = vindx;
1169
1170 A->vals = (doublereal *) G_calloc(ldim * A->cols, sizeof(doublereal));
1171 A->is_init = 1;
1172
1173 return 0;
1174
1175}
1176
1177
1178#if defined(HAVE_LIBBLAS)
1179
1180/*!
1181 * \fn double G_vector_norm_euclid (vec_struct *vc)
1182 *
1183 * \brief Calculates euclidean norm
1184 *
1185 * Calculates the euclidean norm of a row or column vector, using BLAS
1186 * routine dnrm2_().
1187 *
1188 * \param vc
1189 * \return double
1190 */
1191
1192double G_vector_norm_euclid(vec_struct * vc)
1193{
1194 integer incr, Nval;
1195 doublereal *startpt;
1196
1197 if (!vc->is_init)
1198 G_fatal_error(_("Matrix is not initialised"));
1199
1200 if (vc->type == ROWVEC_) {
1201 Nval = (integer) vc->cols;
1202 incr = (integer) vc->ldim;
1203 if (vc->v_indx < 0)
1204 startpt = vc->vals;
1205 else
1206 startpt = vc->vals + vc->v_indx;
1207 }
1208 else {
1209 Nval = (integer) vc->rows;
1210 incr = 1;
1211 if (vc->v_indx < 0)
1212 startpt = vc->vals;
1213 else
1214 startpt = vc->vals + vc->v_indx * vc->ldim;
1215 }
1216
1217 /* Call the BLAS routine dnrm2_() */
1218 return (double)f77_dnrm2(&Nval, startpt, &incr);
1219}
1220
1221#else /* defined(HAVE_LIBBLAS) */
1222#warning G_vector_norm_euclid() not compiled; requires BLAS library
1223#endif /* defined(HAVE_LIBBLAS) */
1224
1225
1226/*!
1227 * \fn double G_vector_norm_maxval (vec_struct *vc, int vflag)
1228 *
1229 * \brief Calculates maximum value
1230 *
1231 * Calculates the maximum value of a row or column vector.
1232 * The vflag setting defines which value to be calculated:
1233 * vflag:
1234 * 1 Indicates maximum value<br>
1235 * -1 Indicates minimum value<br>
1236 * 0 Indicates absolute value [???]
1237 *
1238 * \param vc
1239 * \param vflag
1240 * \return double
1241 */
1242
1243double G_vector_norm_maxval(vec_struct * vc, int vflag)
1244{
1245 doublereal xval, *startpt, *curpt;
1246 double cellval;
1247 int ncells, incr;
1248
1249 if (!vc->is_init)
1250 G_fatal_error(_("Matrix is not initialised"));
1251
1252 if (vc->type == ROWVEC_) {
1253 ncells = (integer) vc->cols;
1254 incr = (integer) vc->ldim;
1255 if (vc->v_indx < 0)
1256 startpt = vc->vals;
1257 else
1258 startpt = vc->vals + vc->v_indx;
1259 }
1260 else {
1261 ncells = (integer) vc->rows;
1262 incr = 1;
1263 if (vc->v_indx < 0)
1264 startpt = vc->vals;
1265 else
1266 startpt = vc->vals + vc->v_indx * vc->ldim;
1267 }
1268
1269 xval = *startpt;
1270 curpt = startpt;
1271
1272 while (ncells > 0) {
1273 if (curpt != startpt) {
1274 switch (vflag) {
1275
1276 case MAX_POS:
1277 {
1278 if (*curpt > xval)
1279 xval = *curpt;
1280 break;
1281 }
1282
1283 case MAX_NEG:
1284 {
1285 if (*curpt < xval)
1286 xval = *curpt;
1287 break;
1288 }
1289
1290 case MAX_ABS:
1291 {
1292 cellval = (double)(*curpt);
1293 if (hypot(cellval, cellval) > (double)xval)
1294 xval = *curpt;
1295 }
1296 } /* switch */
1297 } /* if(curpt != startpt) */
1298
1299 curpt += incr;
1300 ncells--;
1301 }
1302
1303 return (double)xval;
1304}
1305
1306
1307/*!
1308 * \fn double G_vector_norm1 (vec_struct *vc)
1309 *
1310 * \brief Calculates the 1-norm of a vector
1311 *
1312 * Calculates the 1-norm of a vector
1313 *
1314 * \param vc
1315 * \return double
1316 */
1317
1318double G_vector_norm1(vec_struct * vc)
1319{
1320 double result = 0.0;
1321 int idx;
1322 int i;
1323
1324 if (!vc->is_init) {
1325 G_warning(_("Matrix is not initialised"));
1326 return 0.0 / 0.0; /* NaN */
1327 }
1328
1329 idx = (vc->v_indx > 0) ? vc->v_indx : 0;
1330
1331 if (vc->type == ROWVEC_) {
1332 for (i = 0; i < vc->cols; i++)
1333 result += fabs(G_matrix_get_element(vc, idx, i));
1334 }
1335 else {
1336 for (i = 0; i < vc->rows; i++)
1337 result += fabs(G_matrix_get_element(vc, i, idx));
1338 }
1339
1340 return result;
1341}
1342
1343/*!
1344 * \fn vec_struct *G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out)
1345 *
1346 * \brief Calculates the vector product
1347 *
1348 * Calculates the vector product of two vectors
1349 *
1350 * \param v1
1351 * \param v2
1352 * \return vec_struct
1353 */
1354
1355vec_struct *G_vector_product (vec_struct *v1, vec_struct *v2, vec_struct *out)
1356{
1357 int idx1, idx2, idx0;
1358 int i;
1359
1360 if (!out->is_init) {
1361 G_warning (_("Output vector is uninitialized"));
1362 return NULL;
1363 }
1364
1365 if (v1->type != v2->type) {
1366 G_warning (_("Vectors are not of the same type"));
1367 return NULL;
1368 }
1369
1370 if (v1->type != out->type) {
1371 G_warning (_("Output vector is not the same type as others"));
1372 return NULL;
1373 }
1374
1375 if (v1->type == MATRIX_) {
1376 G_warning (_("Matrices not allowed"));
1377 return NULL;
1378 }
1379
1380 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) ||
1381 (v1->type == COLVEC_ && v1->rows != v2->rows))
1382 {
1383 G_warning (_("Vectors have differing dimensions"));
1384 return NULL;
1385 }
1386
1387 if ((v1->type == ROWVEC_ && v1->cols != out->cols) ||
1388 (v1->type == COLVEC_ && v1->rows != out->rows))
1389 {
1390 G_warning (_("Output vector has incorrect dimension"));
1391 return NULL;
1392 }
1393
1394#if defined(HAVE_LAPACK) && defined(HAVE_LIBBLAS)
1395 f77_dhad (v1->cols, 1.0, v1->vals, 1, v2->vals, 1, 0.0, out->vals, 1.0);
1396#else
1397 idx1 = (v1->v_indx > 0) ? v1->v_indx : 0;
1398 idx2 = (v2->v_indx > 0) ? v2->v_indx : 0;
1399 idx0 = (out->v_indx > 0) ? out->v_indx : 0;
1400
1401 if (v1->type == ROWVEC_) {
1402 for (i = 0; i < v1->cols; i++)
1403 G_matrix_set_element(out, idx0, i,
1404 G_matrix_get_element(v1, idx1, i) *
1405 G_matrix_get_element(v2, idx2, i));
1406 } else {
1407 for (i = 0; i < v1->rows; i++)
1408 G_matrix_set_element(out, i, idx0,
1409 G_matrix_get_element(v1, i, idx1) *
1410 G_matrix_get_element(v2, i, idx2));
1411 }
1412#endif
1413
1414 return out;
1415}
1416
1417
1418/*!
1419 * \fn vec_struct *G_vector_copy (const vec_struct *vc1, int comp_flag)
1420 *
1421 * \brief Returns a vector copied from <b>vc1</b>. Underlying structure
1422 * is preserved unless DO_COMPACT flag.
1423 *
1424 * \param vc1
1425 * \param comp_flag
1426 * \return vec_struct
1427 */
1428
1429vec_struct *G_vector_copy(const vec_struct * vc1, int comp_flag)
1430{
1431 vec_struct *tmp_arry;
1432 int incr1, incr2;
1433 doublereal *startpt1, *startpt2, *curpt1, *curpt2;
1434 int cnt;
1435
1436 if (!vc1->is_init) {
1437 G_warning(_("Vector structure is not initialised"));
1438 return NULL;
1439 }
1440
1441 tmp_arry = (vec_struct *) G_malloc(sizeof(vec_struct));
1442
1443 if (comp_flag == DO_COMPACT) {
1444 if (vc1->type == ROWVEC_) {
1445 tmp_arry->rows = 1;
1446 tmp_arry->cols = vc1->cols;
1447 tmp_arry->ldim = 1;
1448 tmp_arry->type = ROWVEC_;
1449 tmp_arry->v_indx = 0;
1450 }
1451 else if (vc1->type == COLVEC_) {
1452 tmp_arry->rows = vc1->rows;
1453 tmp_arry->cols = 1;
1454 tmp_arry->ldim = vc1->ldim;
1455 tmp_arry->type = COLVEC_;
1456 tmp_arry->v_indx = 0;
1457 }
1458 else {
1459 G_warning("Type is not vector.");
1460 return NULL;
1461 }
1462 }
1463 else if (comp_flag == NO_COMPACT) {
1464 tmp_arry->v_indx = vc1->v_indx;
1465 tmp_arry->rows = vc1->rows;
1466 tmp_arry->cols = vc1->cols;
1467 tmp_arry->ldim = vc1->ldim;
1468 tmp_arry->type = vc1->type;
1469 }
1470 else {
1471 G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n");
1472 return NULL;
1473 }
1474
1475 tmp_arry->vals = (doublereal *) G_calloc(tmp_arry->ldim * tmp_arry->cols,
1476 sizeof(doublereal));
1477 if (comp_flag == DO_COMPACT) {
1478 if (tmp_arry->type == ROWVEC_) {
1479 startpt1 = tmp_arry->vals;
1480 startpt2 = vc1->vals + vc1->v_indx;
1481 curpt1 = startpt1;
1482 curpt2 = startpt2;
1483 incr1 = 1;
1484 incr2 = vc1->ldim;
1485 cnt = vc1->cols;
1486 }
1487 else if (tmp_arry->type == COLVEC_) {
1488 startpt1 = tmp_arry->vals;
1489 startpt2 = vc1->vals + vc1->v_indx * vc1->ldim;
1490 curpt1 = startpt1;
1491 curpt2 = startpt2;
1492 incr1 = 1;
1493 incr2 = 1;
1494 cnt = vc1->rows;
1495 }
1496 else {
1497 G_warning("Structure type is not vector.");
1498 return NULL;
1499 }
1500 }
1501 else if (comp_flag == NO_COMPACT) {
1502 startpt1 = tmp_arry->vals;
1503 startpt2 = vc1->vals;
1504 curpt1 = startpt1;
1505 curpt2 = startpt2;
1506 incr1 = 1;
1507 incr2 = 1;
1508 cnt = vc1->ldim * vc1->cols;
1509 }
1510 else {
1511 G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n");
1512 return NULL;
1513 }
1514
1515 while (cnt > 0) {
1516 memcpy(curpt1, curpt2, sizeof(doublereal));
1517 curpt1 += incr1;
1518 curpt2 += incr2;
1519 cnt--;
1520 }
1521
1522 tmp_arry->is_init = 1;
1523
1524 return tmp_arry;
1525}
1526
1527
1528/*!
1529 * \fn int G_matrix_read (FILE *fp, mat_struct *out)
1530 *
1531 * \brief Read a matrix from a file stream
1532 *
1533 * Populates matrix structure <b>out</b> with matrix read from file
1534 * stream <b>fp</b>. Matrix <b>out</b> is automatically initialized.
1535 * Returns -1 on error and 0 on success.
1536 *
1537 * \param fp
1538 * \param out
1539 * \return int
1540 */
1541
1542int G_matrix_read(FILE * fp, mat_struct * out)
1543{
1544 char buff[100];
1545 int rows, cols;
1546 int i, j, row;
1547 double val;
1548
1549 /* skip comments */
1550 for (;;) {
1551 if (!G_getl(buff, sizeof(buff), fp))
1552 return -1;
1553 if (buff[0] != '#')
1554 break;
1555 }
1556
1557 if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) {
1558 G_warning(_("Input format error"));
1559 return -1;
1560 }
1561
1562 G_matrix_set(out, rows, cols, rows);
1563
1564 for (i = 0; i < rows; i++) {
1565 if (fscanf(fp, "row%d:", &row) != 1 || row != i) {
1566 G_warning(_("Input format error"));
1567 return -1;
1568 }
1569 for (j = 0; j < cols; j++) {
1570 if (fscanf(fp, "%lf:", &val) != 1) {
1571 G_warning(_("Input format error"));
1572 return -1;
1573 }
1574
1575 G_matrix_set_element(out, i, j, val);
1576 }
1577 }
1578
1579 return 0;
1580}
1581
1582
1583/*!
1584 * \fn mat_struct *G_matrix_resize(mat_struct *in, int rows, int cols)
1585 *
1586 * \brief Resizes a matrix
1587 *
1588 * Resizes a matrix
1589 *
1590 * \param A
1591 * \param rows
1592 * \param cols
1593 * \return mat_struct
1594 */
1595
1596mat_struct *G_matrix_resize(mat_struct *in, int rows, int cols)
1597{
1598 mat_struct *matrix;
1599 matrix = G_matrix_init(rows, cols, rows);
1600 int i, j, p, index = 0;
1601 for (i = 0; i < rows; i++)
1602 for (j = 0; j < cols; j++)
1603 G_matrix_set_element(matrix, i, j, G_matrix_get_element(in, i, j));
1604/* matrix->vals[index++] = in->vals[i + j * cols];*/
1605
1606 int old_size = in->rows * in->cols;
1607 int new_size = rows * cols;
1608
1609 if (new_size > old_size)
1610 for (p = old_size; p < new_size; p++)
1611 G_matrix_set_element(matrix, i, j, 0.0);
1612
1613 return (matrix);
1614}
1615
1616
1617/*!
1618 * \fn int G_matrix_read_stdin (mat_struct *out)
1619 *
1620 * \brief Read a matrix from standard input
1621 *
1622 * Populates matrix <b>out</b> with matrix read from stdin. Matrix
1623 * <b>out</b> is automatically initialized. Returns -1 on failure or 0
1624 * on success.
1625 *
1626 * \param out
1627 * \return int
1628 */
1629
1630int G_matrix_stdin(mat_struct * out)
1631{
1632 return G_matrix_read(stdin, out);
1633}
1634
1635
1636/*!
1637 * \fn int G_matrix_eigen_sort (vec_struct *d, mat_struct *m)
1638 *
1639 * \brief Sort eigenvectors according to eigenvalues
1640 *
1641 * Sort eigenvectors according to eigenvalues. Returns 0.
1642 *
1643 * \param d
1644 * \param m
1645 * \return int
1646 */
1647
1648int G_matrix_eigen_sort(vec_struct * d, mat_struct * m)
1649{
1650 mat_struct tmp;
1651 int i, j;
1652 int idx;
1653
1654 G_matrix_set(&tmp, m->rows + 1, m->cols, m->ldim + 1);
1655
1656 idx = (d->v_indx > 0) ? d->v_indx : 0;
1657
1658 /* concatenate (vertically) m and d into tmp */
1659 for (i = 0; i < m->cols; i++) {
1660 for (j = 0; j < m->rows; j++)
1661 G_matrix_set_element(&tmp, j + 1, i,
1662 G_matrix_get_element(m, j, i));
1663 if (d->type == ROWVEC_)
1664 G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, idx, i));
1665 else
1666 G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, i, idx));
1667 }
1668
1669 /* sort the combined matrix */
1670 qsort(tmp.vals, tmp.cols, tmp.ldim * sizeof(doublereal), egcmp);
1671
1672 /* split tmp into m and d */
1673 for (i = 0; i < m->cols; i++) {
1674 for (j = 0; j < m->rows; j++)
1675 G_matrix_set_element(m, j, i,
1676 G_matrix_get_element(&tmp, j + 1, i));
1677 if (d->type == ROWVEC_)
1678 G_matrix_set_element(d, idx, i, G_matrix_get_element(&tmp, 0, i));
1679 else
1680 G_matrix_set_element(d, i, idx, G_matrix_get_element(&tmp, 0, i));
1681 }
1682
1683 G_free(tmp.vals);
1684
1685 return 0;
1686}
1687
1688
1689static int egcmp(const void *pa, const void *pb)
1690{
1691 double a = *(doublereal * const)pa;
1692 double b = *(doublereal * const)pb;
1693
1694 if (a > b)
1695 return 1;
1696 if (a < b)
1697 return -1;
1698
1699 return 0;
1700}
1701
1702
1703#endif /* HAVE_BLAS && HAVE_LAPACK && HAVE_G2C */
void G_free(void *buf)
Free allocated memory.
Definition: alloc.c:149
#define NULL
Definition: ccmath.h:32
if(!DBFLoadRecord(psDBF, hEntity)) return NULL
double b
int G_getl(char *buf, int n, FILE *fd)
Gets a line of text from a file.
Definition: getl.c:31
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
void G_message(const char *msg,...)
Print a message to stderr.
Definition: gis/error.c:90
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
vec_struct * G_matvect_get_column(mat_struct *mt, int col)
Retrieve a column of the matrix to a vector structure.
Definition: la.c:782
mat_struct * G_matrix_copy(const mat_struct *A)
Copy a matrix.
Definition: la.c:145
vec_struct * G_vector_sub(vec_struct *v1, vec_struct *v2, vec_struct *out)
Subtract two vectors.
Definition: la.c:1065
double G_vector_norm_maxval(vec_struct *vc, int vflag)
Calculates maximum value.
Definition: la.c:1243
mat_struct * G__matrix_add(mat_struct *mt1, mat_struct *mt2, const double c1, const double c2)
General add/subtract/scalar multiply routine.
Definition: la.c:278
vec_struct * G_vector_product(vec_struct *v1, vec_struct *v2, vec_struct *out)
Calculates the vector product.
Definition: la.c:1355
mat_struct * G_matrix_subtract(mat_struct *mt1, mat_struct *mt2)
Subtract two matricies.
Definition: la.c:197
int G_matvect_retrieve_matrix(vec_struct *vc)
Revert a vector to matrix.
Definition: la.c:916
double G_vector_norm_euclid(vec_struct *vc)
Calculates euclidean norm.
Definition: la.c:1192
vec_struct * G_vector_init(int cells, int ldim, vtype vt)
Initialize a vector structure.
Definition: la.c:995
mat_struct * G_matrix_inverse(mat_struct *mt)
Returns the matrix inverse.
Definition: la.c:608
double G_matrix_get_element(mat_struct *mt, int rowval, int colval)
Retrieve value of the (i,j)th element.
Definition: la.c:759
mat_struct * G_matrix_scale(mat_struct *mt1, const double c)
Scale a matrix by a scalar value.
Definition: la.c:257
int G_matrix_stdin(mat_struct *out)
Definition: la.c:1630
mat_struct * G_matrix_resize(mat_struct *in, int rows, int cols)
Resizes a matrix.
Definition: la.c:1596
void G_matrix_print(mat_struct *mt)
Print out a matrix.
Definition: la.c:687
mat_struct * G_matrix_scalar_mul(double scalar, mat_struct *matrix, mat_struct *out)
Calculates the scalar-matrix multiplication.
Definition: la.c:214
vec_struct * G_matvect_product(mat_struct *A, vec_struct *b, vec_struct *out)
Calculates the matrix-vector product.
Definition: la.c:941
mat_struct * G_matrix_add(mat_struct *mt1, mat_struct *mt2)
Adds two matricies.
Definition: la.c:178
int G_matvect_extract_vector(mat_struct *mt, vtype vt, int indx)
Convert matrix to vector.
Definition: la.c:866
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.
Definition: la.c:480
mat_struct * G_matrix_init(int rows, int cols, int ldim)
Initialize a matrix structure.
Definition: la.c:56
double G_vector_norm1(vec_struct *vc)
Calculates the 1-norm of a vector.
Definition: la.c:1318
int G_matrix_read(FILE *fp, mat_struct *out)
Read a matrix from a file stream.
Definition: la.c:1542
vec_struct * G_matvect_get_row(mat_struct *mt, int row)
Retrieve a row of the matrix to a vector structure.
Definition: la.c:823
int G_matrix_eigen_sort(vec_struct *d, mat_struct *m)
Sort eigenvectors according to eigenvalues.
Definition: la.c:1648
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.
Definition: la.c:1429
void G_vector_free(vec_struct *v)
Free an allocated vector structure.
Definition: la.c:1042
void G_matrix_free(mat_struct *mt)
Free up allocated matrix.
Definition: la.c:667
mat_struct * G_matrix_product(mat_struct *mt1, mat_struct *mt2)
Returns product of two matricies.
Definition: la.c:356
int G_vector_set(vec_struct *A, int cells, int ldim, vtype vt, int vindx)
Set parameters for vector structure.
Definition: la.c:1139
mat_struct * G_matrix_transpose(mat_struct *mt)
Transpose a matrix.
Definition: la.c:410
int G_matrix_set_element(mat_struct *mt, int rowval, int colval, double val)
Set the value of the (i, j)th element.
Definition: la.c:726
int G_matrix_set(mat_struct *A, int rows, int cols, int ldim)
Set parameters for an initialized matrix.
Definition: la.c:114
int G_matrix_zero(mat_struct *A)
Clears (or resets) the matrix values to 0.
Definition: la.c:88