1 /* Program usage: mpiexec -n 1 maros1 [-help] [all TAO options] */ 2 3 /* ---------------------------------------------------------------------- 4 TODO Explain maros example 5 ---------------------------------------------------------------------- */ 6 7 #include <petsctao.h> 8 9 static char help[]=""; 10 11 /* 12 User-defined application context - contains data needed by the 13 application-provided call-back routines, FormFunction(), 14 FormGradient(), and FormHessian(). 15 */ 16 17 /* 18 x,d in R^n 19 f in R 20 bin in R^mi 21 beq in R^me 22 Aeq in R^(me x n) 23 Ain in R^(mi x n) 24 H in R^(n x n) 25 min f=(1/2)*x'*H*x + d'*x 26 s.t. Aeq*x == beq 27 Ain*x >= bin 28 */ 29 typedef struct { 30 char name[32]; 31 PetscInt n; /* Length x */ 32 PetscInt me; /* number of equality constraints */ 33 PetscInt mi; /* number of inequality constraints */ 34 PetscInt m; /* me+mi */ 35 Mat Aeq,Ain,H; 36 Vec beq,bin,d; 37 } AppCtx; 38 39 /* -------- User-defined Routines --------- */ 40 41 PetscErrorCode InitializeProblem(AppCtx*); 42 PetscErrorCode DestroyProblem(AppCtx *); 43 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *); 44 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*); 45 PetscErrorCode FormInequalityConstraints(Tao,Vec,Vec,void*); 46 PetscErrorCode FormEqualityConstraints(Tao,Vec,Vec,void*); 47 PetscErrorCode FormInequalityJacobian(Tao,Vec,Mat,Mat, void*); 48 PetscErrorCode FormEqualityJacobian(Tao,Vec,Mat,Mat, void*); 49 50 PetscErrorCode main(int argc,char **argv) 51 { 52 PetscMPIInt size; 53 Vec x; /* solution */ 54 KSP ksp; 55 PC pc; 56 Vec ceq,cin; 57 PetscBool flg; /* A return value when checking for use options */ 58 Tao tao; /* Tao solver context */ 59 TaoConvergedReason reason; 60 AppCtx user; /* application context */ 61 62 /* Initialize TAO,PETSc */ 63 PetscFunctionBeginUser; 64 PetscCall(PetscInitialize(&argc,&argv,(char *)0,help)); 65 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 66 /* Specify default parameters for the problem, check for command-line overrides */ 67 PetscCall(PetscStrncpy(user.name,"HS21",sizeof(user.name))); 68 PetscCall(PetscOptionsGetString(NULL,NULL,"-cutername",user.name,sizeof(user.name),&flg)); 69 70 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name)); 71 PetscCall(InitializeProblem(&user)); 72 PetscCall(VecDuplicate(user.d,&x)); 73 PetscCall(VecDuplicate(user.beq,&ceq)); 74 PetscCall(VecDuplicate(user.bin,&cin)); 75 PetscCall(VecSet(x,1.0)); 76 77 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 78 PetscCall(TaoSetType(tao,TAOIPM)); 79 PetscCall(TaoSetSolution(tao,x)); 80 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user)); 81 PetscCall(TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user)); 82 PetscCall(TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user)); 83 PetscCall(TaoSetInequalityBounds(tao,user.bin,NULL)); 84 PetscCall(TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user)); 85 PetscCall(TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user)); 86 PetscCall(TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user)); 87 PetscCall(TaoGetKSP(tao,&ksp)); 88 PetscCall(KSPGetPC(ksp,&pc)); 89 PetscCall(PCSetType(pc,PCLU)); 90 /* 91 This algorithm produces matrices with zeros along the diagonal therefore we need to use 92 SuperLU which does partial pivoting 93 */ 94 PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU)); 95 PetscCall(KSPSetType(ksp,KSPPREONLY)); 96 PetscCall(TaoSetTolerances(tao,0,0,0)); 97 98 PetscCall(TaoSetFromOptions(tao)); 99 PetscCall(TaoSolve(tao)); 100 PetscCall(TaoGetConvergedReason(tao,&reason)); 101 if (reason < 0) { 102 PetscCall(PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason])); 103 } else { 104 PetscCall(PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason])); 105 } 106 107 PetscCall(DestroyProblem(&user)); 108 PetscCall(VecDestroy(&x)); 109 PetscCall(VecDestroy(&ceq)); 110 PetscCall(VecDestroy(&cin)); 111 PetscCall(TaoDestroy(&tao)); 112 113 PetscCall(PetscFinalize()); 114 return 0; 115 } 116 117 PetscErrorCode InitializeProblem(AppCtx *user) 118 { 119 PetscViewer loader; 120 MPI_Comm comm; 121 PetscInt nrows,ncols,i; 122 PetscScalar one = 1.0; 123 char filebase[128]; 124 char filename[128]; 125 126 PetscFunctionBegin; 127 comm = PETSC_COMM_WORLD; 128 PetscCall(PetscStrncpy(filebase,user->name,sizeof(filebase))); 129 PetscCall(PetscStrlcat(filebase,"/",sizeof(filebase))); 130 PetscCall(PetscStrncpy(filename,filebase,sizeof(filename))); 131 PetscCall(PetscStrlcat(filename,"f",sizeof(filename))); 132 PetscCall(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 133 134 PetscCall(VecCreate(comm,&user->d)); 135 PetscCall(VecLoad(user->d,loader)); 136 PetscCall(PetscViewerDestroy(&loader)); 137 PetscCall(VecGetSize(user->d,&nrows)); 138 PetscCall(VecSetFromOptions(user->d)); 139 user->n = nrows; 140 141 PetscCall(PetscStrncpy(filename,filebase,sizeof(filename))); 142 PetscCall(PetscStrlcat(filename,"H",sizeof(filename))); 143 PetscCall(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 144 145 PetscCall(MatCreate(comm,&user->H)); 146 PetscCall(MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows)); 147 PetscCall(MatLoad(user->H,loader)); 148 PetscCall(PetscViewerDestroy(&loader)); 149 PetscCall(MatGetSize(user->H,&nrows,&ncols)); 150 PetscCheck(nrows == user->n,comm,PETSC_ERR_ARG_SIZ,"H: nrows != n"); 151 PetscCheck(ncols == user->n,comm,PETSC_ERR_ARG_SIZ,"H: ncols != n"); 152 PetscCall(MatSetFromOptions(user->H)); 153 154 PetscCall(PetscStrncpy(filename,filebase,sizeof(filename))); 155 PetscCall(PetscStrlcat(filename,"Aeq",sizeof(filename))); 156 PetscCall(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 157 PetscCall(MatCreate(comm,&user->Aeq)); 158 PetscCall(MatLoad(user->Aeq,loader)); 159 PetscCall(PetscViewerDestroy(&loader)); 160 PetscCall(MatGetSize(user->Aeq,&nrows,&ncols)); 161 PetscCheck(ncols == user->n,comm,PETSC_ERR_ARG_SIZ,"Aeq ncols != H nrows"); 162 PetscCall(MatSetFromOptions(user->Aeq)); 163 user->me = nrows; 164 165 PetscCall(PetscStrncpy(filename,filebase,sizeof(filename))); 166 PetscCall(PetscStrlcat(filename,"Beq",sizeof(filename))); 167 PetscCall(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 168 PetscCall(VecCreate(comm,&user->beq)); 169 PetscCall(VecLoad(user->beq,loader)); 170 PetscCall(PetscViewerDestroy(&loader)); 171 PetscCall(VecGetSize(user->beq,&nrows)); 172 PetscCheck(nrows == user->me,comm,PETSC_ERR_ARG_SIZ,"Aeq nrows != Beq n"); 173 PetscCall(VecSetFromOptions(user->beq)); 174 175 user->mi = user->n; 176 /* Ain = eye(n,n) */ 177 PetscCall(MatCreate(comm,&user->Ain)); 178 PetscCall(MatSetType(user->Ain,MATAIJ)); 179 PetscCall(MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi)); 180 181 PetscCall(MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL)); 182 PetscCall(MatSeqAIJSetPreallocation(user->Ain,1,NULL)); 183 184 for (i=0;i<user->mi;i++) PetscCall(MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES)); 185 PetscCall(MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY)); 186 PetscCall(MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY)); 187 PetscCall(MatSetFromOptions(user->Ain)); 188 189 /* bin = [0,0 ... 0]' */ 190 PetscCall(VecCreate(comm,&user->bin)); 191 PetscCall(VecSetType(user->bin,VECMPI)); 192 PetscCall(VecSetSizes(user->bin,PETSC_DECIDE,user->mi)); 193 PetscCall(VecSet(user->bin,0.0)); 194 PetscCall(VecSetFromOptions(user->bin)); 195 user->m = user->me + user->mi; 196 PetscFunctionReturn(0); 197 } 198 199 PetscErrorCode DestroyProblem(AppCtx *user) 200 { 201 PetscFunctionBegin; 202 PetscCall(MatDestroy(&user->H)); 203 PetscCall(MatDestroy(&user->Aeq)); 204 PetscCall(MatDestroy(&user->Ain)); 205 PetscCall(VecDestroy(&user->beq)); 206 PetscCall(VecDestroy(&user->bin)); 207 PetscCall(VecDestroy(&user->d)); 208 PetscFunctionReturn(0); 209 } 210 211 PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx) 212 { 213 AppCtx *user = (AppCtx*)ctx; 214 PetscScalar xtHx; 215 216 PetscFunctionBegin; 217 PetscCall(MatMult(user->H,x,g)); 218 PetscCall(VecDot(x,g,&xtHx)); 219 PetscCall(VecDot(x,user->d,f)); 220 *f += 0.5*xtHx; 221 PetscCall(VecAXPY(g,1.0,user->d)); 222 PetscFunctionReturn(0); 223 } 224 225 PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx) 226 { 227 PetscFunctionBegin; 228 PetscFunctionReturn(0); 229 } 230 231 PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx) 232 { 233 AppCtx *user = (AppCtx*)ctx; 234 235 PetscFunctionBegin; 236 PetscCall(MatMult(user->Ain,x,ci)); 237 PetscFunctionReturn(0); 238 } 239 240 PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx) 241 { 242 AppCtx *user = (AppCtx*)ctx; 243 244 PetscFunctionBegin; 245 PetscCall(MatMult(user->Aeq,x,ce)); 246 PetscCall(VecAXPY(ce,-1.0,user->beq)); 247 PetscFunctionReturn(0); 248 } 249 250 PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx) 251 { 252 PetscFunctionBegin; 253 PetscFunctionReturn(0); 254 } 255 256 PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx) 257 { 258 PetscFunctionBegin; 259 PetscFunctionReturn(0); 260 } 261 262 /*TEST 263 264 build: 265 requires: !complex 266 267 test: 268 requires: superlu 269 localrunfiles: HS21 270 271 TEST*/ 272