1 /* 2 This provides a thin wrapper around LMVM matrices in order to use their MatLMVMSolve 3 methods as preconditioner applications in KSP solves. 4 */ 5 6 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 7 #include <petsc/private/matimpl.h> 8 9 typedef struct { 10 Vec xwork, ywork; 11 IS inactive; 12 Mat B; 13 PetscBool allocated; 14 } PC_LMVM; 15 16 /*@ 17 PCLMVMSetMatLMVM - Replaces the LMVM matrix inside the preconditioner with 18 the one provided by the user. 19 20 Input Parameters: 21 + pc - An LMVM preconditioner 22 - B - An LMVM-type matrix (MATLDFP, MATLBFGS, MATLSR1, MATLBRDN, MATLMBRDN, MATLSBRDN) 23 24 Level: intermediate 25 @*/ 26 PetscErrorCode PCLMVMSetMatLMVM(PC pc, Mat B) 27 { 28 PC_LMVM *ctx = (PC_LMVM*)pc->data; 29 PetscErrorCode ierr; 30 PetscBool same; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 34 PetscValidHeaderSpecific(B, MAT_CLASSID, 2); 35 ierr = PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);CHKERRQ(ierr); 36 if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 37 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &same);CHKERRQ(ierr); 38 if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Matrix must be an LMVM-type."); 39 ierr = MatDestroy(&ctx->B);CHKERRQ(ierr); 40 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 41 ctx->B = B; 42 PetscFunctionReturn(0); 43 } 44 45 /*@ 46 PCLMVMGetMatLMVM - Returns a pointer to the underlying LMVM matrix. 47 48 Input Parameters: 49 . pc - An LMVM preconditioner 50 51 Output Parameters: 52 . B - LMVM matrix inside the preconditioner 53 54 Level: intermediate 55 @*/ 56 PetscErrorCode PCLMVMGetMatLMVM(PC pc, Mat *B) 57 { 58 PC_LMVM *ctx = (PC_LMVM*)pc->data; 59 PetscErrorCode ierr; 60 PetscBool same; 61 62 PetscFunctionBegin; 63 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 64 ierr = PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);CHKERRQ(ierr); 65 if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 66 *B = ctx->B; 67 PetscFunctionReturn(0); 68 } 69 70 /*@ 71 PCLMVMSetIS - Sets the index sets that reduce the PC application. 72 73 Input Parameters: 74 + pc - An LMVM preconditioner 75 - inactive - Index set defining the variables removed from the problem 76 77 Level: intermediate 78 79 .seealso: MatLMVMUpdate() 80 @*/ 81 PetscErrorCode PCLMVMSetIS(PC pc, IS inactive) 82 { 83 PC_LMVM *ctx = (PC_LMVM*)pc->data; 84 PetscErrorCode ierr; 85 PetscBool same; 86 87 PetscFunctionBegin; 88 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 89 PetscValidHeaderSpecific(inactive, IS_CLASSID, 2); 90 ierr = PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);CHKERRQ(ierr); 91 if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 92 ierr = PCLMVMClearIS(pc);CHKERRQ(ierr); 93 ierr = PetscObjectReference((PetscObject)inactive);CHKERRQ(ierr); 94 ctx->inactive = inactive; 95 PetscFunctionReturn(0); 96 } 97 98 /*@ 99 PCLMVMClearIS - Removes the inactive variable index set. 100 101 Input Parameters: 102 . pc - An LMVM preconditioner 103 104 Level: intermediate 105 106 .seealso: MatLMVMUpdate() 107 @*/ 108 PetscErrorCode PCLMVMClearIS(PC pc) 109 { 110 PC_LMVM *ctx = (PC_LMVM*)pc->data; 111 PetscErrorCode ierr; 112 PetscBool same; 113 114 PetscFunctionBegin; 115 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 116 ierr = PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &same);CHKERRQ(ierr); 117 if (!same) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC must be a PCLMVM type."); 118 if (ctx->inactive) { 119 ierr = ISDestroy(&ctx->inactive);CHKERRQ(ierr); 120 } 121 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode PCApply_LMVM(PC pc,Vec x,Vec y) 125 { 126 PC_LMVM *ctx = (PC_LMVM*)pc->data; 127 PetscErrorCode ierr; 128 Vec xsub, ysub; 129 130 PetscFunctionBegin; 131 if (ctx->inactive) { 132 ierr = VecZeroEntries(ctx->xwork);CHKERRQ(ierr); 133 ierr = VecGetSubVector(ctx->xwork, ctx->inactive, &xsub);CHKERRQ(ierr); 134 ierr = VecCopy(x, xsub);CHKERRQ(ierr); 135 ierr = VecRestoreSubVector(ctx->xwork, ctx->inactive, &xsub);CHKERRQ(ierr); 136 } else { 137 ierr = VecCopy(x, ctx->xwork);CHKERRQ(ierr); 138 } 139 ierr = MatSolve(ctx->B, ctx->xwork, ctx->ywork);CHKERRQ(ierr); 140 if (ctx->inactive) { 141 ierr = VecGetSubVector(ctx->ywork, ctx->inactive, &ysub);CHKERRQ(ierr); 142 ierr = VecCopy(ysub, y);CHKERRQ(ierr); 143 ierr = VecRestoreSubVector(ctx->ywork, ctx->inactive, &ysub);CHKERRQ(ierr); 144 } else { 145 ierr = VecCopy(ctx->ywork, y);CHKERRQ(ierr); 146 } 147 PetscFunctionReturn(0); 148 } 149 150 static PetscErrorCode PCReset_LMVM(PC pc) 151 { 152 PC_LMVM *ctx = (PC_LMVM*)pc->data; 153 PetscErrorCode ierr; 154 155 PetscFunctionBegin; 156 if (ctx->xwork) { 157 ierr = VecDestroy(&ctx->xwork);CHKERRQ(ierr); 158 } 159 if (ctx->ywork) { 160 ierr = VecDestroy(&ctx->ywork);CHKERRQ(ierr); 161 } 162 PetscFunctionReturn(0); 163 } 164 165 static PetscErrorCode PCSetUp_LMVM(PC pc) 166 { 167 PC_LMVM *ctx = (PC_LMVM*)pc->data; 168 PetscErrorCode ierr; 169 PetscInt n, N; 170 PetscBool allocated; 171 172 PetscFunctionBegin; 173 ierr = MatLMVMIsAllocated(ctx->B, &allocated);CHKERRQ(ierr); 174 if (!allocated) { 175 ierr = MatCreateVecs(pc->mat, &ctx->xwork, &ctx->ywork);CHKERRQ(ierr); 176 ierr = VecGetLocalSize(ctx->xwork, &n);CHKERRQ(ierr); 177 ierr = VecGetSize(ctx->xwork, &N);CHKERRQ(ierr); 178 ierr = MatSetSizes(ctx->B, n, n, N, N);CHKERRQ(ierr); 179 ierr = MatLMVMAllocate(ctx->B, ctx->xwork, ctx->ywork);CHKERRQ(ierr); 180 } else { 181 ierr = MatCreateVecs(ctx->B, &ctx->xwork, &ctx->ywork);CHKERRQ(ierr); 182 } 183 PetscFunctionReturn(0); 184 } 185 186 static PetscErrorCode PCView_LMVM(PC pc,PetscViewer viewer) 187 { 188 PC_LMVM *ctx = (PC_LMVM*)pc->data; 189 PetscErrorCode ierr; 190 PetscBool iascii; 191 192 PetscFunctionBegin; 193 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 194 if (iascii && ctx->B->assembled) { 195 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 196 ierr = MatView(ctx->B, viewer);CHKERRQ(ierr); 197 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 198 } 199 PetscFunctionReturn(0); 200 } 201 202 static PetscErrorCode PCSetFromOptions_LMVM(PetscOptionItems* PetscOptionsObject, PC pc) 203 { 204 PC_LMVM *ctx = (PC_LMVM*)pc->data; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 ierr = MatSetFromOptions(ctx->B);CHKERRQ(ierr); 209 PetscFunctionReturn(0); 210 } 211 212 static PetscErrorCode PCDestroy_LMVM(PC pc) 213 { 214 PC_LMVM *ctx = (PC_LMVM*)pc->data; 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 if (ctx->inactive) { 219 ierr = ISDestroy(&ctx->inactive);CHKERRQ(ierr); 220 } 221 if (pc->setupcalled) { 222 ierr = VecDestroy(&ctx->xwork);CHKERRQ(ierr); 223 ierr = VecDestroy(&ctx->ywork);CHKERRQ(ierr); 224 } 225 ierr = MatDestroy(&ctx->B);CHKERRQ(ierr); 226 ierr = PetscFree(pc->data);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 /*MC 231 PCLMVM - Creates a preconditioner around an LMVM matrix. Options for the 232 underlying LMVM matrix can be access with the "-pc_lmvm_" prefix. 233 234 Level: intermediate 235 236 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), 237 PC, MATLMVM, PCLMVMUpdate(), PCLMVMSetMatLMVM(), PCLMVMGetMatLMVM() 238 M*/ 239 PETSC_EXTERN PetscErrorCode PCCreate_LMVM(PC pc) 240 { 241 PetscErrorCode ierr; 242 PC_LMVM *ctx; 243 244 PetscFunctionBegin; 245 ierr = PetscNewLog(pc,&ctx);CHKERRQ(ierr); 246 pc->data = (void*)ctx; 247 248 pc->ops->reset = PCReset_LMVM; 249 pc->ops->setup = PCSetUp_LMVM; 250 pc->ops->destroy = PCDestroy_LMVM; 251 pc->ops->view = PCView_LMVM; 252 pc->ops->apply = PCApply_LMVM; 253 pc->ops->setfromoptions = PCSetFromOptions_LMVM; 254 pc->ops->applysymmetricleft = NULL; 255 pc->ops->applysymmetricright = NULL; 256 pc->ops->applytranspose = NULL; 257 pc->ops->applyrichardson = NULL; 258 pc->ops->presolve = NULL; 259 pc->ops->postsolve = NULL; 260 261 ierr = PCSetReusePreconditioner(pc, PETSC_TRUE);CHKERRQ(ierr); 262 263 ierr = MatCreate(PetscObjectComm((PetscObject)pc), &ctx->B);CHKERRQ(ierr); 264 ierr = MatSetType(ctx->B, MATLMVMBFGS);CHKERRQ(ierr); 265 ierr = PetscObjectIncrementTabLevel((PetscObject)ctx->B, (PetscObject)pc, 1);CHKERRQ(ierr); 266 ierr = MatSetOptionsPrefix(ctx->B, "pc_lmvm_");CHKERRQ(ierr); 267 PetscFunctionReturn(0); 268 } 269 270 271 272 273 274 275