1 #include <petsctao.h> /*I "petsctao.h" I*/ 2 #include <petsc/private/taoimpl.h> 3 #include <petscsnes.h> 4 #include <petscdmshell.h> 5 6 /* 7 For finited difference computations of the Hessian, we use PETSc's SNESComputeJacobianDefault 8 */ 9 static PetscErrorCode Fsnes(SNES snes,Vec X,Vec G,void* ctx) 10 { 11 Tao tao = (Tao)ctx; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(tao,TAO_CLASSID,4); 15 CHKERRQ(TaoComputeGradient(tao,X,G)); 16 PetscFunctionReturn(0); 17 } 18 19 /*@C 20 TaoDefaultComputeGradient - computes the gradient using finite differences. 21 22 Collective on Tao 23 24 Input Parameters: 25 + tao - the Tao context 26 . X - compute gradient at this point 27 - dummy - not used 28 29 Output Parameter: 30 . G - Gradient Vector 31 32 Options Database Key: 33 + -tao_fd_gradient - activates TaoDefaultComputeGradient() 34 - -tao_fd_delta <delta> - change in X used to calculate finite differences 35 36 Level: advanced 37 38 Notes: 39 This routine is slow and expensive, and is not currently optimized 40 to take advantage of sparsity in the problem. Although 41 TaoDefaultComputeGradient is not recommended for general use 42 in large-scale applications, It can be useful in checking the 43 correctness of a user-provided gradient. Use the tao method TAOTEST 44 to get an indication of whether your gradient is correct. 45 This finite difference gradient evaluation can be set using the routine TaoSetGradient() or by using the command line option -tao_fd_gradient 46 47 .seealso: TaoSetGradient() 48 @*/ 49 PetscErrorCode TaoDefaultComputeGradient(Tao tao,Vec Xin,Vec G,void *dummy) 50 { 51 Vec X; 52 PetscScalar *g; 53 PetscReal f, f2; 54 PetscInt low,high,N,i; 55 PetscBool flg; 56 PetscReal h=.5*PETSC_SQRT_MACHINE_EPSILON; 57 58 PetscFunctionBegin; 59 CHKERRQ(PetscOptionsGetReal(((PetscObject)tao)->options,((PetscObject)tao)->prefix,"-tao_fd_delta",&h,&flg)); 60 CHKERRQ(VecDuplicate(Xin,&X)); 61 CHKERRQ(VecCopy(Xin,X)); 62 CHKERRQ(VecGetSize(X,&N)); 63 CHKERRQ(VecGetOwnershipRange(X,&low,&high)); 64 CHKERRQ(VecSetOption(X,VEC_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE)); 65 CHKERRQ(VecGetArray(G,&g)); 66 for (i=0;i<N;i++) { 67 CHKERRQ(VecSetValue(X,i,-h,ADD_VALUES)); 68 CHKERRQ(VecAssemblyBegin(X)); 69 CHKERRQ(VecAssemblyEnd(X)); 70 CHKERRQ(TaoComputeObjective(tao,X,&f)); 71 CHKERRQ(VecSetValue(X,i,2.0*h,ADD_VALUES)); 72 CHKERRQ(VecAssemblyBegin(X)); 73 CHKERRQ(VecAssemblyEnd(X)); 74 CHKERRQ(TaoComputeObjective(tao,X,&f2)); 75 CHKERRQ(VecSetValue(X,i,-h,ADD_VALUES)); 76 CHKERRQ(VecAssemblyBegin(X)); 77 CHKERRQ(VecAssemblyEnd(X)); 78 if (i>=low && i<high) { 79 g[i-low]=(f2-f)/(2.0*h); 80 } 81 } 82 CHKERRQ(VecRestoreArray(G,&g)); 83 CHKERRQ(VecDestroy(&X)); 84 PetscFunctionReturn(0); 85 } 86 87 /*@C 88 TaoDefaultComputeHessian - Computes the Hessian using finite differences. 89 90 Collective on Tao 91 92 Input Parameters: 93 + tao - the Tao context 94 . V - compute Hessian at this point 95 - dummy - not used 96 97 Output Parameters: 98 + H - Hessian matrix (not altered in this routine) 99 - B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 100 101 Options Database Key: 102 . -tao_fd_hessian - activates TaoDefaultComputeHessian() 103 104 Level: advanced 105 106 Notes: 107 This routine is slow and expensive, and is not currently optimized 108 to take advantage of sparsity in the problem. Although 109 TaoDefaultComputeHessian() is not recommended for general use 110 in large-scale applications, It can be useful in checking the 111 correctness of a user-provided Hessian. 112 113 .seealso: TaoSetHessian(), TaoDefaultComputeHessianColor(), SNESComputeJacobianDefault(), TaoSetGradient(), TaoDefaultComputeGradient() 114 @*/ 115 PetscErrorCode TaoDefaultComputeHessian(Tao tao,Vec V,Mat H,Mat B,void *dummy) 116 { 117 SNES snes; 118 DM dm; 119 120 PetscFunctionBegin; 121 CHKERRQ(PetscInfo(tao,"TAO Using finite differences w/o coloring to compute Hessian matrix\n")); 122 CHKERRQ(SNESCreate(PetscObjectComm((PetscObject)H),&snes)); 123 CHKERRQ(SNESSetFunction(snes,NULL,Fsnes,tao)); 124 CHKERRQ(SNESGetDM(snes,&dm)); 125 CHKERRQ(DMShellSetGlobalVector(dm,V)); 126 CHKERRQ(SNESSetUp(snes)); 127 if (H) { 128 PetscInt n,N; 129 130 CHKERRQ(VecGetSize(V,&N)); 131 CHKERRQ(VecGetLocalSize(V,&n)); 132 CHKERRQ(MatSetSizes(H,n,n,N,N)); 133 CHKERRQ(MatSetUp(H)); 134 } 135 if (B && B != H) { 136 PetscInt n,N; 137 138 CHKERRQ(VecGetSize(V,&N)); 139 CHKERRQ(VecGetLocalSize(V,&n)); 140 CHKERRQ(MatSetSizes(B,n,n,N,N)); 141 CHKERRQ(MatSetUp(B)); 142 } 143 CHKERRQ(SNESComputeJacobianDefault(snes,V,H,B,NULL)); 144 CHKERRQ(SNESDestroy(&snes)); 145 PetscFunctionReturn(0); 146 } 147 148 /*@C 149 TaoDefaultComputeHessianColor - Computes the Hessian using colored finite differences. 150 151 Collective on Tao 152 153 Input Parameters: 154 + tao - the Tao context 155 . V - compute Hessian at this point 156 - ctx - the PetscColoring object (must be of type MatFDColoring) 157 158 Output Parameters: 159 + H - Hessian matrix (not altered in this routine) 160 - B - newly computed Hessian matrix to use with preconditioner (generally the same as H) 161 162 Level: advanced 163 164 .seealso: TaoSetHessian(), TaoDefaultComputeHessian(),SNESComputeJacobianDefaultColor(), TaoSetGradient() 165 @*/ 166 PetscErrorCode TaoDefaultComputeHessianColor(Tao tao,Vec V,Mat H,Mat B,void *ctx) 167 { 168 MatFDColoring coloring = (MatFDColoring)ctx; 169 170 PetscFunctionBegin; 171 PetscValidHeaderSpecific(coloring,MAT_FDCOLORING_CLASSID,5); 172 CHKERRQ(PetscInfo(tao,"TAO computing matrix using finite differences Hessian and coloring\n")); 173 CHKERRQ(MatFDColoringApply(B,coloring,V,ctx)); 174 if (H != B) { 175 CHKERRQ(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); 176 CHKERRQ(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); 177 } 178 PetscFunctionReturn(0); 179 } 180 181 PetscErrorCode TaoDefaultComputeHessianMFFD(Tao tao,Vec X,Mat H,Mat B,void *ctx) 182 { 183 PetscInt n,N; 184 PetscBool assembled; 185 186 PetscFunctionBegin; 187 PetscCheck(!B || B == H,PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Preconditioning Hessian matrix"); 188 CHKERRQ(MatAssembled(H, &assembled)); 189 if (!assembled) { 190 CHKERRQ(VecGetSize(X,&N)); 191 CHKERRQ(VecGetLocalSize(X,&n)); 192 CHKERRQ(MatSetSizes(H,n,n,N,N)); 193 CHKERRQ(MatSetType(H,MATMFFD)); 194 CHKERRQ(MatSetUp(H)); 195 CHKERRQ(MatMFFDSetFunction(H,(PetscErrorCode (*)(void*,Vec,Vec))TaoComputeGradient,tao)); 196 } 197 CHKERRQ(MatMFFDSetBase(H,X,NULL)); 198 CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 199 CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 200 PetscFunctionReturn(0); 201 } 202