GRASS Programmer's Manual
6.4.2(2012)
|
00001 00002 /****************************************************************************** 00003 * la.c 00004 * wrapper modules for linear algebra problems 00005 * linking to BLAS / LAPACK (and others?) 00006 00007 * @Copyright David D.Gray <ddgray@armadce.demon.co.uk> 00008 * 26th. Sep. 2000 00009 * Last updated: 00010 * 2006-11-23 00011 00012 * This file is part of GRASS GIS. It is free software. You can 00013 * redistribute it and/or modify it under the terms of 00014 * the GNU General Public License as published by the Free Software 00015 * Foundation; either version 2 of the License, or (at your option) 00016 * any later version. 00017 00018 * This program is distributed in the hope that it will be useful, 00019 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00020 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00021 * GNU General Public License for more details. 00022 00023 ******************************************************************************/ 00024 00025 #include <stdio.h> /* needed here for ifdef/else */ 00026 #include <stdlib.h> 00027 #include <string.h> 00028 #include <math.h> 00029 00030 00031 /******** 00032 ******** only compile this LAPACK/BLAS wrapper file if g2c.h is present! 00033 ********/ 00034 #include <grass/config.h> 00035 00036 #if defined(HAVE_LIBLAPACK) && defined(HAVE_LIBBLAS) && defined(HAVE_G2C_H) 00037 00038 #include <grass/gis.h> 00039 #include <grass/glocale.h> 00040 #include <grass/la.h> 00041 00042 00043 static int egcmp(const void *pa, const void *pb); 00044 00045 00059 mat_struct *G_matrix_init(int rows, int cols, int ldim) 00060 { 00061 mat_struct *tmp_arry; 00062 00063 if (rows < 1 || cols < 1 || ldim < rows) { 00064 G_warning(_("Matrix dimensions out of range")); 00065 return NULL; 00066 } 00067 00068 tmp_arry = (mat_struct *) G_malloc(sizeof(mat_struct)); 00069 tmp_arry->rows = rows; 00070 tmp_arry->cols = cols; 00071 tmp_arry->ldim = ldim; 00072 tmp_arry->type = MATRIX_; 00073 tmp_arry->v_indx = -1; 00074 00075 tmp_arry->vals = (doublereal *) G_calloc(ldim * cols, sizeof(doublereal)); 00076 tmp_arry->is_init = 1; 00077 00078 return tmp_arry; 00079 } 00080 00081 00091 int G_matrix_zero(mat_struct * A) 00092 { 00093 if (!A->vals) 00094 return 0; 00095 00096 memset(A->vals, 0, sizeof(A->vals)); 00097 00098 return 1; 00099 } 00100 00101 00117 int G_matrix_set(mat_struct * A, int rows, int cols, int ldim) 00118 { 00119 if (rows < 1 || cols < 1 || ldim < 0) { 00120 G_warning(_("Matrix dimensions out of range")); 00121 return -1; 00122 } 00123 00124 A->rows = rows; 00125 A->cols = cols; 00126 A->ldim = ldim; 00127 A->type = MATRIX_; 00128 A->v_indx = -1; 00129 00130 A->vals = (doublereal *) G_calloc(ldim * cols, sizeof(doublereal)); 00131 A->is_init = 1; 00132 00133 return 0; 00134 } 00135 00136 00148 mat_struct *G_matrix_copy(const mat_struct * A) 00149 { 00150 mat_struct *B; 00151 00152 if (!A->is_init) { 00153 G_warning(_("Matrix is not initialised fully.")); 00154 return NULL; 00155 } 00156 00157 if ((B = G_matrix_init(A->rows, A->cols, A->ldim)) == NULL) { 00158 G_warning(_("Unable to allocate space for matrix copy")); 00159 return NULL; 00160 } 00161 00162 memcpy(&B->vals[0], &A->vals[0], A->cols * A->ldim * sizeof(doublereal)); 00163 00164 return B; 00165 } 00166 00167 00181 mat_struct *G_matrix_add(mat_struct * mt1, mat_struct * mt2) 00182 { 00183 return G__matrix_add(mt1, mt2, 1, 1); 00184 } 00185 00186 00200 mat_struct *G_matrix_subtract(mat_struct * mt1, mat_struct * mt2) 00201 { 00202 return G__matrix_add(mt1, mt2, 1, -1); 00203 } 00204 00205 00219 mat_struct *G_matrix_scale(mat_struct * mt1, const double c) 00220 { 00221 return G__matrix_add(mt1, NULL, c, 0); 00222 } 00223 00224 00240 mat_struct *G__matrix_add(mat_struct * mt1, mat_struct * mt2, const double c1, 00241 const double c2) 00242 { 00243 mat_struct *mt3; 00244 int i, j; /* loop variables */ 00245 00246 if (c1 == 0) { 00247 G_warning(_("First scalar multiplier must be non-zero")); 00248 return NULL; 00249 } 00250 00251 if (c2 == 0) { 00252 if (!mt1->is_init) { 00253 G_warning(_("One or both input matrices uninitialised")); 00254 return NULL; 00255 } 00256 } 00257 00258 else { 00259 00260 if (!((mt1->is_init) && (mt2->is_init))) { 00261 G_warning(_("One or both input matrices uninitialised")); 00262 return NULL; 00263 } 00264 00265 if (mt1->rows != mt2->rows || mt1->cols != mt2->cols) { 00266 G_warning(_("Matrix order does not match")); 00267 return NULL; 00268 } 00269 } 00270 00271 if ((mt3 = G_matrix_init(mt1->rows, mt1->cols, mt1->ldim)) == NULL) { 00272 G_warning(_("Unable to allocate space for matrix sum")); 00273 return NULL; 00274 } 00275 00276 if (c2 == 0) { 00277 00278 for (i = 0; i < mt3->rows; i++) { 00279 for (j = 0; j < mt3->cols; j++) { 00280 mt3->vals[i + mt3->ldim * j] = 00281 c1 * mt1->vals[i + mt1->ldim * j]; 00282 } 00283 } 00284 } 00285 00286 else { 00287 00288 for (i = 0; i < mt3->rows; i++) { 00289 for (j = 0; j < mt3->cols; j++) { 00290 mt3->vals[i + mt3->ldim * j] = 00291 c1 * mt1->vals[i + mt1->ldim * j] + c2 * mt2->vals[i + 00292 mt2-> 00293 ldim * 00294 j]; 00295 } 00296 } 00297 } 00298 00299 return mt3; 00300 } 00301 00302 00303 #if defined(HAVE_LIBBLAS) 00304 00318 mat_struct *G_matrix_product(mat_struct * mt1, mat_struct * mt2) 00319 { 00320 mat_struct *mt3; 00321 doublereal unity = 1, zero = 0; 00322 integer rows, cols, interdim, lda, ldb; 00323 integer1 no_trans = 'n'; 00324 00325 if (!((mt1->is_init) || (mt2->is_init))) { 00326 G_warning(_("One or both input matrices uninitialised")); 00327 return NULL; 00328 } 00329 00330 if (mt1->cols != mt2->rows) { 00331 G_warning(_("Matrix order does not match")); 00332 return NULL; 00333 } 00334 00335 if ((mt3 = G_matrix_init(mt1->rows, mt2->cols, mt1->ldim)) == NULL) { 00336 G_warning(_("Unable to allocate space for matrix product")); 00337 return NULL; 00338 } 00339 00340 /* Call the driver */ 00341 00342 rows = (integer) mt1->rows; 00343 interdim = (integer) mt1->cols; 00344 cols = (integer) mt2->cols; 00345 00346 lda = (integer) mt1->ldim; 00347 ldb = (integer) mt2->ldim; 00348 00349 f77_dgemm(&no_trans, &no_trans, &rows, &cols, &interdim, &unity, 00350 mt1->vals, &lda, mt2->vals, &ldb, &zero, mt3->vals, &lda); 00351 00352 return mt3; 00353 } 00354 00355 #else /* defined(HAVE_LIBBLAS) */ 00356 #warning G_matrix_product() not compiled; requires BLAS library 00357 #endif /* defined(HAVE_LIBBLAS) */ 00358 00372 mat_struct *G_matrix_transpose(mat_struct * mt) 00373 { 00374 mat_struct *mt1; 00375 int ldim, ldo; 00376 doublereal *dbo, *dbt, *dbx, *dby; 00377 int cnt, cnt2; 00378 00379 /* Word align the workspace blocks */ 00380 if (mt->cols % 2 == 0) 00381 ldim = mt->cols; 00382 else 00383 ldim = mt->cols + 1; 00384 00385 mt1 = G_matrix_init(mt->cols, mt->rows, ldim); 00386 00387 /* Set initial values for reading arrays */ 00388 dbo = &mt->vals[0]; 00389 dbt = &mt1->vals[0]; 00390 ldo = mt->ldim; 00391 00392 for (cnt = 0; cnt < mt->cols; cnt++) { 00393 dbx = dbo; 00394 dby = dbt; 00395 00396 for (cnt2 = 0; cnt2 < ldo - 1; cnt2++) { 00397 *dby = *dbx; 00398 dby += ldim; 00399 dbx++; 00400 } 00401 00402 *dby = *dbx; 00403 00404 if (cnt < mt->cols - 1) { 00405 dbo += ldo; 00406 dbt++; 00407 } 00408 } 00409 00410 return mt1; 00411 } 00412 00413 00414 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) 00415 00439 /*** NOT YET COMPLETE: only some solutions' options available ***/ 00440 00441 int 00442 G_matrix_LU_solve(const mat_struct * mt1, mat_struct ** xmat0, 00443 const mat_struct * bmat, mat_type mtype) 00444 { 00445 mat_struct *wmat, *xmat, *mtx; 00446 00447 if (mt1->is_init == 0 || bmat->is_init == 0) { 00448 G_warning(_("Input: one or both data matrices uninitialised")); 00449 return -1; 00450 } 00451 00452 if (mt1->rows != mt1->cols || mt1->rows < 1) { 00453 G_warning(_("Principal matrix is not properly dimensioned")); 00454 return -1; 00455 } 00456 00457 if (bmat->cols < 1) { 00458 G_warning(_("Input: you must have at least one array to solve")); 00459 return -1; 00460 } 00461 00462 /* Now create solution matrix by copying the original coefficient matrix */ 00463 if ((xmat = G_matrix_copy(bmat)) == NULL) { 00464 G_warning(_("Could not allocate space for solution matrix")); 00465 return -1; 00466 } 00467 00468 /* Create working matrix for the coefficient array */ 00469 if ((mtx = G_matrix_copy(mt1)) == NULL) { 00470 G_warning(_("Could not allocate space for working matrix")); 00471 return -1; 00472 } 00473 00474 /* Copy the contents of the data matrix, to preserve the 00475 original information 00476 */ 00477 if ((wmat = G_matrix_copy(bmat)) == NULL) { 00478 G_warning(_("Could not allocate space for working matrix")); 00479 return -1; 00480 } 00481 00482 /* Now call appropriate LA driver to solve equations */ 00483 switch (mtype) { 00484 00485 case NONSYM: 00486 { 00487 integer *perm, res_info; 00488 integer num_eqns, nrhs, lda, ldb; 00489 00490 perm = (integer *) G_malloc(wmat->rows); 00491 00492 /* Set fields to pass to fortran routine */ 00493 num_eqns = (integer) mt1->rows; 00494 nrhs = (integer) wmat->cols; 00495 lda = (integer) mt1->ldim; 00496 ldb = (integer) wmat->ldim; 00497 00498 /* Call LA driver */ 00499 f77_dgesv(&num_eqns, &nrhs, mtx->vals, &lda, perm, wmat->vals, 00500 &ldb, &res_info); 00501 00502 /* Copy the results from the modified data matrix, taking account 00503 of pivot permutations ??? 00504 */ 00505 00506 /* 00507 for(indx1 = 0; indx1 < num_eqns; indx1++) { 00508 iperm = perm[indx1]; 00509 ptin = &wmat->vals[0] + indx1; 00510 ptout = &xmat->vals[0] + iperm; 00511 00512 for(indx2 = 0; indx2 < nrhs - 1; indx2++) { 00513 *ptout = *ptin; 00514 ptin += wmat->ldim; 00515 ptout += xmat->ldim; 00516 } 00517 00518 *ptout = *ptin; 00519 } 00520 */ 00521 00522 memcpy(xmat->vals, wmat->vals, 00523 wmat->cols * wmat->ldim * sizeof(doublereal)); 00524 00525 /* Free temp arrays */ 00526 G_free(perm); 00527 G_matrix_free(wmat); 00528 G_matrix_free(mtx); 00529 00530 if (res_info > 0) { 00531 G_warning(_("Matrix (or submatrix is singular). Solution undetermined")); 00532 return 1; 00533 } 00534 else if (res_info < 0) { 00535 G_warning(_("Problem in LA routine.")); 00536 return -1; 00537 } 00538 break; 00539 } 00540 default: 00541 { 00542 G_warning(_("Procedure not yet available for selected matrix type")); 00543 return -1; 00544 } 00545 } /* end switch */ 00546 00547 *xmat0 = xmat; 00548 00549 return 0; 00550 } 00551 00552 #else /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */ 00553 #warning G_matrix_LU_solve() not compiled; requires BLAS and LAPACK libraries 00554 #endif /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */ 00555 00556 #if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) 00557 00570 mat_struct *G_matrix_inverse(mat_struct * mt) 00571 { 00572 mat_struct *mt0, *res; 00573 int i, j, k; /* loop */ 00574 00575 if (mt->rows != mt->cols) { 00576 G_warning(_("Matrix is not square. Cannot determine inverse")); 00577 return NULL; 00578 } 00579 00580 if ((mt0 = G_matrix_init(mt->rows, mt->rows, mt->ldim)) == NULL) { 00581 G_warning(_("Unable to allocate space for matrix")); 00582 return NULL; 00583 } 00584 00585 /* Set `B' matrix to unit matrix */ 00586 for (i = 0; i < mt0->rows - 1; i++) { 00587 mt0->vals[i + i * mt0->ldim] = 1.0; 00588 00589 for (j = i + 1; j < mt0->cols; j++) { 00590 mt0->vals[i + j * mt0->ldim] = mt0->vals[j + i * mt0->ldim] = 0.0; 00591 } 00592 } 00593 00594 mt0->vals[mt0->rows - 1 + (mt0->rows - 1) * mt0->ldim] = 1.0; 00595 00596 /* Solve system */ 00597 if ((k = G_matrix_LU_solve(mt, &res, mt0, NONSYM)) == 1) { 00598 G_warning(_("Matrix is singular")); 00599 G_matrix_free(mt0); 00600 return NULL; 00601 } 00602 else if (k < 0) { 00603 G_warning(_("Problem in LA procedure.")); 00604 G_matrix_free(mt0); 00605 return NULL; 00606 } 00607 else { 00608 G_matrix_free(mt0); 00609 return res; 00610 } 00611 } 00612 00613 #else /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */ 00614 #warning G_matrix_inverse() not compiled; requires BLAS and LAPACK libraries 00615 #endif /* defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK) */ 00616 00617 00629 void G_matrix_free(mat_struct * mt) 00630 { 00631 if (mt->is_init) 00632 G_free(mt->vals); 00633 00634 G_free(mt); 00635 } 00636 00637 00649 void G_matrix_print(mat_struct * mt) 00650 { 00651 int i, j; 00652 char buf[64], numbuf[64]; 00653 00654 for (i = 0; i < mt->rows; i++) { 00655 strcpy(buf, ""); 00656 00657 for (j = 0; j < mt->cols; j++) { 00658 00659 sprintf(numbuf, "%14.6f", G_matrix_get_element(mt, i, j)); 00660 strcat(buf, numbuf); 00661 if (j < mt->cols - 1) 00662 strcat(buf, ", "); 00663 } 00664 00665 G_message("%s", buf); 00666 } 00667 00668 fprintf(stderr, "\n"); 00669 } 00670 00671 00688 int G_matrix_set_element(mat_struct * mt, int rowval, int colval, double val) 00689 { 00690 if (!mt->is_init) { 00691 G_warning(_("Element array has not been allocated")); 00692 return -1; 00693 } 00694 00695 if (rowval >= mt->rows || colval >= mt->cols || rowval < 0 || colval < 0) { 00696 G_warning(_("Specified element is outside array bounds")); 00697 return -1; 00698 } 00699 00700 mt->vals[rowval + colval * mt->ldim] = (doublereal) val; 00701 00702 return 0; 00703 } 00704 00705 00721 double G_matrix_get_element(mat_struct * mt, int rowval, int colval) 00722 { 00723 double val; 00724 00725 /* Should do some checks, but this would require an error control 00726 system: later? */ 00727 00728 return (val = (double)mt->vals[rowval + colval * mt->ldim]); 00729 } 00730 00731 00744 vec_struct *G_matvect_get_column(mat_struct * mt, int col) 00745 { 00746 int i; /* loop */ 00747 vec_struct *vc1; 00748 00749 if (col < 0 || col >= mt->cols) { 00750 G_warning(_("Specified matrix column index is outside range")); 00751 return NULL; 00752 } 00753 00754 if (!mt->is_init) { 00755 G_warning(_("Matrix is not initialised")); 00756 return NULL; 00757 } 00758 00759 if ((vc1 = G_vector_init(mt->rows, mt->ldim, CVEC)) == NULL) { 00760 G_warning(_("Could not allocate space for vector structure")); 00761 return NULL; 00762 } 00763 00764 for (i = 0; i < mt->rows; i++) 00765 G_matrix_set_element((mat_struct *) vc1, i, 0, 00766 G_matrix_get_element(mt, i, col)); 00767 00768 return vc1; 00769 } 00770 00771 00785 vec_struct *G_matvect_get_row(mat_struct * mt, int row) 00786 { 00787 int i; /* loop */ 00788 vec_struct *vc1; 00789 00790 if (row < 0 || row >= mt->cols) { 00791 G_warning(_("Specified matrix row index is outside range")); 00792 return NULL; 00793 } 00794 00795 if (!mt->is_init) { 00796 G_warning(_("Matrix is not initialised")); 00797 return NULL; 00798 } 00799 00800 if ((vc1 = G_vector_init(mt->cols, mt->ldim, RVEC)) == NULL) { 00801 G_warning(_("Could not allocate space for vector structure")); 00802 return NULL; 00803 } 00804 00805 for (i = 0; i < mt->cols; i++) 00806 G_matrix_set_element((mat_struct *) vc1, 0, i, 00807 G_matrix_get_element(mt, row, i)); 00808 00809 return vc1; 00810 } 00811 00812 00828 int G_matvect_extract_vector(mat_struct * mt, vtype vt, int indx) 00829 { 00830 if (vt == RVEC && indx >= mt->rows) { 00831 G_warning(_("Specified row index is outside range")); 00832 return -1; 00833 } 00834 00835 else if (vt == CVEC && indx >= mt->cols) { 00836 G_warning(_("Specified column index is outside range")); 00837 return -1; 00838 } 00839 00840 switch (vt) { 00841 00842 case RVEC: 00843 { 00844 mt->type = ROWVEC_; 00845 mt->v_indx = indx; 00846 } 00847 00848 case CVEC: 00849 { 00850 mt->type = COLVEC_; 00851 mt->v_indx = indx; 00852 } 00853 00854 default: 00855 { 00856 G_warning(_("Unknown vector type.")); 00857 return -1; 00858 } 00859 00860 } 00861 00862 return 0; 00863 00864 } 00865 00866 00878 int G_matvect_retrieve_matrix(vec_struct * vc) 00879 { 00880 /* We have to take the integrity of the vector structure 00881 largely on trust 00882 */ 00883 00884 vc->type = MATRIX_; 00885 vc->v_indx = -1; 00886 00887 return 0; 00888 } 00889 00890 00905 vec_struct *G_vector_init(int cells, int ldim, vtype vt) 00906 { 00907 vec_struct *tmp_arry; 00908 00909 if ((cells < 1) || (vt == RVEC && ldim < 1) 00910 || (vt == CVEC && ldim < cells) || ldim < 0) { 00911 G_warning("Vector dimensions out of range."); 00912 return NULL; 00913 } 00914 00915 tmp_arry = (vec_struct *) G_malloc(sizeof(vec_struct)); 00916 00917 if (vt == RVEC) { 00918 tmp_arry->rows = 1; 00919 tmp_arry->cols = cells; 00920 tmp_arry->ldim = ldim; 00921 tmp_arry->type = ROWVEC_; 00922 } 00923 00924 else if (vt == CVEC) { 00925 tmp_arry->rows = cells; 00926 tmp_arry->cols = 1; 00927 tmp_arry->ldim = ldim; 00928 tmp_arry->type = COLVEC_; 00929 } 00930 00931 tmp_arry->v_indx = 0; 00932 00933 tmp_arry->vals = (doublereal *) G_calloc(ldim * tmp_arry->cols, 00934 sizeof(doublereal)); 00935 tmp_arry->is_init = 1; 00936 00937 return tmp_arry; 00938 } 00939 00940 00952 void G_vector_free(vec_struct * v) 00953 { 00954 if (v->is_init) 00955 G_free(v->vals); 00956 00957 G_free(v); 00958 } 00959 00960 00975 vec_struct *G_vector_sub(vec_struct * v1, vec_struct * v2, vec_struct * out) 00976 { 00977 int idx1, idx2, idx0; 00978 int i; 00979 00980 if (!out->is_init) { 00981 G_warning(_("Output vector is uninitialized")); 00982 return NULL; 00983 } 00984 00985 if (v1->type != v2->type) { 00986 G_warning(_("Vectors are not of the same type")); 00987 return NULL; 00988 } 00989 00990 if (v1->type != out->type) { 00991 G_warning(_("Output vector is of incorrect type")); 00992 return NULL; 00993 } 00994 00995 if (v1->type == MATRIX_) { 00996 G_warning(_("Matrices not allowed")); 00997 return NULL; 00998 } 00999 01000 if ((v1->type == ROWVEC_ && v1->cols != v2->cols) || 01001 (v1->type == COLVEC_ && v1->rows != v2->rows)) { 01002 G_warning(_("Vectors have differing dimensions")); 01003 return NULL; 01004 } 01005 01006 if ((v1->type == ROWVEC_ && v1->cols != out->cols) || 01007 (v1->type == COLVEC_ && v1->rows != out->rows)) { 01008 G_warning(_("Output vector has incorrect dimension")); 01009 return NULL; 01010 } 01011 01012 idx1 = (v1->v_indx > 0) ? v1->v_indx : 0; 01013 idx2 = (v2->v_indx > 0) ? v2->v_indx : 0; 01014 idx0 = (out->v_indx > 0) ? out->v_indx : 0; 01015 01016 if (v1->type == ROWVEC_) { 01017 for (i = 0; i < v1->cols; i++) 01018 G_matrix_set_element(out, idx0, i, 01019 G_matrix_get_element(v1, idx1, i) - 01020 G_matrix_get_element(v2, idx2, i)); 01021 } 01022 else { 01023 for (i = 0; i < v1->rows; i++) 01024 G_matrix_set_element(out, i, idx0, 01025 G_matrix_get_element(v1, i, idx1) - 01026 G_matrix_get_element(v2, i, idx2)); 01027 } 01028 01029 return out; 01030 } 01031 01032 01050 int G_vector_set(vec_struct * A, int cells, int ldim, vtype vt, int vindx) 01051 { 01052 if ((cells < 1) || (vt == RVEC && ldim < 1) 01053 || (vt == CVEC && ldim < cells) || ldim < 0) { 01054 G_warning(_("Vector dimensions out of range")); 01055 return -1; 01056 } 01057 01058 if ((vt == RVEC && vindx >= A->cols) || (vt == CVEC && vindx >= A->rows)) { 01059 G_warning(_("Row/column out of range")); 01060 return -1; 01061 } 01062 01063 if (vt == RVEC) { 01064 A->rows = 1; 01065 A->cols = cells; 01066 A->ldim = ldim; 01067 A->type = ROWVEC_; 01068 } 01069 else { 01070 A->rows = cells; 01071 A->cols = 1; 01072 A->ldim = ldim; 01073 A->type = COLVEC_; 01074 } 01075 01076 if (vindx < 0) 01077 A->v_indx = 0; 01078 else 01079 A->v_indx = vindx; 01080 01081 A->vals = (doublereal *) G_calloc(ldim * A->cols, sizeof(doublereal)); 01082 A->is_init = 1; 01083 01084 return 0; 01085 01086 } 01087 01088 01089 #if defined(HAVE_LIBBLAS) 01090 01103 double G_vector_norm_euclid(vec_struct * vc) 01104 { 01105 integer incr, Nval; 01106 doublereal *startpt; 01107 01108 if (!vc->is_init) 01109 G_fatal_error(_("Matrix is not initialised")); 01110 01111 if (vc->type == ROWVEC_) { 01112 Nval = (integer) vc->cols; 01113 incr = (integer) vc->ldim; 01114 if (vc->v_indx < 0) 01115 startpt = vc->vals; 01116 else 01117 startpt = vc->vals + vc->v_indx; 01118 } 01119 else { 01120 Nval = (integer) vc->rows; 01121 incr = 1; 01122 if (vc->v_indx < 0) 01123 startpt = vc->vals; 01124 else 01125 startpt = vc->vals + vc->v_indx * vc->ldim; 01126 } 01127 01128 /* Call the BLAS routine dnrm2_() */ 01129 return (double)f77_dnrm2(&Nval, startpt, &incr); 01130 } 01131 01132 #else /* defined(HAVE_LIBBLAS) */ 01133 #warning G_vector_norm_euclid() not compiled; requires BLAS library 01134 #endif /* defined(HAVE_LIBBLAS) */ 01135 01136 01154 double G_vector_norm_maxval(vec_struct * vc, int vflag) 01155 { 01156 doublereal xval, *startpt, *curpt; 01157 double cellval; 01158 int ncells, incr; 01159 01160 if (!vc->is_init) 01161 G_fatal_error(_("Matrix is not initialised")); 01162 01163 if (vc->type == ROWVEC_) { 01164 ncells = (integer) vc->cols; 01165 incr = (integer) vc->ldim; 01166 if (vc->v_indx < 0) 01167 startpt = vc->vals; 01168 else 01169 startpt = vc->vals + vc->v_indx; 01170 } 01171 else { 01172 ncells = (integer) vc->rows; 01173 incr = 1; 01174 if (vc->v_indx < 0) 01175 startpt = vc->vals; 01176 else 01177 startpt = vc->vals + vc->v_indx * vc->ldim; 01178 } 01179 01180 xval = *startpt; 01181 curpt = startpt; 01182 01183 while (ncells > 0) { 01184 if (curpt != startpt) { 01185 switch (vflag) { 01186 01187 case MAX_POS: 01188 { 01189 if (*curpt > xval) 01190 xval = *curpt; 01191 break; 01192 } 01193 01194 case MAX_NEG: 01195 { 01196 if (*curpt < xval) 01197 xval = *curpt; 01198 break; 01199 } 01200 01201 case MAX_ABS: 01202 { 01203 cellval = (double)(*curpt); 01204 if (hypot(cellval, cellval) > (double)xval) 01205 xval = *curpt; 01206 } 01207 } /* switch */ 01208 } /* if(curpt != startpt) */ 01209 01210 curpt += incr; 01211 ncells--; 01212 } 01213 01214 return (double)xval; 01215 } 01216 01217 01229 double G_vector_norm1(vec_struct * vc) 01230 { 01231 double result = 0.0; 01232 int idx; 01233 int i; 01234 01235 if (!vc->is_init) { 01236 G_warning(_("Matrix is not initialised")); 01237 return 0.0 / 0.0; /* NaN */ 01238 } 01239 01240 idx = (vc->v_indx > 0) ? vc->v_indx : 0; 01241 01242 if (vc->type == ROWVEC_) { 01243 for (i = 0; i < vc->cols; i++) 01244 result += fabs(G_matrix_get_element(vc, idx, i)); 01245 } 01246 else { 01247 for (i = 0; i < vc->rows; i++) 01248 result += fabs(G_matrix_get_element(vc, i, idx)); 01249 } 01250 01251 return result; 01252 } 01253 01254 01266 vec_struct *G_vector_copy(const vec_struct * vc1, int comp_flag) 01267 { 01268 vec_struct *tmp_arry; 01269 int incr1, incr2; 01270 doublereal *startpt1, *startpt2, *curpt1, *curpt2; 01271 int cnt; 01272 01273 if (!vc1->is_init) { 01274 G_warning(_("Vector structure is not initialised")); 01275 return NULL; 01276 } 01277 01278 tmp_arry = (vec_struct *) G_malloc(sizeof(vec_struct)); 01279 01280 if (comp_flag == DO_COMPACT) { 01281 if (vc1->type == ROWVEC_) { 01282 tmp_arry->rows = 1; 01283 tmp_arry->cols = vc1->cols; 01284 tmp_arry->ldim = 1; 01285 tmp_arry->type = ROWVEC_; 01286 tmp_arry->v_indx = 0; 01287 } 01288 else if (vc1->type == COLVEC_) { 01289 tmp_arry->rows = vc1->rows; 01290 tmp_arry->cols = 1; 01291 tmp_arry->ldim = vc1->ldim; 01292 tmp_arry->type = COLVEC_; 01293 tmp_arry->v_indx = 0; 01294 } 01295 else { 01296 G_warning("Type is not vector."); 01297 return NULL; 01298 } 01299 } 01300 else if (comp_flag == NO_COMPACT) { 01301 tmp_arry->v_indx = vc1->v_indx; 01302 tmp_arry->rows = vc1->rows; 01303 tmp_arry->cols = vc1->cols; 01304 tmp_arry->ldim = vc1->ldim; 01305 tmp_arry->type = vc1->type; 01306 } 01307 else { 01308 G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n"); 01309 return NULL; 01310 } 01311 01312 tmp_arry->vals = (doublereal *) G_calloc(tmp_arry->ldim * tmp_arry->cols, 01313 sizeof(doublereal)); 01314 if (comp_flag == DO_COMPACT) { 01315 if (tmp_arry->type == ROWVEC_) { 01316 startpt1 = tmp_arry->vals; 01317 startpt2 = vc1->vals + vc1->v_indx; 01318 curpt1 = startpt1; 01319 curpt2 = startpt2; 01320 incr1 = 1; 01321 incr2 = vc1->ldim; 01322 cnt = vc1->cols; 01323 } 01324 else if (tmp_arry->type == COLVEC_) { 01325 startpt1 = tmp_arry->vals; 01326 startpt2 = vc1->vals + vc1->v_indx * vc1->ldim; 01327 curpt1 = startpt1; 01328 curpt2 = startpt2; 01329 incr1 = 1; 01330 incr2 = 1; 01331 cnt = vc1->rows; 01332 } 01333 else { 01334 G_warning("Structure type is not vector."); 01335 return NULL; 01336 } 01337 } 01338 else if (comp_flag == NO_COMPACT) { 01339 startpt1 = tmp_arry->vals; 01340 startpt2 = vc1->vals; 01341 curpt1 = startpt1; 01342 curpt2 = startpt2; 01343 incr1 = 1; 01344 incr2 = 1; 01345 cnt = vc1->ldim * vc1->cols; 01346 } 01347 else { 01348 G_warning("Copy method must be specified: [DO,NO]_COMPACT.\n"); 01349 return NULL; 01350 } 01351 01352 while (cnt > 0) { 01353 memcpy(curpt1, curpt2, sizeof(doublereal)); 01354 curpt1 += incr1; 01355 curpt2 += incr2; 01356 cnt--; 01357 } 01358 01359 tmp_arry->is_init = 1; 01360 01361 return tmp_arry; 01362 } 01363 01364 01379 int G_matrix_read(FILE * fp, mat_struct * out) 01380 { 01381 char buff[100]; 01382 int rows, cols; 01383 int i, j, row; 01384 double val; 01385 01386 /* skip comments */ 01387 for (;;) { 01388 if (!G_getl(buff, sizeof(buff), fp)) 01389 return -1; 01390 if (buff[0] != '#') 01391 break; 01392 } 01393 01394 if (sscanf(buff, "Matrix: %d by %d", &rows, &cols) != 2) { 01395 G_warning(_("Input format error")); 01396 return -1; 01397 } 01398 01399 G_matrix_set(out, rows, cols, rows); 01400 01401 for (i = 0; i < rows; i++) { 01402 if (fscanf(fp, "row%d:", &row) != 1 || row != i) { 01403 G_warning(_("Input format error")); 01404 return -1; 01405 } 01406 for (j = 0; j < cols; j++) { 01407 if (fscanf(fp, "%lf:", &val) != 1) { 01408 G_warning(_("Input format error")); 01409 return -1; 01410 } 01411 01412 G_matrix_set_element(out, i, j, val); 01413 } 01414 } 01415 01416 return 0; 01417 } 01418 01419 01433 int G_matrix_stdin(mat_struct * out) 01434 { 01435 return G_matrix_read(stdin, out); 01436 } 01437 01438 01451 int G_matrix_eigen_sort(vec_struct * d, mat_struct * m) 01452 { 01453 mat_struct tmp; 01454 int i, j; 01455 int idx; 01456 01457 G_matrix_set(&tmp, m->rows + 1, m->cols, m->ldim + 1); 01458 01459 idx = (d->v_indx > 0) ? d->v_indx : 0; 01460 01461 /* concatenate (vertically) m and d into tmp */ 01462 for (i = 0; i < m->cols; i++) { 01463 for (j = 0; j < m->rows; j++) 01464 G_matrix_set_element(&tmp, j + 1, i, 01465 G_matrix_get_element(m, j, i)); 01466 if (d->type == ROWVEC_) 01467 G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, idx, i)); 01468 else 01469 G_matrix_set_element(&tmp, 0, i, G_matrix_get_element(d, i, idx)); 01470 } 01471 01472 /* sort the combined matrix */ 01473 qsort(tmp.vals, tmp.cols, tmp.ldim * sizeof(doublereal), egcmp); 01474 01475 /* split tmp into m and d */ 01476 for (i = 0; i < m->cols; i++) { 01477 for (j = 0; j < m->rows; j++) 01478 G_matrix_set_element(m, j, i, 01479 G_matrix_get_element(&tmp, j + 1, i)); 01480 if (d->type == ROWVEC_) 01481 G_matrix_set_element(d, idx, i, G_matrix_get_element(&tmp, 0, i)); 01482 else 01483 G_matrix_set_element(d, i, idx, G_matrix_get_element(&tmp, 0, i)); 01484 } 01485 01486 G_free(tmp.vals); 01487 01488 return 0; 01489 } 01490 01491 01492 static int egcmp(const void *pa, const void *pb) 01493 { 01494 double a = *(doublereal * const)pa; 01495 double b = *(doublereal * const)pb; 01496 01497 if (a > b) 01498 return 1; 01499 if (a < b) 01500 return -1; 01501 01502 return 0; 01503 } 01504 01505 01506 #endif /* HAVE_BLAS && HAVE_LAPACK && HAVE_G2C */