testgradluaf.c

sqp程序包。用sqp算法实现非线性约束的优化求解

#define  X
#include "o8gene.h"
#undef   X

/* **************************************************************************** */
/*                      testprogram for egradf,egradh,egradg                    */
/* **************************************************************************** */
void main(void) {

    void    setup0(void);
    void    setup (void);
    void    ef    (DOUBLE x[],DOUBLE *fx);
    void    egradf(DOUBLE x[],DOUBLE gradf[]);
    void    eh    (INTEGER i,DOUBLE x[],DOUBLE *hxi);
    void    egradh(INTEGER i,DOUBLE x[],DOUBLE gradhi[]);
    void    eg    (INTEGER i,DOUBLE x[],DOUBLE *gxi);
    void    egradg(INTEGER i,DOUBLE x[],DOUBLE gradgi[]);
    DOUBLE  vecnor(INTEGER n,DOUBLE x[]);
    void matdr2(DOUBLE a[][NRESM+1],INTEGER me,INTEGER ne,INTEGER ma,INTEGER na,
                char head[],FILE *lognum,LOGICAL fix,INTEGER nd1,INTEGER nd2);
    
    static INTEGER  i,j,k,countc;
    
    static DOUBLE   yy[NX+1],err[NX+1][NRESM+1],errtol,term;
    static DOUBLE   gnorm[NRESM+1];
    static DOUBLE   one = 1.;
    static LOGICAL  errfnd;
    static char     head[41];
    static FILE     *buf;
    
    epsmac=1.e-16;    
    setup0();
    
    nres = nh+ng;
    
    setup();
    
    countc = 0;
    
    printf("testprogram for egradf,egradh,egradg\n");
    printf("required input in standard format\n");
    printf("output in file testgrad.res\n");
    printf("input: stepsize for num. diff. = ");
    scanf ("%le", &epsx);
    printf("we check precision relative to the norm of \n");
    printf("the presumably correct analytic gradient\n");
    printf("input feasible error tolerance = ");
    scanf ("%le", &errtol);
    printf("output of true gradients desired? (1/0) ");
    scanf ("%i", &te1);
    
    buf = fopen("testgrad.res","w");
    
    L100:

    printf("input of x desired (1/0)(otherwise taken from setup) ");
    scanf ("%i", &inx);
    
    if ( inx ) {
        for (i = 1 ; i <= n ; i++) {
            printf("x[%3i] = ? ", i);
            scanf ("%le", &x[i]);
        }
    } else {
        if ( countc >= 1 ) {
            printf("x taken from setup\n");
            printf("results written in testgrad.res\n");
            fclose(buf);
            exit(0);
        }
    }
    ef(x,&fx);
    
    if ( gunit[1][0] != 1 ) {
    
        egradf(x,gradf);
        
    } else {
        for (j = 1 ; j <= n ; j++) {
            gradf[j] = 0.e0;
        }
        gradf[gunit[2][0]] = gunit[3][0];
    }
    gnorm[0] = max(one,vecnor(n,gradf));
    for (i = 1 ; i <= nh ; i++) {
    
        eh(i,x,&res[i]);
        
        if ( gunit[1][i] != 1 ) {
        
            egradh(i,x,yy);
            
            for (j = 1 ; j <= n ; j++) {
                gres[j][i] = yy[j];
            }
            gnorm[i] = max(one,vecnor(n,yy));
        } else {
            for (j = 1 ; j <= n ; j++) {
                gres[j][i] = 0.e0;
            }
            gnorm[i] = max(fabs((DOUBLE)gunit[3][i]),one);
            gres[gunit[2][i]][i] = gunit[3][i];
        }
    }
    for (i = 1 ; i <= ng ; i++) {
    
        eg(i,x,&res[i+nh]);
        
        if ( gunit[1][i+nh] != 1 ) {
        
            egradg(i,x,yy);
            
            for (j = 1 ; j <= n ; j++) {
                gres[j][i+nh] = yy[j];
            }
            gnorm[i+nh] = max(one,vecnor(n,yy));
        } else {
            for (j = 1 ; j <= n ; j++) {
                gres[j][i+nh] = 0.;
            }
            gres[gunit[2][i+nh]][i+nh] = gunit[3][i+nh];
            gnorm[i+nh] = max(fabs((DOUBLE)gunit[3][i+nh]),one);
        }
    }
    if ( te1 ) {
        strcpy(head,"your analytic  gradients of constraints");
        fprintf(buf,"analytic gradient of f\n");
        for (i = 1 ; i <= n ; i+= 4) {
            fprintf(buf,"    %4i   ", i);
            for (j = 0 ; j <= min(n-i,3) ; j++) {
                fprintf(buf,"  %13.6e", gradf[i+j]);
            }
            fprintf(buf,"\n");
        }
        matdr2(gres,NX,NRESM,n,nh+ng,head,buf,FALSE,0,7);
    }
    for (i = 1 ; i <= n ; i++) {
        for (j = 1 ; j <= n ; j++) {
            x0[j] = x[j];
        }
        term  = epsx*max(one,fabs(x0[i]));
        x0[i] = x0[i]+term;
         
        ef(x0,&fx0);
         
        err[i][0] = ((fx0-fx)/term-gradf[i])/gnorm[0];
        for (j = 1 ; j <= nh ; j++) {
         
            eh(j,x0,&res0[j]);
            
            err[i][j] = ((res0[j]-res[j])/term-gres[i][j])/gnorm[j];
        }
        for (j = 1 ; j <= ng ; j++) {
         
            eg(j,x0,&res0[j+nh]);
            
            err[i][j+nh]= ((res0[j+nh]-res[j+nh])/term-gres[i][j+nh])/gnorm[j+nh];
        }
    }
    fprintf(buf,"protocol of errors\n");
    for (i = 0 ; i <= nres ; i++) {
        errfnd = FALSE;
        for (j = 1 ; j <= n ; j++) {
            if ( fabs(err[j][i]) > errtol ) {
                errfnd = TRUE;
            }
        }
        if ( errfnd ) {
            fprintf(buf,"gradient nr. %i\n", i);
            for (k = 1 ; k <= n ; k+= 4) {
                fprintf(buf,"    %4i   ", k);
                for (j = 0 ; j <= min(n-k,3) ; j++) {
                    fprintf(buf,"  %13.6e", err[k+j][i]);
                }
                fprintf(buf,"\n");
            }
            fprintf(buf,"norm of analytic gradient  %10.4e\n\n", gnorm[i]);
        } else {
            fprintf(buf,"grad. nr. %4i : no errors found\n", i);
        }
    }
    countc = countc+1;
    
    goto L100;
}

