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