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: 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*/ 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 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 81 /* Specify default parameters for the problem, check for command-line overrides */ 82 ierr = PetscStrncpy(user.name,"HS21",8);CHKERRQ(ierr); 83 ierr = PetscOptionsGetString(NULL,NULL,"-cutername",user.name,24,&flg);CHKERRQ(ierr); 84 85 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n---- MAROS Problem %s -----\n",user.name);CHKERRQ(ierr); 86 ierr = InitializeProblem(&user);CHKERRQ(ierr); 87 ierr = VecDuplicate(user.d,&x);CHKERRQ(ierr); 88 ierr = VecDuplicate(user.beq,&ceq);CHKERRQ(ierr); 89 ierr = VecDuplicate(user.bin,&cin);CHKERRQ(ierr); 90 ierr = VecSet(x,1.0);CHKERRQ(ierr); 91 92 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 93 ierr = TaoSetType(tao,TAOIPM);CHKERRQ(ierr); 94 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr); 95 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr); 96 ierr = TaoSetEqualityConstraintsRoutine(tao,ceq,FormEqualityConstraints,(void*)&user);CHKERRQ(ierr); 97 ierr = TaoSetInequalityConstraintsRoutine(tao,cin,FormInequalityConstraints,(void*)&user);CHKERRQ(ierr); 98 ierr = TaoSetInequalityBounds(tao,user.bin,NULL);CHKERRQ(ierr); 99 ierr = TaoSetJacobianEqualityRoutine(tao,user.Aeq,user.Aeq,FormEqualityJacobian,(void*)&user);CHKERRQ(ierr); 100 ierr = TaoSetJacobianInequalityRoutine(tao,user.Ain,user.Ain,FormInequalityJacobian,(void*)&user);CHKERRQ(ierr); 101 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr); 102 ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr); 103 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 104 ierr = PCSetType(pc,PCLU);CHKERRQ(ierr); 105 /* 106 This algorithm produces matrices with zeros along the diagonal therefore we need to use 107 SuperLU which does partial pivoting 108 */ 109 ierr = PCFactorSetMatSolverType(pc,MATSOLVERSUPERLU);CHKERRQ(ierr); 110 ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr); 111 ierr = TaoSetTolerances(tao,0,0,0);CHKERRQ(ierr); 112 113 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 114 ierr = TaoSolve(tao);CHKERRQ(ierr); 115 ierr = TaoGetConvergedReason(tao,&reason);CHKERRQ(ierr); 116 if (reason < 0) { 117 ierr = PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge due to %s.\n",TaoConvergedReasons[reason]);CHKERRQ(ierr); 118 } else { 119 ierr = PetscPrintf(MPI_COMM_WORLD, "Optimization completed with status %s.\n",TaoConvergedReasons[reason]);CHKERRQ(ierr); 120 } 121 122 ierr = DestroyProblem(&user);CHKERRQ(ierr); 123 ierr = VecDestroy(&x);CHKERRQ(ierr); 124 ierr = VecDestroy(&ceq);CHKERRQ(ierr); 125 ierr = VecDestroy(&cin);CHKERRQ(ierr); 126 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 127 128 ierr = PetscFinalize(); 129 return ierr; 130 } 131 132 PetscErrorCode InitializeProblem(AppCtx *user) 133 { 134 PetscErrorCode ierr; 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]; 141 142 PetscFunctionBegin; 143 comm = PETSC_COMM_WORLD; 144 ierr = PetscStrncpy(filebase,user->name,sizeof(filebase));CHKERRQ(ierr); 145 ierr = PetscStrlcat(filebase,"/",sizeof(filebase));CHKERRQ(ierr); 146 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr); 147 ierr = PetscStrlcat(filename,"f",sizeof(filename));CHKERRQ(ierr); 148 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);CHKERRQ(ierr); 149 150 ierr = VecCreate(comm,&user->d);CHKERRQ(ierr); 151 ierr = VecLoad(user->d,loader);CHKERRQ(ierr); 152 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr); 153 ierr = VecGetSize(user->d,&nrows);CHKERRQ(ierr); 154 ierr = VecSetFromOptions(user->d);CHKERRQ(ierr); 155 user->n = nrows; 156 157 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr); 158 ierr = PetscStrlcat(filename,"H",sizeof(filename));CHKERRQ(ierr); 159 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);CHKERRQ(ierr); 160 161 ierr = MatCreate(comm,&user->H);CHKERRQ(ierr); 162 ierr = MatSetSizes(user->H,PETSC_DECIDE,PETSC_DECIDE,nrows,nrows);CHKERRQ(ierr); 163 ierr = MatLoad(user->H,loader);CHKERRQ(ierr); 164 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr); 165 ierr = MatGetSize(user->H,&nrows,&ncols);CHKERRQ(ierr); 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 ierr = MatSetFromOptions(user->H);CHKERRQ(ierr); 169 170 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr); 171 ierr = PetscStrlcat(filename,"Aeq",sizeof(filename));CHKERRQ(ierr); 172 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader); 173 if (ierr) { 174 user->Aeq = NULL; 175 user->me = 0; 176 } else { 177 ierr = MatCreate(comm,&user->Aeq);CHKERRQ(ierr); 178 ierr = MatLoad(user->Aeq,loader);CHKERRQ(ierr); 179 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr); 180 ierr = MatGetSize(user->Aeq,&nrows,&ncols);CHKERRQ(ierr); 181 if (ncols != user->n) SETERRQ(comm,0,"Aeq ncols != H nrows\n"); 182 ierr = MatSetFromOptions(user->Aeq);CHKERRQ(ierr); 183 user->me = nrows; 184 } 185 186 ierr = PetscStrncpy(filename,filebase,sizeof(filename));CHKERRQ(ierr); 187 ierr = PetscStrlcat(filename,"Beq",sizeof(filename));CHKERRQ(ierr); 188 ierr = PetscViewerBinaryOpen(comm,filename,FILE_MODE_READ,&loader);CHKERRQ(ierr); 189 if (ierr) { 190 user->beq = 0; 191 } else { 192 ierr = VecCreate(comm,&user->beq);CHKERRQ(ierr); 193 ierr = VecLoad(user->beq,loader);CHKERRQ(ierr); 194 ierr = PetscViewerDestroy(&loader);CHKERRQ(ierr); 195 ierr = VecGetSize(user->beq,&nrows);CHKERRQ(ierr); 196 if (nrows != user->me) SETERRQ(comm,0,"Aeq nrows != Beq n\n"); 197 ierr = VecSetFromOptions(user->beq);CHKERRQ(ierr); 198 } 199 200 user->mi = user->n; 201 /* Ain = eye(n,n) */ 202 ierr = MatCreate(comm,&user->Ain);CHKERRQ(ierr); 203 ierr = MatSetType(user->Ain,MATAIJ);CHKERRQ(ierr); 204 ierr = MatSetSizes(user->Ain,PETSC_DECIDE,PETSC_DECIDE,user->mi,user->mi);CHKERRQ(ierr); 205 206 ierr = MatMPIAIJSetPreallocation(user->Ain,1,NULL,0,NULL);CHKERRQ(ierr); 207 ierr = MatSeqAIJSetPreallocation(user->Ain,1,NULL);CHKERRQ(ierr); 208 209 for (i=0;i<user->mi;i++) { 210 ierr = MatSetValues(user->Ain,1,&i,1,&i,&one,INSERT_VALUES);CHKERRQ(ierr); 211 } 212 ierr = MatAssemblyBegin(user->Ain,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 213 ierr = MatAssemblyEnd(user->Ain,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 214 ierr = MatSetFromOptions(user->Ain);CHKERRQ(ierr); 215 216 /* bin = [0,0 ... 0]' */ 217 ierr = VecCreate(comm,&user->bin);CHKERRQ(ierr); 218 ierr = VecSetType(user->bin,VECMPI);CHKERRQ(ierr); 219 ierr = VecSetSizes(user->bin,PETSC_DECIDE,user->mi);CHKERRQ(ierr); 220 ierr = VecSet(user->bin,0.0);CHKERRQ(ierr); 221 ierr = VecSetFromOptions(user->bin);CHKERRQ(ierr); 222 user->m = user->me + user->mi; 223 PetscFunctionReturn(0); 224 } 225 226 PetscErrorCode DestroyProblem(AppCtx *user) 227 { 228 PetscErrorCode ierr; 229 230 PetscFunctionBegin; 231 ierr = MatDestroy(&user->H);CHKERRQ(ierr); 232 ierr = MatDestroy(&user->Aeq);CHKERRQ(ierr); 233 ierr = MatDestroy(&user->Ain);CHKERRQ(ierr); 234 ierr = VecDestroy(&user->beq);CHKERRQ(ierr); 235 ierr = VecDestroy(&user->bin);CHKERRQ(ierr); 236 ierr = VecDestroy(&user->d);CHKERRQ(ierr); 237 PetscFunctionReturn(0); 238 } 239 PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx) 240 { 241 AppCtx *user = (AppCtx*)ctx; 242 PetscScalar xtHx; 243 PetscErrorCode ierr; 244 245 PetscFunctionBegin; 246 ierr = MatMult(user->H,x,g);CHKERRQ(ierr); 247 ierr = VecDot(x,g,&xtHx);CHKERRQ(ierr); 248 ierr = VecDot(x,user->d,f);CHKERRQ(ierr); 249 *f += 0.5*xtHx; 250 ierr = VecAXPY(g,1.0,user->d);CHKERRQ(ierr); 251 PetscFunctionReturn(0); 252 } 253 254 PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx) 255 { 256 PetscFunctionBegin; 257 PetscFunctionReturn(0); 258 } 259 260 PetscErrorCode FormInequalityConstraints(Tao tao, Vec x, Vec ci, void *ctx) 261 { 262 AppCtx *user = (AppCtx*)ctx; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 ierr = MatMult(user->Ain,x,ci);CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 PetscErrorCode FormEqualityConstraints(Tao tao, Vec x, Vec ce,void *ctx) 271 { 272 AppCtx *user = (AppCtx*)ctx; 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 ierr = MatMult(user->Aeq,x,ce);CHKERRQ(ierr); 277 ierr = VecAXPY(ce,-1.0,user->beq);CHKERRQ(ierr); 278 PetscFunctionReturn(0); 279 } 280 281 PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx) 282 { 283 PetscFunctionBegin; 284 PetscFunctionReturn(0); 285 } 286 287 PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx) 288 { 289 PetscFunctionBegin; 290 PetscFunctionReturn(0); 291 } 292 293 294 /*TEST 295 296 build: 297 requires: !complex 298 299 test: 300 requires: superlu 301 localrunfiles: HS21 302 303 TEST*/ 304