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