Changeset 12183


Ignore:
Timestamp:
05/03/12 09:54:55 (13 years ago)
Author:
Mathieu Morlighem
Message:

added options and GSL support

Location:
issm/trunk-jpl/src/c/modules/Krigingx
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp

    r12177 r12183  
    99#include "../../objects/objects.h"
    1010#include "../modules.h"
     11#include <gsl/gsl_linalg.h>
    1112
    12 int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp){
     13int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options){
    1314
    1415        /*output*/
    … …  
    4950
    5051                /*Solve the three linear systems*/
    51                 MumpsSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
    52                 MumpsSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
    53                 MumpsSolve(&GinvZ, Gamma,observations,n_obs); // Gamma^-1 Z
     52                GslSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
     53                GslSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
     54                GslSolve(&GinvZ, Gamma,observations,n_obs); // Gamma^-1 Z
    5455
    5556                /*Prepare predictor*/
    … …  
    103104}
    104105
    105 void MumpsSolve(double** Xdouble,double* Adouble,double* Bdouble,int n){
    106 #ifdef _HAVE_PETSC_
     106void GslSolve(double** pX,double* A,double* B,int n){
     107#ifdef _HAVE_GSL_
    107108
    108         /*Intermediaries*/
    109         Mat A = NULL;
    110         Vec B = NULL;
    111         Vec X = NULL;
    112         KSP ksp = NULL;
    113         PC  pc = NULL;
     109        /*GSL Matrices and vectors: */
     110        int              s;
     111        gsl_matrix_view  a;
     112        gsl_vector_view  b;
     113        gsl_vector      *x = NULL;
     114        gsl_permutation *p = NULL;
    114115
    115         /*Create parameters and select MUMPS solver*/
    116         PetscOptionsSetFromOptions();
    117         PetscOptionsClear();
    118         PetscOptionsInsertMultipleString("-mat_type mpiaij \
    119                                 -ksp_type preonly \
    120                                 -pc_type lu \
    121                                 -pc_factor_mat_solver_package mumps \
    122                                 -mat_mumps_icntl_14 120 \
    123                                 -pc_factor_shift_positive_definite true");
     116        /*A will be modified by LU decomposition. Use copy*/
     117        double* Acopy = (double*)xmalloc(n*n*sizeof(double));
     118        memcpy(Acopy,A,n*n*sizeof(double));
    124119
    125         /*Create indexing and ones*/
    126         int* idx=(int*)xmalloc(n*sizeof(int));
    127         for(int i=0;i<n;i++) idx[i]=i;
     120        /*Initialize gsl matrices and vectors: */
     121        a = gsl_matrix_view_array (Acopy,n,n);
     122        b = gsl_vector_view_array (B,n);
     123        x = gsl_vector_alloc (n);
    128124
    129         /*Convert matrix to Dense PETSc matrix*/
    130         MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
    131         MatSetValues(A,n,idx,n,idx,Adouble,INSERT_VALUES);
    132         MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    133         MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
     125        /*Run LU and solve: */
     126        p = gsl_permutation_alloc (n);
     127        gsl_linalg_LU_decomp (&a.matrix, p, &s);
     128        gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
    134129
    135         /*Convert vector to PETSc vector*/
    136         VecCreateSeq(PETSC_COMM_SELF,n,&B);
    137         VecSetValues(B,n,idx,Bdouble,INSERT_VALUES);
    138         VecAssemblyBegin(B);
    139         VecAssemblyEnd(B);
     130        //printf ("x = \n");
     131        //gsl_vector_fprintf (stdout, x, "%g");
    140132
    141         /*Allocate output*/
    142         VecCreateSeq(PETSC_COMM_SELF,n,&X);
     133        /*Copy result*/
     134        double* X = (double*)xmalloc(n*sizeof(double));
     135        memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
    143136
    144         /*Go Solve*/
    145         KSPCreate(MPI_COMM_WORLD,&ksp);
    146         KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
    147         KSPSetFromOptions(ksp);
    148         KSPGetPC(ksp,&pc);
    149         PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
    150         KSPSolve(ksp,B,X);
    151         KSPFree(&ksp);
    152 
    153         /*Serialize vector*/
    154         VecToMPISerial(Xdouble,X);
     137        /*Clean up and assign output pointer*/
     138        xfree((void**)&Acopy);
     139        gsl_permutation_free(p);
     140        gsl_vector_free(x);
     141        *pX=X;
    155142#else
    156         error_("PETSc support required");
     143        error_("GSL support required");
    157144#endif
    158145}

  • issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h

    r12174 r12183  
    99#include "../../toolkits/toolkits.h"
    1010
    11 int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp);
     11int  Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options);
    1212void SemiVariogram(double* gamma,double x1,double y1,double x2,double y2);
    13 void MumpsSolve(double** X,double* A,double* B,int n); //Solve X for the linear system AX = B
     13void GslSolve(double** pX,double* A,double* B,int n);
    1414
    1515#endif /* _KRIGINGX_H */
Note: See TracChangeset for help on using the changeset viewer.