/* **************************************************************************** */
/*                      euclidean norm of x , avoid overflow                    */
/* **************************************************************************** */
DOUBLE vecnor(INTEGER n,DOUBLE x[]) {
    
    static INTEGER  i;
    static DOUBLE   xm,h;
    
    if ( n <= 0 ) {

        return 0.e0;
    }
    xm = fabs(x[1]);
    for (i = 2 ; i <= n ; i++) {
        xm = max(xm,fabs(x[i]));
    }
    if ( xm == 0.e0 ) {

        return 0.e0;
        
    } else {
        h = 0.e0;
        for (i = 1 ; i <= n ; i++) {
            h = h+pow(x[i]/xm,2);
        }
        return xm*sqrt(h);
    }
}

/* **************************************************************************** */
/* subprogram for structured output of a matrix                                 */
/* uses a run time computed  format string                                      */
/* **************************************************************************** */
void matdr2(DOUBLE a[][NRESM+1],INTEGER me,INTEGER ne,INTEGER ma,INTEGER na,
            char head[],FILE *lognum,LOGICAL fix,INTEGER nd1,INTEGER nd2) {
            
    /* prints the matrix a on unit lognum using f-format if fix = TRUE          */
    /* and e-format otherwise with nd2 digits behind the point                  */
    /* if fix = FALSE nd1 is not used and treated as nd1 = 0                    */
    /* me and ne are the dimensions of the actual parameter a as defined        */
    /* in the calling program unit and ma , na the dimensions actually used     */
    /* ma<=me and na<=ne of course                                              */
    /* output is for 80-column-devices. otherwise the parameter pwid below      */
    /* may be changed                                                           */

    static INTEGER       anz,i,j,ju,jo,width,spacl,spacr;
    static const INTEGER pwid = 80;
    
    static char form1[22];

    
    if ( ma > me || na > ne || nd1 < 0 || nd2 <= 0 ) {
        fprintf(lognum,"%s\n",head);
        fprintf(lognum,"erroneous call of matdr2\n");
        
        return;
    }
    if ( ma > 999 || na > 999 ) {
        fprintf(lognum,"%s\n",head);
        fprintf(lognum,"called matdru with dim too large\n");
         
        return;
    }
    if ( fix ) {
        width = nd1+nd2+3;
    } else {
        width = nd2+7;
    }
    fprintf(lognum,"\n %s\n", head);
    
    anz   = (pwid-11)/(width+1);
    spacl = (width+1-3)/2;
    spacr = width+1-3-spacl;
    
    form1[0] = '\0';
    for (i = 1 ; i <= spacl ; i++) strcat(form1," ");
    strcat(form1,"%3i");
    for (i = 1 ; i <= spacr ; i++) strcat(form1," ");
    
    jo = 0;
    while  ( jo < na ) {
        ju = jo+1;
        jo = min(ju+anz-1,na);
        fprintf(lognum," row/column");
        for (j = ju ; j <= jo ; j++) {
            fprintf(lognum,form1,j);
            if ((j-ju+1) % anz == 0 || j == jo) fprintf(lognum,"\n");
        }
        for (i = 1 ; i <= ma ; i++) {
            fprintf(lognum,"    %4i    ", i);
            for (j = ju ; j <= jo ; j++) {
                if ( fix ) {
                fprintf(lognum,'%*.*f',width,nd2, a[i][j]);
                }
                else {
                fprintf(lognum,'%*.*e',width,nd2,a[i][j] );
                }
                if ((j-ju+1) % anz == 0 || j == jo) fprintf(lognum,"\n");
            }
        }
    }