subroutine da_check_gradient(grid, config_flags, cv_size, xhat, cv, pdx, itertest, xbx, be, iv, y, re, j_cost) !----------------------------------------------------------------------- ! Purpose: Gradient test ! Adopted from grtest.f90 of GSI by Xin Zhang , September, 2011 !----------------------------------------------------------------------- implicit none type(domain), intent(inout) :: grid type(grid_config_rec_type), intent(inout) :: config_flags real, intent(in) :: pdx integer, intent(in ) :: itertest integer, intent(in) :: cv_size ! Size of cv array. real, intent(inout) :: xhat(1:cv_size) ! control variable (local). real, intent(inout) :: cv(1:cv_size) ! control variable (local). type (xbx_type), intent(inout) :: xbx ! Header & non-gridded vars. type (be_type), intent(in) :: be ! background error structure. type (iv_type), intent(inout) :: iv ! ob. increment vector. type (y_type), intent(inout) :: y ! y = h (xa) type (y_type), intent(inout) :: re ! residual (o-a) structure. type (j_type), intent(inout) :: j_cost ! cost function real :: xdir(1:cv_size) , grad(1:cv_size), yhat(1:cv_size) real :: zfy,zf0,zdf0,za,zfa,zdfa real :: zabuf(itertest),zfabuf(itertest),ztf2buf(itertest) real :: ZT1,ZB,ZFB,ztco,ZTC1,ZT2,ZTC1A,ZTC2,ZTF2 real :: ZAL,ZFAL,ZBL,ZFBL,ZTF2L real :: ZTC00,ZTC02,ZTC10,ZTC12 real :: ZERMIN,ZT1TST,ZREF integer :: jp_start, jp_end ! Start/end indices of Jp. #ifdef CLOUD_CV integer :: mz(11) integer :: sz(9) #else integer :: mz(7) integer :: sz(5) #endif integer :: ibest,idig integer :: ii, jj call da_trace_entry("da_check_gradient") #ifdef CLOUD_CV mz = (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%v6%mz, be%v7%mz,be%v8%mz,be%v9%mz, be%alpha%mz, be % ne /) sz = (/ be%cv%size1, be%cv%size2, be%cv%size3, be%cv%size4, be%cv%size5, be%cv%size6, be%cv%size7,be%cv%size8,be%cv%size9/) #else mz = (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz, be % ne /) sz = (/ be%cv%size1, be%cv%size2, be%cv%size3, be%cv%size4, be%cv%size5 /) #endif jp_start = be % cv % size_jb + be % cv % size_je + 1 jp_end = be % cv % size_jb + be % cv % size_je + be % cv % size_jp call da_message((/' gradient test starting'/)) if (pdx<=EPSILON(pdx)) then if (rootproc) write(6,*)'grtest, pdx=',pdx write(6,*)'grtest: pdx too small',pdx return endif !---------------------------------------------------------------------------- ! [1] Initial point !---------------------------------------------------------------------------- call da_initialize_cv (cv_size, cv) call da_initialize_cv (cv_size, xhat) ! Initialize cv values with random data: call random_number(xhat(:)) xhat(:) = xhat(:) - 0.5 if (rootproc) write(6,*)'grtest: use random_number(xhat)' yhat=xhat call da_calculate_j(1, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & be%cv%size_jl, xbx, be, iv, yhat, cv, re, y, j_cost, grad, grid, config_flags) zfy = j_cost%total !---------------------------------------------------------------------------- ! [1.1] Define perturbation direction ZH !---------------------------------------------------------------------------- call da_message((/' The test direction is the opposite of the gradient '/)) xdir = -1.0 * grad !---------------------------------------------------------------------------- ! [1.2] Set function f value and derivative at origin !---------------------------------------------------------------------------- zf0=zfy zdf0=da_dot_cv(cv_size,grad,xdir,grid,mz,jp_start,jp_end) if (rootproc) write(6,*)'grtest: F(0)=',zf0,' DF(0)=',zdf0 IF (ZDF0>0.0) write(6,*) 'GRTEST Warning, DF should be negative' IF (ABS(ZDF0) < SQRT(EPSILON(ZDF0))) THEN if (rootproc) write(6,*) 'GRTEST WARNING, DERIVATIVE IS TOO SMALL' ENDIF !---------------------------------------------------------------------------- ! [2] Loop on test point !---------------------------------------------------------------------------- ztf2buf(1)=0.0 DO jj=1,itertest za=pdx*(10.0**(jj-1)) if (rootproc) write(6,*)'grtest iter=',jj,' alpha=',za !---------------------------------------------------------------------------- ! [2.1] Compute f and df at new point y=x+a.h !---------------------------------------------------------------------------- do ii=1,cv_size yhat(ii) = xhat(ii) + za * xdir(ii) end do call da_calculate_j(1, 1, cv_size, be%cv%size_jb, be%cv%size_je, be%cv%size_jp, & be%cv%size_jl, xbx, be, iv, yhat, cv, re, y, j_cost, grad, grid, config_flags) zfy = j_cost%total zfa=zfy zdfa=da_dot_cv(cv_size,grad,xdir,grid,mz,jp_start,jp_end) if (rootproc) write(6,*)'grtest: alpha=',za,' F(a)=',zfa,' DF(a)=',zdfa zabuf(jj)=za zfabuf(jj)=zfa !---------------------------------------------------------------------------- ! [2.2] Quantity TC0=f(a)/f(0)-1 !---------------------------------------------------------------------------- ! if f is continuous then TC0->1 at origin, ! at least linearly with a. IF (ABS(zf0)<=TINY(zf0)) THEN ! do not compute T1 in this unlikely case if (rootproc) write(6,*) 'grtest: Warning: zf0 is suspiciously small.' if (rootproc) write(6,*) 'grtest: F(a)-F(0)=',zfa-zf0 ELSE ztco=zfa/zf0-1.0 if (rootproc) write(6,*)'grtest: continuity TC0=',ztco ENDIF !---------------------------------------------------------------------------- ! f(a)-f(0) ! [2.3] Quantity T1=----------- ! a.df(0) !---------------------------------------------------------------------------- ! if df is the gradient then T1->1 at origin, ! linearly with a. T1 is undefined if df(0)=0. IF (ABS(za*zdf0)<=SQRT(TINY(zf0))) THEN if (rootproc) write(6,*)'grtest: Warning: could not compute ',& & 'gradient test T1, a.df(0)=',za*zdf0 ELSE zt1=(zfa-zf0)/(za*zdf0) if (rootproc) write(6,*)'grtest: gradient T1=',zt1 ENDIF !---------------------------------------------------------------------------- ! [2.4] Quantity TC1=( f(a)-f(0)-a.df(0) )/a !---------------------------------------------------------------------------- ! if df is the gradient and df is continuous, ! then TC1->0 linearly with a. ZTC1=(ZFA-ZF0-ZA*ZDF0)/ZA if (rootproc) write(6,*)'grtest: grad continuity TC1=',ZTC1 !---------------------------------------------------------------------------- ! [2.5] Quantity T2=( f(a)-f(0)-a.df(0) )*2/a**2 !---------------------------------------------------------------------------- ! if d2f exists then T2 -> d2f(0) linearly with a. ZT2=(ZFA-ZF0-ZA*ZDF0)*2.0/(ZA**2) if (rootproc) write(6,*)'grtest: second derivative T2=',ZT2 !---------------------------------------------------------------------------- ! [2.6] Quantity TC1A=df(a)-df(0) !---------------------------------------------------------------------------- ! if df is the gradient in a and df is continuous, ! then TC1A->0 linearly with a. ZTC1A=ZDFA-ZDF0 if (rootproc) write(6,*)'grtest: a-grad continuity TC1A=',ZTC1A !---------------------------------------------------------------------------- ! [2.7] Quantity TC2=( 2(f(0)-f(a))+ a(df(0)+df(a))/a**2 !---------------------------------------------------------------------------- ! if f is exactly quadratic, then TC2=0, always: numerically ! it has to -> 0 when a is BIG. Otherwise TC2->0 linearly for ! small a is trivially implied by TC1A and T2. ZTC2=(2.0*(ZF0-ZFA)+ZA*(ZDF0+ZDFA))/(ZA**2) if (rootproc) write(6,*)'grtest: quadraticity TC2=',ZTC2 !---------------------------------------------------------------------------- ! 2 f(0)-f(b) f(a)-f(b) ! [2.8] Quantity TF2=---( --------- + --------- ) ! a b a-b !---------------------------------------------------------------------------- ! if 0, a and b are distinct and f is quadratic then ! TF2=d2f, always. The estimate is most stable when a,b are big. ! This test works only after two loops, but it is immune against ! gradient bugs. IF (jj>=2) THEN ZB =ZABUF (jj-1) ZFB=ZFABUF(jj-1) ZTF2=2.0/ZA*((ZF0-ZFB)/ZB+(ZFA-ZFB)/(ZA-ZB)) if (rootproc) write(6,*)'grtest: convexity ZTF2=',ZTF2 ztf2buf(jj)=ztf2 ENDIF ! End loop ENDDO !---------------------------------------------------------------------------- ! [3] Comment on the results !---------------------------------------------------------------------------- ! TC0(0)/TC0(2)<.011 -> df looks continuous ! item with (T1<1 and 1-T1 is min) = best grad test item ! reldif(TF2(last),TF2(last-1)) = precision on quadraticity !---------------------------------------------------------------------------- ! 3.1 Fundamental checks !---------------------------------------------------------------------------- if (rootproc) then write(6,*) 'GRTEST: TENTATIVE CONCLUSIONS :' ZTC00=ABS(zfabuf(1)-zf0) ZTC02=ABS(zfabuf(3)-zf0) IF( ZTC00/zabuf(1) <= 1.5*(ZTC02/zabuf(3)) )THEN write(6,*) 'GRTEST: function f looks continous.' ELSE write(6,*) 'GRTEST: WARNING f does not look continuous',& & ' (perhaps truncation problem)' ENDIF !---------------------------------------------------------------------------- ! 3.2 Gradient quality !---------------------------------------------------------------------------- IF (ABS(zdf0)<=SQRT(TINY(zf0))) THEN write(6,*) 'GRTEST: The gradient is 0, which is unusual !' ZTC10=ABS(zfabuf(1)-zf0) ZTC12=ABS(zfabuf(3)-zf0) IF( ZTC10/zabuf(1)**2 <= 1.1*ZTC12/zabuf(3)**2)THEN write(6,*)'GRTEST: The gradient looks good anyway.' ENDIF ELSE ! Find best gradient test index ZERMIN=HUGE(0.0) ibest=-1 DO jj=1,itertest ZT1TST=(zfabuf(jj)-zf0)/(zabuf(jj)*zdf0) ZT1TST=ABS(ZT1TST-1.0) IF (ZT1TST