* Jacobi checker SUBROUTINE CHKJAC(FDF,N,M,X,DF,F,DX,W) ************************************************************************ * Check Jacobian from * SUBROUTINE FDF(N,M,X,DF,F) * REAL*8 X(N),DF(M,N),F(M) * * Input parameters * FDF : Name of the subroutine. * N : Integer. Number of independent variables. Not changed. * M : Integer. Number of functions. Not changed. * X : REAL*8 ARRAY of length N. Values of independent variables. * D : REAL*8. Steplength. Not changed. * W : REAL*8 ARRAY of length 10 + N + M*(N+3) * * Output parameters * DF : REAL*8 ARRAY of dimension M*N. Jacobian from FDF. * F : REAL*8 ARRAY of length M. Function vector from FDF. * W : REAL*8 ARRAY. Results: * W(1) : Max element in |DF| * W(2,5,6) : [Max deviation for forward approximation * row (function) number * column (x_j) number ] * W(3,7,8) : Backward approximation. * W(4,9,10): Extrapolated approximation. * * Hans Bruun Nielsen, IMM, DTU. 00.09.17 ************************************************************************ INTEGER N,M, IDAMAX, ND,NFB,NFF,NJ,NX, I,J,LW DOUBLE PRECISION X(N),DF(M,N),F(M),DX,W(*), DJ,H EXTERNAL FDF INTRINSIC ABS,DBLE,MAX,MIN * * If BLAS is available, then remove the * in col. 1 of the next line * and delete from line 112 to the end of the file * EXTERNAL DAXPY,DCOPY,DSCAL,IDAMAX C C ... Simple checks IF ((N .LT. 1) .OR. (M .LT. 1)) THEN WRITE(7,'(A)') 'N and M must be positive' RETURN ENDIF IF (DX .EQ. 0D0) THEN WRITE(7,'(A)') 'DX must be nonzero' RETURN ENDIF C ... Splitting of workspace NX = 11 NFF = NX + N NFB = NFF + M ND = NFB + M NJ = ND + M LW = NJ + M*N - 1 DO 10 I = 1,LW 10 W(I) = 0D0 C ... First call and max|DF| CALL FDF(N,M,X,DF,F) DO 20 J = 1, N 20 W(1) = MAX(W(1), ABS(DF(IDAMAX(M,DF(1,J),1),J))) C ... Difference approximations CALL DCOPY(N, X,1, W(NX),1) DO 30 J = 1,N C ... Forward difference W(NX-1+J) = X(J) + DX H = W(NX-1+J) - X(J) CALL FDF(N,M,W(NX),W(NJ),W(NFF)) CALL DAXPY(M, -1D0,F,1, W(NFF),1) CALL DSCAL(M, 1D0/H, W(NFF),1) CALL DCOPY(M, W(NFF),1, W(ND),1) CALL DAXPY(M, -1D0, DF(1,J),1, W(ND),1) I = IDAMAX(M, W(ND),1) DJ = W(ND-1+I) IF (ABS(DJ) .GT. ABS(W(2))) THEN W(2) = DJ W(5) = DBLE(I) W(6) = DBLE(J) ENDIF C ... Backward difference W(NX-1+J) = X(J) - DX/2 H = W(NX-1+J) - X(J) CALL FDF(N,M,W(NX),W(NJ),W(NFB)) CALL DAXPY(M, -1D0,F,1, W(NFB),1) CALL DSCAL(M, 1D0/H, W(NFB),1) CALL DCOPY(M, W(NFB),1, W(ND),1) CALL DAXPY(M, -1D0, DF(1,J),1, W(ND),1) I = IDAMAX(M, W(ND),1) DJ = W(ND-1+I) IF (ABS(DJ) .GT. ABS(W(3))) THEN W(3) = DJ W(7) = DBLE(I) W(8) = DBLE(J) ENDIF C ... Extrapolated difference CALL DCOPY(M, W(NFF),1, W(ND),1) CALL DAXPY(M, 2D0,W(NFB),1, W(ND),1) CALL DSCAL(M, 1D0/3D0, W(ND),1) CALL DAXPY(M, -1D0, DF(1,J),1, W(ND),1) I = IDAMAX(M, W(ND),1) DJ = W(ND-1+I) IF (ABS(DJ) .GT. ABS(W(4))) THEN W(4) = DJ W(9) = DBLE(I) W(10) = DBLE(J) ENDIF C ... Reset Jth element in x_copy W(NX-1+J) = X(J) 30 CONTINUE RETURN *-------------------------------- End of CHKJAC ------------------ END ************************************************************************ * CHKJAC calls the BLAS Level 1 subroutines and functions * DSCAL DCOPY DAXPY IDAMAX * Below we give the version obtained from * http://www.netlib.org/blas/blas.tgz ************************************************************************ subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end subroutine dcopy(n,dx,incx,dy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 i = mp1,n,7 dy(i) = dx(i) dy(i + 1) = dx(i + 1) dy(i + 2) = dx(i + 2) dy(i + 3) = dx(i + 3) dy(i + 4) = dx(i + 4) dy(i + 5) = dx(i + 5) dy(i + 6) = dx(i + 6) 50 continue return end subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dmax integer i,incx,ix,n c idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end