GRASS Programmer's Manual  6.4.2(2012)
la.c
Go to the documentation of this file.
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 */
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines