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 int 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 223 } 224 225 PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx) 226 { 227 PetscFunctionBegin; 228 PetscFunctionReturn(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 248 } 249 250 PetscErrorCode FormInequalityJacobian(Tao tao, Vec x, Mat JI, Mat JIpre, void *ctx) 251 { 252 PetscFunctionBegin; 253 PetscFunctionReturn(PETSC_SUCCESS); 254 } 255 256 PetscErrorCode FormEqualityJacobian(Tao tao, Vec x, Mat JE, Mat JEpre, void *ctx) 257 { 258 PetscFunctionBegin; 259 PetscFunctionReturn(PETSC_SUCCESS); 260 } 261 262 /*TEST 263 264 build: 265 requires: !complex 266 267 test: 268 requires: !single superlu 269 localrunfiles: HS21 270 271 TEST*/ 272