1 #include <../src/tao/bound/impls/bqnk/bqnk.h> /*I "petsctao.h" I*/ /*I "petscmat.h" I*/ 2 #include <petscksp.h> 3 4 static PetscErrorCode TaoBQNKComputeHessian(Tao tao) 5 { 6 TAO_BNK *bnk = (TAO_BNK *)tao->data; 7 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 8 PetscErrorCode ierr; 9 PetscReal gnorm2, delta; 10 11 PetscFunctionBegin; 12 /* Alias the LMVM matrix into the TAO hessian */ 13 if (tao->hessian) { 14 ierr = MatDestroy(&tao->hessian);CHKERRQ(ierr); 15 } 16 if (tao->hessian_pre) { 17 ierr = MatDestroy(&tao->hessian_pre);CHKERRQ(ierr); 18 } 19 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 20 tao->hessian = bqnk->B; 21 ierr = PetscObjectReference((PetscObject)bqnk->B);CHKERRQ(ierr); 22 tao->hessian_pre = bqnk->B; 23 /* Update the Hessian with the latest solution */ 24 if (bqnk->is_spd) { 25 gnorm2 = bnk->gnorm*bnk->gnorm; 26 if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON; 27 if (bnk->f == 0.0) { 28 delta = 2.0 / gnorm2; 29 } else { 30 delta = 2.0 * PetscAbsScalar(bnk->f) / gnorm2; 31 } 32 ierr = MatLMVMSymBroydenSetDelta(bqnk->B, delta);CHKERRQ(ierr); 33 } 34 ierr = MatLMVMUpdate(tao->hessian, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 35 ierr = MatLMVMResetShift(tao->hessian);CHKERRQ(ierr); 36 /* Prepare the reduced sub-matrices for the inactive set */ 37 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 38 if (bnk->active_idx) { 39 ierr = MatCreateSubMatrixVirtual(tao->hessian, bnk->inactive_idx, bnk->inactive_idx, &bnk->H_inactive);CHKERRQ(ierr); 40 ierr = PCLMVMSetIS(bqnk->pc, bnk->inactive_idx);CHKERRQ(ierr); 41 } else { 42 ierr = PetscObjectReference((PetscObject)tao->hessian);CHKERRQ(ierr); 43 bnk->H_inactive = tao->hessian; 44 ierr = PCLMVMClearIS(bqnk->pc);CHKERRQ(ierr); 45 } 46 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 47 ierr = PetscObjectReference((PetscObject)bnk->H_inactive);CHKERRQ(ierr); 48 bnk->Hpre_inactive = bnk->H_inactive; 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode TaoBQNKComputeStep(Tao tao, PetscBool shift, KSPConvergedReason *ksp_reason, PetscInt *step_type) 53 { 54 TAO_BNK *bnk = (TAO_BNK *)tao->data; 55 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = TaoBNKComputeStep(tao, shift, ksp_reason, step_type);CHKERRQ(ierr); 60 if (*ksp_reason < 0) { 61 /* Krylov solver failed to converge so reset the LMVM matrix */ 62 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 63 ierr = MatLMVMUpdate(bqnk->B, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 64 } 65 PetscFunctionReturn(0); 66 } 67 68 PetscErrorCode TaoSolve_BQNK(Tao tao) 69 { 70 TAO_BNK *bnk = (TAO_BNK *)tao->data; 71 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 72 Mat_LMVM *lmvm = (Mat_LMVM*)bqnk->B->data; 73 Mat_LMVM *J0; 74 Mat_SymBrdn *diag_ctx; 75 PetscBool flg = PETSC_FALSE; 76 PetscErrorCode ierr; 77 78 PetscFunctionBegin; 79 if (!tao->recycle) { 80 ierr = MatLMVMReset(bqnk->B, PETSC_FALSE);CHKERRQ(ierr); 81 lmvm->nresets = 0; 82 if (lmvm->J0) { 83 ierr = PetscObjectBaseTypeCompare((PetscObject)lmvm->J0, MATLMVM, &flg);CHKERRQ(ierr); 84 if (flg) { 85 J0 = (Mat_LMVM*)lmvm->J0->data; 86 J0->nresets = 0; 87 } 88 } 89 flg = PETSC_FALSE; 90 ierr = PetscObjectTypeCompareAny((PetscObject)bqnk->B, &flg, MATLMVMSYMBROYDEN, MATLMVMSYMBADBROYDEN, MATLMVMBFGS, MATLMVMDFP, "");CHKERRQ(ierr); 91 if (flg) { 92 diag_ctx = (Mat_SymBrdn*)lmvm->ctx; 93 J0 = (Mat_LMVM*)diag_ctx->D->data; 94 J0->nresets = 0; 95 } 96 } 97 ierr = (*bqnk->solve)(tao);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode TaoSetUp_BQNK(Tao tao) 102 { 103 TAO_BNK *bnk = (TAO_BNK *)tao->data; 104 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 105 PetscErrorCode ierr; 106 PetscInt n, N; 107 PetscBool is_lmvm, is_sym, is_spd; 108 109 PetscFunctionBegin; 110 ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr); 111 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 112 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 113 ierr = MatSetSizes(bqnk->B, n, n, N, N);CHKERRQ(ierr); 114 ierr = MatLMVMAllocate(bqnk->B,tao->solution,bnk->unprojected_gradient);CHKERRQ(ierr); 115 ierr = PetscObjectBaseTypeCompare((PetscObject)bqnk->B, MATLMVM, &is_lmvm);CHKERRQ(ierr); 116 PetscCheck(is_lmvm,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Matrix must be an LMVM-type"); 117 ierr = MatGetOption(bqnk->B, MAT_SYMMETRIC, &is_sym);CHKERRQ(ierr); 118 PetscCheck(is_sym,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric"); 119 ierr = MatGetOption(bqnk->B, MAT_SPD, &is_spd);CHKERRQ(ierr); 120 ierr = KSPGetPC(tao->ksp, &bqnk->pc);CHKERRQ(ierr); 121 ierr = PCSetType(bqnk->pc, PCLMVM);CHKERRQ(ierr); 122 ierr = PCLMVMSetMatLMVM(bqnk->pc, bqnk->B);CHKERRQ(ierr); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode TaoSetFromOptions_BQNK(PetscOptionItems *PetscOptionsObject,Tao tao) 127 { 128 TAO_BNK *bnk = (TAO_BNK *)tao->data; 129 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = TaoSetFromOptions_BNK(PetscOptionsObject,tao);CHKERRQ(ierr); 134 if (bnk->init_type == BNK_INIT_INTERPOLATION) bnk->init_type = BNK_INIT_DIRECTION; 135 ierr = MatSetOptionsPrefix(bqnk->B, ((PetscObject)tao)->prefix);CHKERRQ(ierr); 136 ierr = MatAppendOptionsPrefix(bqnk->B, "tao_bqnk_");CHKERRQ(ierr); 137 ierr = MatSetFromOptions(bqnk->B);CHKERRQ(ierr); 138 ierr = MatGetOption(bqnk->B, MAT_SPD, &bqnk->is_spd);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode TaoView_BQNK(Tao tao, PetscViewer viewer) 143 { 144 TAO_BNK *bnk = (TAO_BNK*)tao->data; 145 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 146 PetscErrorCode ierr; 147 PetscBool isascii; 148 149 PetscFunctionBegin; 150 ierr = TaoView_BNK(tao, viewer);CHKERRQ(ierr); 151 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 152 if (isascii) { 153 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr); 154 ierr = MatView(bqnk->B, viewer);CHKERRQ(ierr); 155 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr); 156 } 157 PetscFunctionReturn(0); 158 } 159 160 static PetscErrorCode TaoDestroy_BQNK(Tao tao) 161 { 162 TAO_BNK *bnk = (TAO_BNK*)tao->data; 163 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 164 PetscErrorCode ierr; 165 166 PetscFunctionBegin; 167 ierr = MatDestroy(&bnk->Hpre_inactive);CHKERRQ(ierr); 168 ierr = MatDestroy(&bnk->H_inactive);CHKERRQ(ierr); 169 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 170 ierr = PetscFree(bnk->ctx);CHKERRQ(ierr); 171 ierr = TaoDestroy_BNK(tao);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175 PETSC_INTERN PetscErrorCode TaoCreate_BQNK(Tao tao) 176 { 177 TAO_BNK *bnk; 178 TAO_BQNK *bqnk; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 183 tao->ops->solve = TaoSolve_BQNK; 184 tao->ops->setfromoptions = TaoSetFromOptions_BQNK; 185 tao->ops->destroy = TaoDestroy_BQNK; 186 tao->ops->view = TaoView_BQNK; 187 tao->ops->setup = TaoSetUp_BQNK; 188 189 bnk = (TAO_BNK *)tao->data; 190 bnk->computehessian = TaoBQNKComputeHessian; 191 bnk->computestep = TaoBQNKComputeStep; 192 bnk->init_type = BNK_INIT_DIRECTION; 193 194 ierr = PetscNewLog(tao,&bqnk);CHKERRQ(ierr); 195 bnk->ctx = (void*)bqnk; 196 bqnk->is_spd = PETSC_TRUE; 197 198 ierr = MatCreate(PetscObjectComm((PetscObject)tao), &bqnk->B);CHKERRQ(ierr); 199 ierr = PetscObjectIncrementTabLevel((PetscObject)bqnk->B, (PetscObject)tao, 1);CHKERRQ(ierr); 200 ierr = MatSetType(bqnk->B, MATLMVMSR1);CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 /*@ 205 TaoGetLMVMMatrix - Returns a pointer to the internal LMVM matrix. Valid 206 only for quasi-Newton family of methods. 207 208 Input Parameters: 209 . tao - Tao solver context 210 211 Output Parameters: 212 . B - LMVM matrix 213 214 Level: advanced 215 216 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoSetLMVMMatrix() 217 @*/ 218 PetscErrorCode TaoGetLMVMMatrix(Tao tao, Mat *B) 219 { 220 TAO_BNK *bnk = (TAO_BNK*)tao->data; 221 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 222 PetscErrorCode ierr; 223 PetscBool flg = PETSC_FALSE; 224 225 PetscFunctionBegin; 226 ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr); 227 PetscCheck(flg,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms"); 228 *B = bqnk->B; 229 PetscFunctionReturn(0); 230 } 231 232 /*@ 233 TaoSetLMVMMatrix - Sets an external LMVM matrix into the Tao solver. Valid 234 only for quasi-Newton family of methods. 235 236 QN family of methods create their own LMVM matrices and users who wish to 237 manipulate this matrix should use TaoGetLMVMMatrix() instead. 238 239 Input Parameters: 240 + tao - Tao solver context 241 - B - LMVM matrix 242 243 Level: advanced 244 245 .seealso: TAOBQNLS, TAOBQNKLS, TAOBQNKTL, TAOBQNKTR, MATLMVM, TaoGetLMVMMatrix() 246 @*/ 247 PetscErrorCode TaoSetLMVMMatrix(Tao tao, Mat B) 248 { 249 TAO_BNK *bnk = (TAO_BNK*)tao->data; 250 TAO_BQNK *bqnk = (TAO_BQNK*)bnk->ctx; 251 PetscErrorCode ierr; 252 PetscBool flg = PETSC_FALSE; 253 254 PetscFunctionBegin; 255 ierr = PetscObjectTypeCompareAny((PetscObject)tao, &flg, TAOBQNLS, TAOBQNKLS, TAOBQNKTR, TAOBQNKTL, "");CHKERRQ(ierr); 256 PetscCheck(flg,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM Matrix only exists for quasi-Newton algorithms"); 257 ierr = PetscObjectBaseTypeCompare((PetscObject)B, MATLMVM, &flg);CHKERRQ(ierr); 258 PetscCheck(flg,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "Given matrix is not an LMVM matrix"); 259 if (bqnk->B) { 260 ierr = MatDestroy(&bqnk->B);CHKERRQ(ierr); 261 } 262 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 263 bqnk->B = B; 264 PetscFunctionReturn(0); 265 } 266