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