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