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 /*T 12 Concepts: TAO^Solving an unconstrained minimization problem 13 Routines: TaoCreate(); TaoSetType(); 14 Routines: TaoSetSolution(); 15 Routines: TaoSetObjectiveAndGradient(); 16 Routines: TaoSetEqualityConstraintsRoutine(); 17 Routines: TaoSetInequalityConstraintsRoutine(); 18 Routines: TaoSetEqualityJacobianRoutine(); 19 Routines: TaoSetInequalityJacobianRoutine(); 20 Routines: TaoSetHessian(); TaoSetFromOptions(); 21 Routines: TaoGetKSP(); TaoSolve(); 22 Routines: TaoGetConvergedReason();TaoDestroy(); 23 Processors: 1 24 T*/ 25 26 /* 27 User-defined application context - contains data needed by the 28 application-provided call-back routines, FormFunction(), 29 FormGradient(), and FormHessian(). 30 */ 31 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; 53 54 /* -------- User-defined Routines --------- */ 55 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*); 64 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 */ 77 78 /* Initialize TAO,PETSc */ 79 ierr = PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr; 80 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 81 /* Specify default parameters for the problem, check for command-line overrides */ 82 CHKERRQ(PetscStrncpy(user.name,"HS21",sizeof(user.name))); 83 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-cutername",user.name,sizeof(user.name),&flg)); 84 85 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name)); 86 CHKERRQ(InitializeProblem(&user)); 87 CHKERRQ(VecDuplicate(user.d,&x)); 88 CHKERRQ(VecDuplicate(user.beq,&ceq)); 89 CHKERRQ(VecDuplicate(user.bin,&cin)); 90 CHKERRQ(VecSet(x,1.0)); 91 92 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 93 CHKERRQ(TaoSetType(tao,TAOIPM)); 94 CHKERRQ(TaoSetSolution(tao,x)); 95 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*)&user)); 96 CHKERRQ(TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user)); 97 CHKERRQ(TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user)); 98 CHKERRQ(TaoSetInequalityBounds(tao,user.bin,NULL)); 99 CHKERRQ(TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user)); 100 CHKERRQ(TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user)); 101 CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void*)&user)); 102 CHKERRQ(TaoGetKSP(tao,&ksp)); 103 CHKERRQ(KSPGetPC(ksp,&pc)); 104 CHKERRQ(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 CHKERRQ(PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU)); 110 CHKERRQ(KSPSetType(ksp,KSPPREONLY)); 111 CHKERRQ(TaoSetTolerances(tao,0,0,0)); 112 113 CHKERRQ(TaoSetFromOptions(tao)); 114 CHKERRQ(TaoSolve(tao)); 115 CHKERRQ(TaoGetConvergedReason(tao,&reason)); 116 if (reason < 0) { 117 CHKERRQ(PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason])); 118 } else { 119 CHKERRQ(PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason])); 120 } 121 122 CHKERRQ(DestroyProblem(&user)); 123 CHKERRQ(VecDestroy(&x)); 124 CHKERRQ(VecDestroy(&ceq)); 125 CHKERRQ(VecDestroy(&cin)); 126 CHKERRQ(TaoDestroy(&tao)); 127 128 ierr = PetscFinalize(); 129 return ierr; 130 } 131 132 PetscErrorCode InitializeProblem(AppCtx *user) 133 { 134 PetscViewer loader; 135 MPI_Comm comm; 136 PetscInt nrows,ncols,i; 137 PetscScalar one = 1.0; 138 char filebase[128]; 139 char filename[128]; 140 141 PetscFunctionBegin; 142 comm = PETSC_COMM_WORLD; 143 CHKERRQ(PetscStrncpy(filebase,user->name,sizeof(filebase))); 144 CHKERRQ(PetscStrlcat(filebase,"/",sizeof(filebase))); 145 CHKERRQ(PetscStrncpy(filename,filebase,sizeof(filename))); 146 CHKERRQ(PetscStrlcat(filename,"f",sizeof(filename))); 147 CHKERRQ(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 148 149 CHKERRQ(VecCreate(comm,&user->d)); 150 CHKERRQ(VecLoad(user->d,loader)); 151 CHKERRQ(PetscViewerDestroy(&loader)); 152 CHKERRQ(VecGetSize(user->d,&nrows)); 153 CHKERRQ(VecSetFromOptions(user->d)); 154 user->n = nrows; 155 156 CHKERRQ(PetscStrncpy(filename,filebase,sizeof(filename))); 157 CHKERRQ(PetscStrlcat(filename,"H",sizeof(filename))); 158 CHKERRQ(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 159 160 CHKERRQ(MatCreate(comm,&user->H)); 161 CHKERRQ(MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows)); 162 CHKERRQ(MatLoad(user->H,loader)); 163 CHKERRQ(PetscViewerDestroy(&loader)); 164 CHKERRQ(MatGetSize(user->H,&nrows,&ncols)); 165 PetscCheck(nrows == user->n,comm,PETSC_ERR_ARG_SIZ,"H: nrows != n"); 166 PetscCheck(ncols == user->n,comm,PETSC_ERR_ARG_SIZ,"H: ncols != n"); 167 CHKERRQ(MatSetFromOptions(user->H)); 168 169 CHKERRQ(PetscStrncpy(filename,filebase,sizeof(filename))); 170 CHKERRQ(PetscStrlcat(filename,"Aeq",sizeof(filename))); 171 CHKERRQ(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 172 CHKERRQ(MatCreate(comm,&user->Aeq)); 173 CHKERRQ(MatLoad(user->Aeq,loader)); 174 CHKERRQ(PetscViewerDestroy(&loader)); 175 CHKERRQ(MatGetSize(user->Aeq,&nrows,&ncols)); 176 PetscCheck(ncols == user->n,comm,PETSC_ERR_ARG_SIZ,"Aeq ncols != H nrows"); 177 CHKERRQ(MatSetFromOptions(user->Aeq)); 178 user->me = nrows; 179 180 CHKERRQ(PetscStrncpy(filename,filebase,sizeof(filename))); 181 CHKERRQ(PetscStrlcat(filename,"Beq",sizeof(filename))); 182 CHKERRQ(PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader)); 183 CHKERRQ(VecCreate(comm,&user->beq)); 184 CHKERRQ(VecLoad(user->beq,loader)); 185 CHKERRQ(PetscViewerDestroy(&loader)); 186 CHKERRQ(VecGetSize(user->beq,&nrows)); 187 PetscCheck(nrows == user->me,comm,PETSC_ERR_ARG_SIZ,"Aeq nrows != Beq n"); 188 CHKERRQ(VecSetFromOptions(user->beq)); 189 190 user->mi = user->n; 191 /* Ain = eye(n,n) */ 192 CHKERRQ(MatCreate(comm,&user->Ain)); 193 CHKERRQ(MatSetType(user->Ain,MATAIJ)); 194 CHKERRQ(MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi)); 195 196 CHKERRQ(MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL)); 197 CHKERRQ(MatSeqAIJSetPreallocation(user->Ain,1,NULL)); 198 199 for (i=0;i<user->mi;i++) CHKERRQ(MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES)); 200 CHKERRQ(MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY)); 201 CHKERRQ(MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY)); 202 CHKERRQ(MatSetFromOptions(user->Ain)); 203 204 /* bin = [0,0 ... 0]' */ 205 CHKERRQ(VecCreate(comm,&user->bin)); 206 CHKERRQ(VecSetType(user->bin,VECMPI)); 207 CHKERRQ(VecSetSizes(user->bin,PETSC_DECIDE,user->mi)); 208 CHKERRQ(VecSet(user->bin,0.0)); 209 CHKERRQ(VecSetFromOptions(user->bin)); 210 user->m = user->me + user->mi; 211 PetscFunctionReturn(0); 212 } 213 214 PetscErrorCode DestroyProblem(AppCtx *user) 215 { 216 PetscFunctionBegin; 217 CHKERRQ(MatDestroy(&user->H)); 218 CHKERRQ(MatDestroy(&user->Aeq)); 219 CHKERRQ(MatDestroy(&user->Ain)); 220 CHKERRQ(VecDestroy(&user->beq)); 221 CHKERRQ(VecDestroy(&user->bin)); 222 CHKERRQ(VecDestroy(&user->d)); 223 PetscFunctionReturn(0); 224 } 225 226 PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx) 227 { 228 AppCtx *user = (AppCtx*)ctx; 229 PetscScalar xtHx; 230 231 PetscFunctionBegin; 232 CHKERRQ(MatMult(user->H,x,g)); 233 CHKERRQ(VecDot(x,g,&xtHx)); 234 CHKERRQ(VecDot(x,user->d,f)); 235 *f += 0.5*xtHx; 236 CHKERRQ(VecAXPY(g,1.0,user->d)); 237 PetscFunctionReturn(0); 238 } 239 240 PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx) 241 { 242 PetscFunctionBegin; 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx) 247 { 248 AppCtx *user = (AppCtx*)ctx; 249 250 PetscFunctionBegin; 251 CHKERRQ(MatMult(user->Ain,x,ci)); 252 PetscFunctionReturn(0); 253 } 254 255 PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx) 256 { 257 AppCtx *user = (AppCtx*)ctx; 258 259 PetscFunctionBegin; 260 CHKERRQ(MatMult(user->Aeq,x,ce)); 261 CHKERRQ(VecAXPY(ce,-1.0,user->beq)); 262 PetscFunctionReturn(0); 263 } 264 265 PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx) 266 { 267 PetscFunctionBegin; 268 PetscFunctionReturn(0); 269 } 270 271 PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx) 272 { 273 PetscFunctionBegin; 274 PetscFunctionReturn(0); 275 } 276 277 /*TEST 278 279 build: 280 requires: !complex 281 282 test: 283 requires: superlu 284 localrunfiles: HS21 285 286 TEST*/ 287