Changeset 17899


Ignore:
Timestamp:
04/30/14 16:12:35 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: m1qn3 now workingsvn diff cores/controlm1qn3_core.cpp

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp

    r17897 r17899  
    33 */
    44
     5#include <config.h>
    56#include "./cores.h"
    67#include "../toolkits/toolkits.h"
    … …  
    1516extern "C" void *ctcabe_; // DIS mode : Conversion
    1617extern "C" void *euclid_; // Scalar product
    17 typedef void (*SimulFunc) (long*, long *, double [], double *, double [], long [], float [], double []);
    18 extern "C" void m1qn3_ (void f(long*, long *, double [], double *, double [], long [], float [], double []),
     18typedef void (*SimulFunc) (long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs);
     19extern "C" void m1qn3_ (void f(long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs),
    1920                        void **, void **, void **,
    2021                        long *, double [], double *, double [], double*, double *,
    2122                        double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *,
    22                         long *, long *, long [], float [], double []
    23                         );
     23                        long *, long *, long [], float [],void* );
    2424
    25 void simul(long* indic,long* n,double x[2],double* f,double g[2],long izs[1],float rzs[1],double dzs[1]){
    26         /*minimization of x^2+cy^4 where c>0 is a parameter*/
    27         double c=dzs[0];
    28         *f=x[0]*x[0]+c*x[1]*x[1]*x[1]*x[1];
    29         g[0]=2.*x[0];
    30         g[1]=4.*c*x[1]*x[1]*x[1];
    31 }
     25/*Cost function prototype*/
     26void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
    3227
    3328void controlm1qn3_core(FemModel* femmodel){
    3429
    3530        /*Intermediaries*/
    36         long   omode; /*m1qn3 output flag*/
    37         double f;     /*cost function value*/
     31        long         omode;
     32        double       f;
     33        int          nsteps,maxiter;
     34        int          intn,num_controls,solution_type;
     35        IssmDouble  *X  = NULL;
     36        IssmDouble  *G  = NULL;
     37        IssmDouble  *XL = NULL;
     38        IssmDouble  *XU = NULL;
    3839
     40        /*Recover some parameters*/
     41        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     42        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     43        femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
     44        femmodel->parameters->SetParam(false,SaveResultsEnum);
     45        maxiter=nsteps*10;
     46
     47        /*Initialize M1QN3 parameters*/
    3948        if(VerboseControl())_printf0_("   Initialize M1QN3 parameters\n");
    4049        SimulFunc costfuncion  = &simul;    /*Cost function address*/
    4150        void**    prosca       = &euclid_;  /*Dot product function (euclid is the default)*/
    4251        char      normtype[3];              /*Norm type: dfn = scalar product defined by prosca*/ strcpy(normtype, "dfn");
    43         double    dzs[1];                   /*Arrays used by m1qn3 subroutines*/
    4452        long      izs[5];                   /*Arrays used by m1qn3 subroutines*/
     53        long      iz[5];                    /*Integer m1qn3 working array of size 5*/
    4554        float     rzs[1];                   /*Arrays used by m1qn3 subroutines*/
    4655        long      impres       = 0;         /*verbosity level*/
    … …  
    4857        long      indic        = 4;         /*compute f and g*/
    4958        long      reverse      = 0;         /*reverse or direct mode*/
    50         long      iz[5];                    /*Integer m1qn3 working array of size 5*/
    5159        long      io           = 6;         /*Channel number for the output*/
    5260
    5361        /*Optimization criterions*/
    54         double    dxmin        = 1.e-10;    /*Resolution for the solution x*/
    55         double    epsrel       = 1.e-5;     /*Gradient stopping criterion in ]0 1[ -> |gk|/|g1| < epsrel*/
    56         long      niter        = 200;       /*Maximum number of iterations*/
    57         long      nsim         = 200;       /*Maximum number of function calls*/
     62        double    dxmin        = 1.e-1;    /*Resolution for the solution x*/
     63        double    epsrel       = 1.e-4;     /*Gradient stopping criterion in ]0 1[ -> |gk|/|g1| < epsrel*/
     64        long      niter        = long(nsteps); /*Maximum number of iterations*/
     65        long      nsim         = long(maxiter);/*Maximum number of function calls*/
     66
     67        /*Get initial guess*/
     68        Vector<IssmDouble> *Xpetsc = NULL;
     69        GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
     70        X = Xpetsc->ToMPISerial();
     71        Xpetsc->GetSize(&intn);
     72        delete Xpetsc;
    5873
    5974        /*Get problem dimension and initialize gradient and initial guess*/
    60         long    n = 2;
    61         double* g = xNew<double>(n);
    62         double* x = xNew<double>(n);
    63         for(int i=0;i<n;i++) x[i]=5.;
    64         dzs[0] = 10; //c = 10 function parameter
     75        long n = long(intn);
     76        G = xNew<double>(n);
    6577
    6678        /*Allocate m1qn3 working arrays (see doc)*/
    … …  
    7082
    7183        if(VerboseControl())_printf0_("   Computing initial solution\n");
    72         simul(&indic,&n,x,&f,g,izs,rzs,&dzs[0]);
     84        _printf0_("\n");
     85        _printf0_("Cost function f(x)   |  List of contributions\n");
     86        _printf0_("_____________________________________________\n");
     87        indic = 0; //no adjoint required
     88        simul(&indic,&n,X,&f,G,izs,rzs,(void*)femmodel);
    7389        double f1=f;
    7490
     91        indic = 4; //adjoint and gradient required
    7592        m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_,
    76                                 &n,x,&f,g,&dxmin,&f1,
     93                                &n,X,&f,G,&dxmin,&f1,
    7794                                &epsrel,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
    78                                 &reverse,&indic,izs,rzs,dzs);
     95                                &reverse,&indic,izs,rzs,(void*)femmodel);
    7996
    8097        switch(int(omode)){
    … …  
    90107        }
    91108
    92         _printf0_(" == Final cost function = "<< f <<"\n");
    93         //_printf0_(" == Final x = ["<<x[0]<<" "<<x[1]<<"]\n");
     109        /*Get solution*/
     110        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     111        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
     112        femmodel->OutputControlsx(&femmodel->results);
     113        femmodel->results->AddObject(new GenericExternalResult<double>(JEnum,f,1,0));
     114
     115        /*Finalize*/
     116        if(VerboseControl()) _printf0_("   preparing final solution\n");
     117        femmodel->parameters->SetParam(true,SaveResultsEnum);
     118        void (*solutioncore)(FemModel*)=NULL;
     119        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     120        solutioncore(femmodel);
    94121
    95122        /*Clean-up and return*/
    96         xDelete<double>(g);
    97         xDelete<double>(x);
     123        xDelete<double>(G);
     124        xDelete<double>(X);
    98125        xDelete<double>(dz);
    99126}
     127
     128/*Cost function definition*/
     129void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){
     130
     131        /*Recover Femmodel*/
     132        double      f;
     133        int         intn,solution_type;
     134        FemModel   *femmodel  = (FemModel*)dzs;
     135
     136        /*Recover responses*/
     137        int         num_responses;
     138        int        *responses = NULL;
     139        femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
     140
     141        /*Update control input*/
     142        SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
     143
     144        /*Recover some parameters*/
     145        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     146
     147        /*Compute solution and adjoint*/
     148        void (*solutioncore)(FemModel*)=NULL;
     149        void (*adjointcore)(FemModel*)=NULL;
     150        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     151        solutioncore(femmodel);
     152
     153        /*Compute objective function*/
     154        femmodel->CostFunctionx(pf);
     155        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
     156
     157        /*Retrieve objective functions independently*/
     158        for(int i=0;i<num_responses;i++){
     159                femmodel->Responsex(&f,EnumToStringx(responses[i]));
     160                _printf0_(" "<<setw(12)<<setprecision(7)<<f);
     161        }
     162        _printf0_("\n");
     163
     164        if(indic==0) return; /*dry run, no gradient required*/
     165
     166        /*Compute Adjoint*/
     167        AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
     168        adjointcore(femmodel);
     169
     170        /*Compute gradient*/
     171        IssmDouble* G2 = NULL;
     172        Gradjx(&G2,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     173        for(long i=0;i<*n;i++) G[i] = -G2[i];
     174
     175        /*Clean-up and return*/
     176        xDelete<int>(responses);
     177        xDelete<IssmDouble>(G2);
     178}
     179
    100180#else
    101181void controlm1qn3_core(FemModel* femmodel){
    102182        _error_("M1QN3 not installed");
    103183}
    104 #endif //_HAVE_TAO_
     184#endif //_HAVE_M1QN3_

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

    r16167 r17899  
    6666        xDelete<int>(control_type);
    6767}
     68void Gradjx(IssmDouble** pgradient,IssmDouble** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
     69
     70        /*output: */
     71        IssmDouble* gradient=NULL;
     72
     73        /*intermediary: */
     74        Vector<IssmDouble>* vec_gradient=NULL;
     75
     76        Gradjx(&vec_gradient,pnorm_list,elements,nodes, vertices,loads,materials,parameters);
     77        gradient=vec_gradient->ToMPISerial();
     78
     79        /*Free ressources:*/
     80        delete vec_gradient;
     81
     82        /*Assign output pointers:*/
     83        *pgradient=gradient;
     84}

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

    r15000 r17899  
    1010/* local prototypes: */
    1111void Gradjx(Vector<IssmDouble>** pgrad_g,IssmDouble** pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters);
     12void Gradjx(IssmDouble** pgrad_g,IssmDouble** pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters);
    1213
    1314#endif  /* _GRADJX_H */
Note: See TracChangeset for help on using the changeset viewer.