Actual source code: maros.c

      1: /* Program usage: mpiexec -n 1 maros1 [-help] [all TAO options] */

      3: /* ----------------------------------------------------------------------
      4: TODO Explain maros example
      5: ---------------------------------------------------------------------- */

      7: #include <petsctao.h>

      9: static  char help[]="";

     11: /*T
     12:    Concepts: TAO^Solving an unconstrained minimization problem
     13:    Routines: TaoCreate(); TaoSetType();
     14:    Routines: TaoSetInitialVector();
     15:    Routines: TaoSetObjectiveAndGradientRoutine();
     16:    Routines: TaoSetEqualityConstraintsRoutine();
     17:    Routines: TaoSetInequalityConstraintsRoutine();
     18:    Routines: TaoSetEqualityJacobianRoutine();
     19:    Routines: TaoSetInequalityJacobianRoutine();
     20:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
     21:    Routines: TaoGetKSP(); TaoSolve();
     22:    Routines: TaoGetConvergedReason();TaoDestroy();
     23:    Processors: 1
     24: T*/

     26: /*
     27:    User-defined application context - contains data needed by the
     28:    application-provided call-back routines, FormFunction(),
     29:    FormGradient(), and FormHessian().
     30: */

     32: /*
     33:    x,d in R^n
     34:    f in R
     35:    bin in R^mi
     36:    beq in R^me
     37:    Aeq in R^(me x n)
     38:    Ain in R^(mi x n)
     39:    H in R^(n x n)
     40:    min f=(1/2)*x'*H*x + d'*x
     41:    s.t.  Aeq*x == beq
     42:          Ain*x >= bin
     43: */
     44: typedef struct {
     45:   char     name[32];
     46:   PetscInt n; /* Length x */
     47:   PetscInt me; /* number of equality constraints */
     48:   PetscInt mi; /* number of inequality constraints */
     49:   PetscInt m;  /* me+mi */
     50:   Mat      Aeq,Ain,H;
     51:   Vec      beq,bin,d;
     52: } AppCtx;

     54: /* -------- User-defined Routines --------- */

     56: PetscErrorCode InitializeProblem(AppCtx*);
     57: PetscErrorCode DestroyProblem(AppCtx *);
     58: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
     59: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
     60: PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*);
     61: PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*);
     62: PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*);
     63: PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*);

     65: PetscErrorCode main(int argc,char **argv)
     66: {
     67:   PetscErrorCode     ierr;                /* used to check for functions returning nonzeros */
     68:   PetscMPIInt        size;
     69:   Vec                x;                   /* solution */
     70:   KSP                ksp;
     71:   PC                 pc;
     72:   Vec                ceq,cin;
     73:   PetscBool          flg;                 /* A return value when checking for use options */
     74:   Tao                tao;                 /* Tao solver context */
     75:   TaoConvergedReason reason;
     76:   AppCtx             user;                /* application context */

     78:   /* Initialize TAO,PETSc */
     79:   PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;
     80:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
     81:   /* Specify default parameters for the problem, check for command-line overrides */
     82:   PetscStrncpy(user.name,"HS21",sizeof(user.name));
     83:   PetscOptionsGetString(NULL,NULL,"-cutername",user.name,sizeof(user.name),&flg);

     85:   PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name);
     86:   InitializeProblem(&user);
     87:   VecDuplicate(user.d,&x);
     88:   VecDuplicate(user.beq,&ceq);
     89:   VecDuplicate(user.bin,&cin);
     90:   VecSet(x,1.0);

     92:   TaoCreate(PETSC_COMM_WORLD,&tao);
     93:   TaoSetType(tao,TAOIPM);
     94:   TaoSetInitialVector(tao,x);
     95:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);
     96:   TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user);
     97:   TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user);
     98:   TaoSetInequalityBounds(tao,user.bin,NULL);
     99:   TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user);
    100:   TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user);
    101:   TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
    102:   TaoGetKSP(tao,&ksp);
    103:   KSPGetPC(ksp,&pc);
    104:   PCSetType(pc,PCLU);
    105:   /*
    106:       This algorithm produces matrices with zeros along the diagonal therefore we need to use
    107:     SuperLU which does partial pivoting
    108:   */
    109:   PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);
    110:   KSPSetType(ksp,KSPPREONLY);
    111:   TaoSetTolerances(tao,0,0,0);

    113:   TaoSetFromOptions(tao);
    114:   TaoSolve(tao);
    115:   TaoGetConvergedReason(tao,&reason);
    116:   if (reason < 0) {
    117:     PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason]);
    118:   } else {
    119:     PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason]);
    120:   }

    122:   DestroyProblem(&user);
    123:   VecDestroy(&x);
    124:   VecDestroy(&ceq);
    125:   VecDestroy(&cin);
    126:   TaoDestroy(&tao);

    128:   PetscFinalize();
    129:   return ierr;
    130: }

    132: PetscErrorCode InitializeProblem(AppCtx *user)
    133: {
    135:   PetscViewer    loader;
    136:   MPI_Comm       comm;
    137:   PetscInt       nrows,ncols,i;
    138:   PetscScalar    one=1.0;
    139:   char           filebase[128];
    140:   char           filename[128];

    143:   comm = PETSC_COMM_WORLD;
    144:   PetscStrncpy(filebase,user->name,sizeof(filebase));
    145:   PetscStrlcat(filebase,"/",sizeof(filebase));
    146:   PetscStrncpy(filename,filebase,sizeof(filename));
    147:   PetscStrlcat(filename,"f",sizeof(filename));
    148:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);

    150:   VecCreate(comm,&user->d);
    151:   VecLoad(user->d,loader);
    152:   PetscViewerDestroy(&loader);
    153:   VecGetSize(user->d,&nrows);
    154:   VecSetFromOptions(user->d);
    155:   user->n = nrows;

    157:   PetscStrncpy(filename,filebase,sizeof(filename));
    158:   PetscStrlcat(filename,"H",sizeof(filename));
    159:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);

    161:   MatCreate(comm,&user->H);
    162:   MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);
    163:   MatLoad(user->H,loader);
    164:   PetscViewerDestroy(&loader);
    165:   MatGetSize(user->H,&nrows,&ncols);
    166:   if (nrows != user->n) SETERRQ(comm,0,"H: nrows != n\n");
    167:   if (ncols != user->n) SETERRQ(comm,0,"H: ncols != n\n");
    168:   MatSetFromOptions(user->H);

    170:   PetscStrncpy(filename,filebase,sizeof(filename));
    171:   PetscStrlcat(filename,"Aeq",sizeof(filename));
    172:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
    173:   if (ierr) {
    174:     user->Aeq = NULL;
    175:     user->me  = 0;
    176:   } else {
    177:     MatCreate(comm,&user->Aeq);
    178:     MatLoad(user->Aeq,loader);
    179:     PetscViewerDestroy(&loader);
    180:     MatGetSize(user->Aeq,&nrows,&ncols);
    181:     if (ncols != user->n) SETERRQ(comm,0,"Aeq ncols != H nrows\n");
    182:     MatSetFromOptions(user->Aeq);
    183:     user->me = nrows;
    184:   }

    186:   PetscStrncpy(filename,filebase,sizeof(filename));
    187:   PetscStrlcat(filename,"Beq",sizeof(filename));
    188:   PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);
    189:   if (ierr) {
    190:     user->beq = 0;
    191:   } else {
    192:     VecCreate(comm,&user->beq);
    193:     VecLoad(user->beq,loader);
    194:     PetscViewerDestroy(&loader);
    195:     VecGetSize(user->beq,&nrows);
    196:     if (nrows != user->me) SETERRQ(comm,0,"Aeq nrows != Beq n\n");
    197:     VecSetFromOptions(user->beq);
    198:   }

    200:   user->mi = user->n;
    201:   /* Ain = eye(n,n) */
    202:   MatCreate(comm,&user->Ain);
    203:   MatSetType(user->Ain,MATAIJ);
    204:   MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);

    206:   MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);
    207:   MatSeqAIJSetPreallocation(user->Ain,1,NULL);

    209:   for (i=0;i<user->mi;i++) {
    210:     MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);
    211:   }
    212:   MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);
    213:   MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);
    214:   MatSetFromOptions(user->Ain);

    216:   /* bin = [0,0 ... 0]' */
    217:   VecCreate(comm,&user->bin);
    218:   VecSetType(user->bin,VECMPI);
    219:   VecSetSizes(user->bin,PETSC_DECIDE,user->mi);
    220:   VecSet(user->bin,0.0);
    221:   VecSetFromOptions(user->bin);
    222:   user->m = user->me + user->mi;
    223:   return(0);
    224: }

    226: PetscErrorCode DestroyProblem(AppCtx *user)
    227: {

    231:   MatDestroy(&user->H);
    232:   MatDestroy(&user->Aeq);
    233:   MatDestroy(&user->Ain);
    234:   VecDestroy(&user->beq);
    235:   VecDestroy(&user->bin);
    236:   VecDestroy(&user->d);
    237:   return(0);
    238: }
    239: PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
    240: {
    241:   AppCtx         *user = (AppCtx*)ctx;
    242:   PetscScalar    xtHx;

    246:   MatMult(user->H,x,g);
    247:   VecDot(x,g,&xtHx);
    248:   VecDot(x,user->d,f);
    249:   *f += 0.5*xtHx;
    250:   VecAXPY(g,1.0,user->d);
    251:   return(0);
    252: }

    254: PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
    255: {
    257:   return(0);
    258: }

    260: PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx)
    261: {
    262:   AppCtx         *user = (AppCtx*)ctx;

    266:   MatMult(user->Ain,x,ci);
    267:   return(0);
    268: }

    270: PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx)
    271: {
    272:   AppCtx         *user = (AppCtx*)ctx;

    276:   MatMult(user->Aeq,x,ce);
    277:   VecAXPY(ce,-1.0,user->beq);
    278:   return(0);
    279: }

    281: PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre,  void *ctx)
    282: {
    284:   return(0);
    285: }

    287: PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx)
    288: {
    290:   return(0);
    291: }

    294: /*TEST

    296:    build:
    297:       requires: !complex

    299:    test:
    300:       requires: superlu
    301:       localrunfiles: HS21

    303: TEST*/