program su4 implicit none integer i1,i2,ia,ix,i,nmax,ik complex kappa(4,4),delta(4,4),h(100,100),u(4,4),v(4,4) complex vf(100,100) real d(100),e(100),t1,rho(4,4),kin real spin,isospin,mu_spin nmax=8 c initial guess for kappa matrix: call init(kappa) print *,' kappa :' do i1=1,4 write(6,12)(real(kappa(i1,i2)),i2=1,4) enddo read*,mu_spin do ik=1,3 c determine delta for the given kappa call v_cf(kappa,delta) print*,' delta :' do i2=1,4 write(6,12)(real(delta(i2,i1)),i1=1,4) end do c generate the HfB matrix h: h=0.d0 h(1,1)= 0.05 + mu_spin/4.0 h(2,2)= 0.05 - mu_spin/4.0 h(3,3)= 0.05 + mu_spin/4.0 h(4,4) = 0.05 - mu_spin/4.0 do i1=1,4 h(i1+4,i1+4)=-h(i1,i1) enddo h(1:4,5:8)= delta(:,:) h(5:8,1:4) = delta(:,:) print*,' hfb matrix' do i2=1,8 write(6,12)(real(h(i2,i1)),i1=1,8) write(16,12)(real(h(i2,i1)),i1=1,8) 12 format(9f8.4) enddo c eigenvalues (d) and eigenvector matrix (vf) of HfB: call eig_c(nmax,h,d,vf) ia=0 c use eigenvectors corresponding to positive eigenvalues to generate u and v matrix: do ix=1,8 if(d(ix).gt.0) then print *,d(ix) ia=ia+1 e(ia)=d(ix) u(1:4,ia)=vf(1:4,ix) v(1:4,ia) = vf(5:8,ix) endif enddo if(ia.ne.4) stop 'trouble 1' c generate the new kappa matrix : kappa=0.0 do ia=1,4 do i1=1,4 kappa(i1,1:4) = kappa(i1,1:4) + u(i1,ia)*conjg(v(1:4,ia)) enddo enddo print*, 'Re( kappa )' do i1=1,4 write(6,12) (real(kappa(i1,i2)),i2=1,4) enddo print*, 'Im( kappa )' do i1=1,4 write(6,12) (aimag(kappa(i1,i2)),i2=1,4) enddo c compute new delta matrix: call v_cf(kappa,delta) c compute energy: t1 = 0.0 do i1=1,4 do i2=1,4 t1=t1+0.5*conjg(kappa(i1,i2))*delta(i1,i2) enddo enddo print*,'pair energy = ' ,t1 print*, ' rho density matrix' do i1=1,4 do i2=1,4 rho(i1,i2)=0.0 do ia=1,4 rho(i1,i2)= rho(i1,i2) + v(i1,ia)*conjg(v(i2,ia)) enddo enddo write(6,12)(rho(i1,i2),i2=1,4) enddo kin = 0.0 do i1=1,4 kin = kin + rho(i1,i1)*h(i1,i1) end do spin = (rho(1,1)+rho(3,3)) - (rho(2,2)+rho(4,4)) isospin = (rho(1,1)+rho(2,2)) - (rho(3,3)+rho(4,4)) c pause enddo print*, 'mu_spin, mu_ispin,e_corr,v_pair,kin,rho_spin,rho_isospin' write(6,12)real(h(1,1)-h(2,2)),real(h(1,1)-h(3,3)), & t1+kin,t1,kin,spin,isospin c write(23,12)real(h(1,1)-h(2,2)),real(h(1,1)-h(3,3)), c & t1+kin,t1,kin,spin,isospin end program su4 c----------------------------------------------------------------------- subroutine init(kappa) implicit none complex kappa(4,4) kappa=0.0 kappa(1,2) = 0.0 kappa(2,1) = 0.0 kappa(1,3) = -1.0 kappa(3,1) = 1.0 kappa(2,4) = 1.0 kappa(4,2) = -1.0 kappa(3,4) = -0.0 kappa(4,3) = 0.0 end subroutine init c------------------------------------------------------------------------ subroutine v_cf(kappa,delta) implicit none complex delta(4,4),kappa(4,4),s1(4,4),s2(4,4),pot(16,16) complex s3(4,4),kvec(16) real v1,v2 integer i,j,i1,i2,j1,j2,if1,ic1,if2,ic2,i3,i4,i3p,i4p do i=1,16 do j=1,16 i1= 1 + (i-1 - mod((i-1),4))/4 i2 = 1 + mod(i-1,4) j1 = 1 + (j-1 - mod((j-1),4))/4 j2 = 1 + mod(j-1,4) v1 = 0.0 v2 = 0.0 if(i1.eq.j1.and.i2.eq.j2)v1 = -1.0 if(i1.eq.j2.and.i2.eq.j1)v2 = 1.0 pot(i,j) = v1 + v2 end do end do c do i2=1,16 c write(6,12)(real(pot(i2,i1)),i1=1,8) c write(6,12)(real(pot(i2,i1)),i1=9,16) c c end do c 12 format(9f8.2) c pause do i1=1,4 do i2=1,4 kvec(4*(i1-1)+i2)=kappa(i1,i2) end do enddo delta = 0.0 do i1=1,4 do i2=1,4 do i3=1,16 delta(i1,i2) = delta(i1,i2) & + 0.5*pot(4*(i1-1)+i2,i3)*kvec(i3) end do end do end do end subroutine v_cf c------------------------------------------------------- c subroutine eig_c(n,a,e,vf) real*8 ar(100,100),ai(100,100),w(100),zr(100,100),zi(100,100), x fv1(100),fv2(100),fm1(2,100) complex a(100,100),vf(100,100),ei real e(100) ei=(0.0,1.0) matz=1 do i1=1,n do i2=1,n ar(i1,i2)=real(a(i1,i2)) ai(i1,i2)=aimag(a(i1,i2)) enddo enddo call ch(100,n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) c print*,'eigenvalues ' c write(6,12) (w(i),i=1,n) e(1:n) = w(1:n) c print*, ' eigenvectors ' do i1= 1,n c write(6,14) (zr(i1,ia),ia=1,n) do ia = 1,n vf(i1,ia) = zr(i1,ia) + ei*zi (i1,ia) enddo enddo c print*, ' ' do i1=1,n c write(6,12) (real(vf(i1,ia)),ia=1,n) 14 format(4(3x,2f8.4)) enddo c print*,' ' 12 format(8f10.5) if(ierr.ne.0) print*,'error code ',ierr end subroutine ch(nm,n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) double precision ar(nm,n),ai(nm,n),w(n),zr(nm,n),zi(nm,n), x fv1(n),fv2(n),fm1(2,n) c integer i,j,n,nm,ierr,matz c c this subroutine calls the recommended sequence of c subroutines from the eigensystem subroutine package (eispack) c to find the eigenvalues and eigenvectors (if desired) c of a complex hermitian matrix. c c on input c c nm must be set to the row dimension of the two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix a=(ar,ai). c c ar and ai contain the real and imaginary parts, c respectively, of the complex hermitian matrix. c c matz is an integer variable set equal to zero if c only eigenvalues are desired. otherwise it is set to c any non-zero integer for both eigenvalues and eigenvectors. c c on output c c w contains the eigenvalues in ascending order. c c zr and zi contain the real and imaginary parts, c respectively, of the eigenvectors if matz is not zero. c c ierr is an integer output variable set equal to an error c completion code described in the documentation for tqlrat c and tql2. the normal completion code is zero. c c fv1, fv2, and fm1 are temporary storage arrays. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (n .le. nm) go to 10 ierr = 10 * n go to 50 c 10 call htridi(nm,n,ar,ai,w,fv1,fv2,fm1) if (matz .ne. 0) go to 20 c .......... find eigenvalues only .......... call tqlrat(n,w,fv2,ierr) go to 50 c .......... find both eigenvalues and eigenvectors .......... 20 do 40 i = 1, n c do 30 j = 1, n zr(j,i) = 0.0d0 30 continue c zr(i,i) = 1.0d0 40 continue c call tql2(nm,n,w,fv1,zr,ierr) if (ierr .ne. 0) go to 50 call htribk(nm,n,ar,ai,fm1,n,zr,zi) 50 return end subroutine tql2(nm,n,d,e,z,ierr) c integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr real*8 d(n),e(n),z(nm,n) real*8 c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag c c this subroutine is a translation of the algol procedure tql2, c num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and c wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). c c this subroutine finds the eigenvalues and eigenvectors c of a symmetric tridiagonal matrix by the ql method. c the eigenvectors of a full symmetric matrix can also c be found if tred2 has been used to reduce this c full matrix to tridiagonal form. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e contains the subdiagonal elements of the input matrix c in its last n-1 positions. e(1) is arbitrary. c c z contains the transformation matrix produced in the c reduction by tred2, if performed. if the eigenvectors c of the tridiagonal matrix are desired, z must contain c the identity matrix. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct but c unordered for indices 1,2,...,ierr-1. c c e has been destroyed. c c z contains orthonormal eigenvectors of the symmetric c tridiagonal (or full) matrix. if an error exit is made, c z contains the eigenvectors associated with the stored c eigenvalues. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e(i-1) = e(i) c f = 0.0d0 tst1 = 0.0d0 e(n) = 0.0d0 c do 240 l = 1, n j = 0 h = dabs(d(l)) + dabs(e(l)) if (tst1 .lt. h) tst1 = h c .......... look for small sub-diagonal element .......... do 110 m = l, n tst2 = tst1 + dabs(e(m)) if (tst2 .eq. tst1) go to 120 c .......... e(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 220 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 l2 = l1 + 1 g = d(l) p = (d(l1) - g) / (2.0d0 * e(l)) r = pythag(p,1.0d0) d(l) = e(l) / (p + dsign(r,p)) d(l1) = e(l) * (p + dsign(r,p)) dl1 = d(l1) h = g - d(l) if (l2 .gt. n) go to 145 c do 140 i = l2, n 140 d(i) = d(i) - h c 145 f = f + h c .......... ql transformation .......... p = d(m) c = 1.0d0 c2 = c el1 = e(l1) s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml c3 = c2 c2 = c s2 = s i = m - ii g = c * e(i) h = c * p r = pythag(p,e(i)) e(i+1) = s * r s = e(i) / r c = p / r p = c * d(i) - s * g d(i+1) = h + s * (c * g + s * d(i)) c .......... form vector .......... do 180 k = 1, n h = z(k,i+1) z(k,i+1) = s * z(k,i) + c * h z(k,i) = c * z(k,i) - s * h 180 continue c 200 continue c p = -s * s2 * c3 * el1 * e(l) / dl1 e(l) = s * p d(l) = c * p tst2 = tst1 + dabs(e(l)) if (tst2 .gt. tst1) go to 130 220 d(l) = d(l) + f 240 continue c .......... order eigenvalues and eigenvectors .......... do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p c do 280 j = 1, n p = z(j,i) z(j,i) = z(j,k) z(j,k) = p 280 continue c 300 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end real*8 function pythag(a,b) real*8 a,b c c finds dsqrt(a**2+b**2) without overflow or destructive underflow c real*8 p,r,s,t,u p = dmax1(dabs(a),dabs(b)) if (p .eq. 0.0d0) go to 20 r = (dmin1(dabs(a),dabs(b))/p)**2 10 continue t = 4.0d0 + r if (t .eq. 4.0d0) go to 20 s = r/t u = 1.0d0 + 2.0d0*s p = u*p r = (s/u)**2 * r go to 10 20 pythag = p return end subroutine tqlrat(n,d,e2,ierr) c integer i,j,l,m,n,ii,l1,mml,ierr real*8 d(n),e2(n) real*8 b,c,f,g,h,p,r,s,t,epslon,pythag c c this subroutine is a translation of the algol procedure tqlrat, c algorithm 464, comm. acm 16, 689(1973) by reinsch. c c this subroutine finds the eigenvalues of a symmetric c tridiagonal matrix by the rational ql method. c c on input c c n is the order of the matrix. c c d contains the diagonal elements of the input matrix. c c e2 contains the squares of the subdiagonal elements of the c input matrix in its last n-1 positions. e2(1) is arbitrary. c c on output c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct and c ordered for indices 1,2,...ierr-1, but may not be c the smallest eigenvalues. c c e2 has been destroyed. c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c ierr = 0 if (n .eq. 1) go to 1001 c do 100 i = 2, n 100 e2(i-1) = e2(i) c f = 0.0d0 t = 0.0d0 e2(n) = 0.0d0 c do 290 l = 1, n j = 0 h = dabs(d(l)) + dsqrt(e2(l)) if (t .gt. h) go to 105 t = h b = epslon(t) c = b * b c .......... look for small squared sub-diagonal element .......... 105 do 110 m = l, n if (e2(m) .le. c) go to 120 c .......... e2(n) is always zero, so there is no exit c through the bottom of the loop .......... 110 continue c 120 if (m .eq. l) go to 210 130 if (j .eq. 30) go to 1000 j = j + 1 c .......... form shift .......... l1 = l + 1 s = dsqrt(e2(l)) g = d(l) p = (d(l1) - g) / (2.0d0 * s) r = pythag(p,1.0d0) d(l) = s / (p + dsign(r,p)) h = g - d(l) c do 140 i = l1, n 140 d(i) = d(i) - h c f = f + h c .......... rational ql transformation .......... g = d(m) if (g .eq. 0.0d0) g = b h = g s = 0.0d0 mml = m - l c .......... for i=m-1 step -1 until l do -- .......... do 200 ii = 1, mml i = m - ii p = g * h r = p + e2(i) e2(i+1) = s * r s = e2(i) / r d(i+1) = h + s * (h + d(i)) g = d(i) - e2(i) / g if (g .eq. 0.0d0) g = b h = g * p / r 200 continue c e2(l) = s * g d(l) = h c .......... guard against underflow in convergence test .......... if (h .eq. 0.0d0) go to 210 if (dabs(e2(l)) .le. dabs(c/h)) go to 210 e2(l) = h * e2(l) if (e2(l) .ne. 0.0d0) go to 130 210 p = d(l) + f c .......... order eigenvalues .......... if (l .eq. 1) go to 250 c .......... for i=l step -1 until 2 do -- .......... do 230 ii = 2, l i = l + 2 - ii if (p .ge. d(i-1)) go to 270 d(i) = d(i-1) 230 continue c 250 i = 1 270 d(i) = p 290 continue c go to 1001 c .......... set error -- no convergence to an c eigenvalue after 30 iterations .......... 1000 ierr = l 1001 return end subroutine htridi(nm,n,ar,ai,d,e,e2,tau) c integer i,j,k,l,n,ii,nm,jp1 real*8 ar(nm,n),ai(nm,n),d(n),e(n),e2(n),tau(2,n) real*8 f,g,h,fi,gi,hh,si,scale,pythag c c this subroutine is a translation of a complex analogue of c the algol procedure tred1, num. math. 11, 181-195(1968) c by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine reduces a complex hermitian matrix c to a real symmetric tridiagonal matrix using c unitary similarity transformations. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c ar and ai contain the real and imaginary parts, c respectively, of the complex hermitian input matrix. c only the lower triangle of the matrix need be supplied. c c on output c c ar and ai contain information about the unitary trans- c formations used in the reduction in their full lower c triangles. their strict upper triangles and the c diagonal of ar are unaltered. c c d contains the diagonal elements of the the tridiagonal matrix. c c e contains the subdiagonal elements of the tridiagonal c matrix in its last n-1 positions. e(1) is set to zero. c c e2 contains the squares of the corresponding elements of e. c e2 may coincide with e if the squares are not needed. c c tau contains further information about the transformations. c c calls pythag for dsqrt(a*a + b*b) . c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c tau(1,n) = 1.0d0 tau(2,n) = 0.0d0 c do 100 i = 1, n 100 d(i) = ar(i,i) c .......... for i=n step -1 until 1 do -- .......... do 300 ii = 1, n i = n + 1 - ii l = i - 1 h = 0.0d0 scale = 0.0d0 if (l .lt. 1) go to 130 c .......... scale row (algol tol then not needed) .......... do 120 k = 1, l 120 scale = scale + dabs(ar(i,k)) + dabs(ai(i,k)) c if (scale .ne. 0.0d0) go to 140 tau(1,l) = 1.0d0 tau(2,l) = 0.0d0 130 e(i) = 0.0d0 e2(i) = 0.0d0 go to 290 c 140 do 150 k = 1, l ar(i,k) = ar(i,k) / scale ai(i,k) = ai(i,k) / scale h = h + ar(i,k) * ar(i,k) + ai(i,k) * ai(i,k) 150 continue c e2(i) = scale * scale * h g = dsqrt(h) e(i) = scale * g f = pythag(ar(i,l),ai(i,l)) c .......... form next diagonal element of matrix t .......... if (f .eq. 0.0d0) go to 160 tau(1,l) = (ai(i,l) * tau(2,i) - ar(i,l) * tau(1,i)) / f si = (ar(i,l) * tau(2,i) + ai(i,l) * tau(1,i)) / f h = h + f * g g = 1.0d0 + g / f ar(i,l) = g * ar(i,l) ai(i,l) = g * ai(i,l) if (l .eq. 1) go to 270 go to 170 160 tau(1,l) = -tau(1,i) si = tau(2,i) ar(i,l) = g 170 f = 0.0d0 c do 240 j = 1, l g = 0.0d0 gi = 0.0d0 c .......... form element of a*u .......... do 180 k = 1, j g = g + ar(j,k) * ar(i,k) + ai(j,k) * ai(i,k) gi = gi - ar(j,k) * ai(i,k) + ai(j,k) * ar(i,k) 180 continue c jp1 = j + 1 if (l .lt. jp1) go to 220 c do 200 k = jp1, l g = g + ar(k,j) * ar(i,k) - ai(k,j) * ai(i,k) gi = gi - ar(k,j) * ai(i,k) - ai(k,j) * ar(i,k) 200 continue c .......... form element of p .......... 220 e(j) = g / h tau(2,j) = gi / h f = f + e(j) * ar(i,j) - tau(2,j) * ai(i,j) 240 continue c hh = f / (h + h) c .......... form reduced a .......... do 260 j = 1, l f = ar(i,j) g = e(j) - hh * f e(j) = g fi = -ai(i,j) gi = tau(2,j) - hh * fi tau(2,j) = -gi c do 260 k = 1, j ar(j,k) = ar(j,k) - f * e(k) - g * ar(i,k) x + fi * tau(2,k) + gi * ai(i,k) ai(j,k) = ai(j,k) - f * tau(2,k) - g * ai(i,k) x - fi * e(k) - gi * ar(i,k) 260 continue c 270 do 280 k = 1, l ar(i,k) = scale * ar(i,k) ai(i,k) = scale * ai(i,k) 280 continue c tau(2,l) = -si 290 hh = d(i) d(i) = ar(i,i) ar(i,i) = hh ai(i,i) = scale * dsqrt(h) 300 continue c return end subroutine htribk(nm,n,ar,ai,tau,m,zr,zi) c integer i,j,k,l,m,n,nm real*8 ar(nm,n),ai(nm,n),tau(2,n),zr(nm,m),zi(nm,m) real*8 h,s,si c c this subroutine is a translation of a complex analogue of c the algol procedure trbak1, num. math. 11, 181-195(1968) c by martin, reinsch, and wilkinson. c handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). c c this subroutine forms the eigenvectors of a complex hermitian c matrix by back transforming those of the corresponding c real symmetric tridiagonal matrix determined by htridi. c c on input c c nm must be set to the row dimension of two-dimensional c array parameters as declared in the calling program c dimension statement. c c n is the order of the matrix. c c ar and ai contain information about the unitary trans- c formations used in the reduction by htridi in their c full lower triangles except for the diagonal of ar. c c tau contains further information about the transformations. c c m is the number of eigenvectors to be back transformed. c c zr contains the eigenvectors to be back transformed c in its first m columns. c c on output c c zr and zi contain the real and imaginary parts, c respectively, of the transformed eigenvectors c in their first m columns. c c note that the last component of each returned vector c is real and that vector euclidean norms are preserved. c c questions and comments should be directed to burton s. garbow, c mathematics and computer science div, argonne national laboratory c c this version dated august 1983. c c ------------------------------------------------------------------ c if (m .eq. 0) go to 200 c .......... transform the eigenvectors of the real symmetric c tridiagonal matrix to those of the hermitian c tridiagonal matrix. .......... do 50 k = 1, n c do 50 j = 1, m zi(k,j) = -zr(k,j) * tau(2,k) zr(k,j) = zr(k,j) * tau(1,k) 50 continue c if (n .eq. 1) go to 200 c .......... recover and apply the householder matrices .......... do 140 i = 2, n l = i - 1 h = ai(i,i) if (h .eq. 0.0d0) go to 140 c do 130 j = 1, m s = 0.0d0 si = 0.0d0 c do 110 k = 1, l s = s + ar(i,k) * zr(k,j) - ai(i,k) * zi(k,j) si = si + ar(i,k) * zi(k,j) + ai(i,k) * zr(k,j) 110 continue c .......... double divisions avoid possible underflow .......... s = (s / h) / h si = (si / h) / h c do 120 k = 1, l zr(k,j) = zr(k,j) - s * ar(i,k) - si * ai(i,k) zi(k,j) = zi(k,j) - si * ar(i,k) + s * ai(i,k) 120 continue c 130 continue c 140 continue c 200 return end real*8 function epslon (x) real*8 x c c estimate unit roundoff in quantities of size x. c real*8 a,b,c,eps c c this program should function properly on all systems c satisfying the following two assumptions, c 1. the base used in representing floating point c numbers is not a power of three. c 2. the quantity a in statement 10 is represented to c the accuracy used in floating point variables c that are stored in memory. c the statement number 10 and the go to 10 are intended to c force optimizing compilers to generate code satisfying c assumption 2. c under these assumptions, it should be true that, c a is not exactly equal to four-thirds, c b has a zero for its last bit or digit, c c is not exactly equal to one, c eps measures the separation of 1.0 from c the next larger floating point number. c the developers of eispack would appreciate being informed c about any systems where these assumptions do not hold. c c this version dated 4/6/83. c a = 4.0d0/3.0d0 10 b = a - 1.0d0 c = b + b + b eps = dabs(c-1.0d0) if (eps .eq. 0.0d0) go to 10 epslon = eps*dabs(x) return end