PGI User Forum
 SearchSearch   MemberlistMemberlist     RegisterRegister   ProfileProfile    Log inLog in 

CUDA-x86.

Compute four precision on GPU but Permission denied why?
Goto page Previous  1, 2, 3, 4  Next
 
Post new topic   Reply to topic    PGI User Forum Forum Index -> Accelerator Programming
View previous topic :: View next topic  
Author Message
SWL_EGGBABY



Joined: 16 Dec 2009
Posts: 31

PostPosted: Wed Mar 17, 2010 12:24 am    Post subject: Compute four precision on GPU but Permission denied why? Reply with quote

module ddfunmod

contains

subroutine ddabrt
implicit none



stop
end subroutine

subroutine ddadd (dda, ddb, ddc)



implicit none
real*8 dda(2), ddb(2), ddc(2)
real*8 e, t1, t2



t1 = dda(1) + ddb(1)
e = t1 - dda(1)
t2 = ((ddb(1) - e) + (dda(1) - (t1 - e))) + dda(2) + ddb(2)



ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine

subroutine ddang (x, y, a)



implicit none
integer i, ix, iy, k, kk, nx, ny
real*8 t1, t2, t3
real*8 a(2), pi(2), x(2), y(2), s0(2), s1(2), s2(2), s3(2), s4(2)
save pi

data pi/ z'400921FB54442D18', z'3CA1A62633145C07'/


if (x(1) .eq. 0.d0 .and. y(1) .eq. 0.d0) then
write (6, 1)
1 format ('*** DDANG: Both arguments are zero.')
call ddabrt
return
endif


if (x(1) .eq. 0.d0) then
if (y(1) .gt. 0.d0) then
call ddmuld (pi, 0.5d0, a)
else
call ddmuld (pi, -0.5d0, a)
endif
goto 120
elseif (y(1) .eq. 0.d0) then
if (x(1) .gt. 0.d0) then
a(1) = 0.d0
a(2) = 0.d0
else
a(1) = pi(1)
a(2) = pi(2)
endif
goto 120
endif



call ddmul (x, x, s0)
call ddmul (y, y, s1)
call ddadd (s0, s1, s2)
call ddsqrt (s2, s3)
call dddiv (x, s3, s1)
call dddiv (y, s3, s2)



call ddqdc (s1, t1)
call ddqdc (s2, t2)
t3 = atan2 (t2, t1)
a(1) = t3
a(2) = 0.d0



if (abs (t1) .le. abs (t2)) then
kk = 1
s0(1) = s1(1)
s0(2) = s1(2)
else
kk = 2
s0(1) = s2(1)
s0(2) = s2(2)
endif


do k = 1, 3
call ddcssn (a, s1, s2)
if (kk .eq. 1) then
call ddsub (s0, s1, s3)
call dddiv (s3, s2, s4)
call ddsub (a, s4, s1)
else
call ddsub (s0, s2, s3)
call dddiv (s3, s1, s4)
call ddadd (a, s4, s1)
endif
a(1) = s1(1)
a(2) = s1(2)
enddo

120 continue

return
end subroutine

subroutine ddcadd (a, b, c)



implicit none
real*8 a(4), b(4), c(4)

call ddadd (a, b, c)
call ddadd (a(3), b(3), c(3))

return
end subroutine

subroutine ddcdiv (a, b, c)



implicit none
real*8 a(4), b(4), c(4), f(2), s0(2), s1(2), s2(2), s3(2), s4(2)

if (b(1) .eq. 0.d0 .and. b(3) .eq. 0.d0) then
write (6, 1)
1 format ('*** DDCDIV: Divisor is zero.')
call ddabrt
return
endif

f(1) = 1.d0
f(2) = 0.d0
call ddmul (a, b, s0)
call ddmul (a(3), b(3), s1)
call ddadd (s0, s1, s2)
call ddsub (s0, s1, s3)
call ddadd (a, a(3), s0)
call ddsub (b, b(3), s1)
call ddmul (s0, s1, s4)
call ddsub (s4, s3, s1)
call ddmul (b, b, s0)
call ddmul (b(3), b(3), s3)
call ddadd (s0, s3, s4)
call dddiv (f, s4, s0)
call ddmul (s2, s0, c)
call ddmul (s1, s0, c(3))

return
end subroutine

subroutine ddceq (a, b)



implicit none
real*8 a(4), b(4)

b(1) = a(1)
b(2) = a(2)
b(3) = a(3)
b(4) = a(4)

return
end subroutine

subroutine ddcmul (a, b, c)



implicit none
real*8 a(4), b(4), c(4), s0(2), s1(2), s2(2), s3(2)

call ddmul (a, b, s0)
call ddmul (a(3), b(3), s1)
call ddmul (a, b(3), s2)
call ddmul (a(3), b, s3)
call ddsub (s0, s1, c)
call ddadd (s2, s3, c(3))

return
end subroutine

subroutine ddcpr (a, b, ic)


implicit none
integer ic
real*8 a(2), b(2)

if (a(1) .lt. b(1)) then
ic = -1
elseif (a(1) .eq. b(1)) then
if (a(2) .lt. b(2)) then
ic = -1
elseif (a(2) .eq. b(2)) then
ic = 0
else
ic = 1
endif
else
ic = 1
endif

return
end subroutine

subroutine ddcpwr (a, n, b)



implicit none
integer j, kk, kn, l1, mn, n, na1, na2, nn
real*8 cl2, t1
parameter (cl2 = 1.4426950408889633d0)
real*8 a(4), b(4), s0(4), s1(4), s2(4), s3(4)

if (a(1) .eq. 0.d0 .and. a(3) .eq. 0.d0) then
if (n .ge. 0) then
b(1) = 0.d0
b(2) = 0.d0
b(3) = 0.d0
b(4) = 0.d0
goto 120
else
write (6, 1)
1 format ('*** DDCPWR: Argument is zero and N is negative or zero.')
call ddabrt
return
endif
endif

nn = abs (n)
if (nn .eq. 0) then
s2(1) = 1.d0
s2(2) = 0.d0
s2(3) = 0.d0
s2(4) = 0.d0
goto 120
elseif (nn .eq. 1) then
s2(1) = a(1)
s2(2) = a(2)
s2(3) = a(3)
s2(4) = a(4)
goto 110
elseif (nn .eq. 2) then
call ddcmul (a, a, s2)
goto 110
endif


t1 = nn
mn = cl2 * log (t1) + 1.d0 + 1.d-14

s0(1) = a(1)
s0(2) = a(2)
s0(3) = a(3)
s0(4) = a(4)
s2(1) = 1.d0
s2(2) = 0.d0
s2(3) = 0.d0
s2(4) = 0.d0
kn = nn


do j = 1, mn
kk = kn / 2
if (kn .ne. 2 * kk) then
call ddcmul (s2, s0, s1)
s2(1) = s1(1)
s2(2) = s1(2)
s2(3) = s1(3)
s2(4) = s2(4)
endif
kn = kk
if (j .lt. mn) then
call ddcmul (s0, s0, s1)
s0(1) = s1(1)
s0(2) = s1(2)
s0(3) = s1(3)
s0(4) = s1(4)
endif
enddo

110 continue

if (n .lt. 0) then
s1(1) = 1.d0
s1(2) = 0.d0
s1(3) = 0.d0
s1(4) = 0.d0
call ddcdiv (s1, s2, s0)
s2(1) = s0(1)
s2(2) = s0(2)
s2(3) = s0(3)
s2(4) = s0(4)
endif

b(1) = s2(1)
b(2) = s2(2)
b(3) = s2(3)
b(4) = s2(4)

120 continue
return
end subroutine

subroutine ddcsqrt (a, b)



implicit none
real*8 a(4), b(4), s0(2), s1(2), s2(2)

if (a(1) .eq. 0.d0 .and. a(3) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
b(3) = 0.d0
b(4) = 0.d0
goto 100
endif

call ddmul (a, a, s0)
call ddmul (a(3), a(3), s1)
call ddadd (s0, s1, s2)
call ddsqrt (s2, s0)
s1(1) = a(1)
s1(2) = a(2)
if (s1(1) .lt. 0.d0) then
s1(1) = - s1(1)
s1(2) = - s1(2)
endif
s1 = abs (s1)
call ddadd (s0, s1, s2)
call ddmuld (s2, 0.5d0, s1)
call ddsqrt (s1, s0)
call ddmuld (s0, 2.d0, s1)
if (a(1) .ge. 0.d0) then
b(1) = s0(1)
b(2) = s0(2)
call dddiv (a(3), s1, b(3))
else
call dddiv (a(3), s1, b)
if (b(1) .lt. 0.d0) then
b(1) = - b(1)
b(2) = - b(2)
endif
b(3) = s0(1)
b(4) = s0(2)
if (a(3) .lt. 0.d0) then
b(3) = - b(3)
b(4) = - b(4)
endif
endif

100 continue
return
end subroutine

subroutine ddcssh (a, x, y)


implicit none
real*8 a(2), f(2), x(2), y(2), s0(2), s1(2), s2(2)

f(1) = 1.d0
f(2) = 0.d0
call ddexp (a, s0)
call dddiv (f, s0, s1)
call ddadd (s0, s1, s2)
call ddmuld (s2, 0.5d0, x)
call ddsub (s0, s1, s2)
call ddmuld (s2, 0.5d0, y)

return
end subroutine

subroutine ddcssn (a, x, y)



implicit none
integer ia, ka, kb, kc, l1
real*8 t1, t2
real*8 a(2), f(2), pi(2), x(2), y(2), s0(2), s1(2), s2(2), s3(2), s4(2), &
s5(2), s6(2)
real*8 cs(2,2,4)
save cs, pi


data cs/ &
z'3FEF6297CFF75CB0', z'3C7562172A361FD3', &
z'3FC8F8B83C69A60B', z'BC626D19B9FF8D82', &
z'3FED906BCF328D46', z'3C7457E610231AC2', &
z'3FD87DE2A6AEA963', z'BC672CEDD3D5A610', &
z'3FEA9B66290EA1A3', z'3C39F630E8B6DAC8', &
z'3FE1C73B39AE68C8', z'3C8B25DD267F6600', &
z'3FE6A09E667F3BCD', z'BC8BDD3413B26456', &
z'3FE6A09E667F3BCD', z'BC8BDD3413B26456' /
data pi/ z'400921FB54442D18', z'3CA1A62633145C07'/



if (a(1) .eq. 0.d0) then
x(1) = 1.d0
x(2) = 0.d0
y(1) = 0.d0
y(2) = 0.d0
goto 120
endif

f(1) = 1.d0
f(2) = 0.d0


call ddmuld (pi, 2.d0, s0)
call dddiv (a, s0, s1)
call ddnint (s1, s2)
call ddsub (s1, s2, s3)

t1 = s3(1)
t2 = 4.d0 * t1
ka = nint (t2)
kb = nint (8.d0 * (t2 - ka))
t1 = (8 * ka + kb) / 32.d0
s1(1) = t1
s1(2) = 0.d0
call ddsub (s3, s1, s2)
call ddmul (s0, s2, s1)


if (s1(1) .eq. 0.d0) then
s0(1) = 0.d0
s0(2) = 0.d0
goto 110
endif
s0(1) = s1(1)
s0(2) = s1(2)
call ddmul (s0, s0, s2)
l1 = 0

100 l1 = l1 + 1
if (l1 .eq. 100) then
write (6, 1)
1 format ('*** DDCSSN: Iteration limit exceeded.')
call ddabrt
return
endif

t2 = - (2.d0 * l1) * (2.d0 * l1 + 1.d0)
call ddmul (s2, s1, s3)
call dddivd (s3, t2, s1)
call ddadd (s1, s0, s3)
s0(1) = s3(1)
s0(2) = s3(2)


if (abs (s1(1)) .gt. 1d-33 * abs (s3(1))) goto 100


110 continue
s1(1) = s0(1)
s1(2) = s0(2)
call ddmul (s0, s0, s2)
call ddsub (f, s2, s3)
call ddsqrt (s3, s0)


kc = abs (kb)
f(1) = 2.
if (kc .eq. 0) then
s2(1) = 1.d0
s2(2) = 0.d0
s3(1) = 0.d0
s3(2) = 0.d0
else
s2(1) = cs(1,1,kc)
s2(2) = cs(2,1,kc)
s3(1) = cs(1,2,kc)
s3(2) = cs(2,2,kc)
endif
if (kb .lt. 0) then
s3(1) = - s3(1)
s3(2) = - s3(2)
endif

call ddmul (s0, s2, s4)
call ddmul (s1, s3, s5)
call ddsub (s4, s5, s6)
call ddmul (s1, s2, s4)
call ddmul (s0, s3, s5)
call ddadd (s4, s5, s1)
s0(1) = s6(1)
s0(2) = s6(2)


if (ka .eq. 0) then
x(1) = s0(1)
x(2) = s0(2)
y(1) = s1(1)
y(2) = s1(2)
elseif (ka .eq. 1) then
x(1) = - s1(1)
x(2) = - s1(2)
y(1) = s0(1)
y(2) = s0(2)
elseif (ka .eq. -1) then
x(1) = s1(1)
x(2) = s1(2)
y(1) = - s0(1)
y(2) = - s0(2)
elseif (ka .eq. 2 .or. ka .eq. -2) then
x(1) = - s0(1)
x(2) = - s0(2)
y(1) = - s1(1)
y(2) = - s1(2)
endif

120 continue
return
end subroutine

subroutine ddcsub (a, b, c)

implicit none
real*8 a(4), b(4), c(4)

call ddsub (a, b, c)
call ddsub (a(3), b(3), c(3))

return
end subroutine

subroutine dddiv (dda, ddb, ddc)


implicit none
real*8 dda(2), ddb(2), ddc(2)
real*8 a1, a2, b1, b2, cona, conb, c11, c2, c21, e, split, s1, s2, &
t1, t2, t11, t12, t21, t22
parameter (split = 134217729.d0)

s1 = dda(1) / ddb(1)

cona = s1 * split
conb = ddb(1) * split
a1 = cona - (cona - s1)
b1 = conb - (conb - ddb(1))
a2 = s1 - a1
b2 = ddb(1) - b1


c11 = s1 * ddb(1)
c21 = (((a1 * b1 - c11) + a1 * b2) + a2 * b1) + a2 * b2

c2 = s1 * ddb(2)


t1 = c11 + c2
e = t1 - c11
t2 = ((c2 - e) + (c11 - (t1 - e))) + c21


t12 = t1 + t2
t22 = t2 - (t12 - t1)


t11 = dda(1) - t12
e = t11 - dda(1)
t21 = ((-t12 - e) + (dda(1) - (t11 - e))) + dda(2) - t22


s2 = (t11 + t21) / ddb(1)


ddc(1) = s1 + s2
ddc(2) = s2 - (ddc(1) - s1)

return
end subroutine

subroutine dddivd (dda, db, ddc)

implicit none
real*8 dda(2), db, ddc(2)
real*8 a1, a2, b1, b2, cona, conb, e, split, t1, t2, t11, t12, t21, t22
parameter (split = 134217729.d0)


t1 = dda(1) / db

cona = t1 * split
conb = db * split
a1 = cona - (cona - t1)
b1 = conb - (conb - db)
a2 = t1 - a1
b2 = db - b1


t12 = t1 * db
t22 = (((a1 * b1 - t12) + a1 * b2) + a2 * b1) + a2 * b2


t11 = dda(1) - t12
e = t11 - dda(1)
t21 = ((-t12 - e) + (dda(1) - (t11 - e))) + dda(2) - t22



t2 = (t11 + t21) / db


ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine

subroutine dddqc (a, b)


implicit none
real*8 a, b(2)

b(1) = a
b(2) = 0.d0
return
end subroutine

subroutine ddeq (a, b)


implicit none
real*8 a(2), b(2)

b(1) = a(1)
b(2) = a(2)
end subroutine

subroutine ddeform (a, n1, n2, b)


implicit none
integer i, ln, m1, n1, n2
parameter (ln = 40)
character*1 b(n1)
character*40 cs
real*8 a(2)

if (n1 .lt. 0 .or. n2 .lt. 0 .or. n1 .gt. 80 .or. n2 .gt. 30 &
.or. n2 .gt. n1 - 8) then
write (6, 1) n1, n2
1 format ('*** DDEFORM: Improper n1, n2 =',2i6)
call ddabrt
endif

call ddoutc (a, cs)
m1 = n1 - n2 - 8

do i = 1, m1
b(i) = ' '
enddo

do i = 1, n2 + 3
b(i+m1) = cs(i+2:i+2)
enddo

do i = 1, 5
b(i+m1+n2+3) = cs(i+35:i+35)
enddo

return
end subroutine

subroutine ddexp (a, b)



implicit none
integer i, ia, l1, na, nq, nz, n1
real*8 t1, t2
parameter (nq = 6)
real*8 a(2), b(2), al2(2), f(2), s0(2), s1(2), s2(2), s3(2), tl
save al2
!>
! Uncomment one of the following two lines, preferably the first if it is
! accepted by the compiler.

data al2 / z'3FE62E42FEFA39EF', z'3C7ABC9E3B39803F'/
! data al2/ 6.9314718055994529D-01, 2.3190468138462996D-17/

! Check for overflows and underflows.

if (abs (a(1)) .ge. 709.d0) then
if (a(1) .gt. 0.d0) then
write (6, 1) a(1)
1 format ('*** DDEXP: Argument is too large',f12.6)
call ddabrt
return
else
call dddqc (0.d0, b)
goto 130
endif
endif

f(1) = 1.d0
f(2) = 0.d0

! Compute the reduced argument A' = A - Log(2) * Nint [A / Log(2)]. Save
! NZ = Nint [A / Log(2)] for correcting the exponent of the final result.

call dddiv (a, al2, s0)
call ddnint (s0, s1)
t1 = s1(1)
nz = t1 + sign (1.d-14, t1)
call ddmul (al2, s1, s2)
call ddsub (a, s2, s0)

! Check if the reduced argument is zero.

if (s0(1) .eq. 0.d0) then
s0(1) = 1.d0
s0(2) = 0.d0
l1 = 0
goto 120
endif

! Divide the reduced argument by 2 ^ NQ.

call dddivd (s0, 2.d0 ** nq, s1)

! Compute Exp using the usual Taylor series.

s2(1) = 1.d0
s2(2) = 0.d0
s3(1) = 1.d0
s3(2) = 0.d0
l1 = 0

100 l1 = l1 + 1
if (l1 .eq. 100) then
write (6, 2)
2 format ('*** DDEXP: Iteration limit exceeded.')
call ddabrt
return
endif

t2 = l1
call ddmul (s2, s1, s0)
call dddivd (s0, t2, s2)
call ddadd (s3, s2, s0)
call ddeq (s0, s3)

! Check for convergence of the series.

if (abs (s2(1)) .gt. 1d-33 * abs (s3(1))) goto 100

! Raise to the (2 ^ NQ)-th power.

do i = 1, nq
call ddmul (s0, s0, s1)
s0(1) = s1(1)
s0(2) = s1(2)
enddo

! Multiply by 2 ^ NZ.

120 call ddmuld (s0, 2.d0 ** nz, b)

! Restore original precision level.

130 continue
return
end subroutine

subroutine ddfform (a, n1, n2, b)

! This routine converts the DD number A to F format, i.e. Fn1.n2.
! B is the output array (type character*1 of size n1). N1 must exceed
! N2 by at least 3, and N1 must not exceed 80. N2 must not exceed 30.

implicit none
integer i, ix, kx, ln, ls, lz, mx, nx, n1, n2
parameter (ln = 40)
real*8 a(2)
character*1 b(n1)
character*40 c
character*40 chr40

if (n1 .lt. 0 .or. n2 .lt. 0 .or. n1 .gt. 80 .or. n2 .gt. 30 &
.or. n1 - n2 .lt. 3) then
write (6, 1) n1, n2
1 format ('*** DDFFORM: Improper n1, n2 =',2i6)
call ddabrt
endif

! Call ddoutc and extract exponent.

call ddoutc (a, c)
ix = dddigin (c(ln-3:ln), 4)

do i = 1, n1
b(i) = ' '
enddo

if (a(1) .ge. 0.d0) then
ls = 0
else
ls = 1
endif
if (ix .ge. 0 .and. a(1) .ne. 0.d0) then
lz = 0
else
lz = 1
endif
mx = max (ix, 0)

! Check for overflow of field length.

if (ls + lz + mx + n2 + 2 .gt. n1) then
do i = 1, n1
b(i) = '*'
enddo

goto 200
endif

! Check if a zero should be output.

if (a(1) .eq. 0. .or. -ix .gt. n2) then
do i = 1, n1 - n2 - 2
b(i) = ' '
enddo

b(n1-n2-1) = '0'
b(n1-n2) = '.'

do i = 1, n2
b(i+n1-n2) = '0'
enddo

goto 200
endif

! Process other cases.

if (a(1) .lt. 0.) b(n1-n2-mx-2) = '-'
if (ix .ge. 0) then
b(n1-n2-ix-1) = c(4:4)
kx = min (ln - 9, ix)

do i = 1, kx
b(i+n1-n2-ix-1) = c(i+5:i+5)
enddo

do i = kx + 1, ix
b(i+n1-n2-ix-1) = '0'
enddo

b(n1-n2) = '.'
kx = max (min (ln - 9 - ix, n2), 0)

do i = 1, kx
b(i+n1-n2) = c(i+ix+5:i+ix+5)
enddo

do i = kx + 1, n2
b(i+n1-n2) = '0'
enddo
else
nx = - ix
b(n1-n2-1) = '0'
b(n1-n2) = '.'

do i = 1, nx - 1
b(i+n1-n2) = '0'
enddo

b(n1-n2+nx) = c(4:4)
kx = min (ln - 8, n2 - nx)

do i = 1, kx
b(i+n1-n2+nx) = c(i+5:i+5)
enddo

do i = kx + 1, n2 - nx
b(i+n1-n2+nx) = '0'
enddo
endif

200 continue

return
end subroutine

subroutine ddinfr (a, b, c)

! Sets B to the integer part of the DD number A and sets C equal to the
! fractional part of A. Note that if A = -3.3, then B = -3 and C = -0.3.

implicit none
integer ic
real*8 a(2), b(2), c(2), con(2), f(2), s0(2), s1(2), t105, t52
parameter (t105 = 2.d0 ** 105, t52 = 2.d0 ** 52)
save con
data con / t105, t52/

! Check if A is zero.

if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
c(1) = 0.d0
c(2) = 0.d0
goto 120
endif

if (a(1) .ge. t105) then
write (6, 1)
1 format ('*** DDINFR: Argument is too large.')
call ddabrt
return
endif

f(1) = 1.d0
f(2) = 0.d0
if (a(1) .gt. 0.d0) then
call ddadd (a, con, s0)
call ddsub (s0, con, b)
call ddcpr (a, b, ic)
if (ic .ge. 0) then
call ddsub (a, b, c)
else
call ddsub (b, f, s1)
b(1) = s1(1)
b(2) = s1(2)
call ddsub (a, b, c)
endif
else
call ddsub (a, con, s0)
call ddadd (s0, con, b)
call ddcpr (a, b, ic)
if (ic .le. 0) then
call ddsub (a, b, c)
else
call ddadd (b, f, s1)
b(1) = s1(1)
b(2) = s1(2)
call ddsub (a, b, c)
endif
endif

120 continue
return
end subroutine

subroutine ddinp (iu, a)

! This routine reads the DD number A from logical unit IU. The input
! value must be placed on a single line of not more than 80 characters.

implicit none
integer iu, ln
parameter (ln = 80)
character*80 cs
real*8 a(2)

read (iu, '(a)', end = 100) cs
call ddinpc (cs, a)
goto 110

100 write (6, 1)
1 format ('*** DDINP: End-of-file encountered.')
call ddabrt

110 return
end subroutine

subroutine ddinpc (a, b)

! Converts the CHARACTER*80 array A into the DD number B.

implicit none
integer i, ib, id, ie, inz, ip, is, ix, k, ln, lnn
parameter (ln = 80)
real*8 bi
character*80 a
character*1 ai
character*10 ca, dig
parameter (dig = '0123456789')
real*8 b(2), f(2), s0(2), s1(2), s2(2)

id = 0
ip = -1
is = 0
inz = 0
s1(1) = 0.d0
s1(2) = 0.d0

do i = 80, 1, -1
if (a(i:i) /= ' ') goto 90
enddo

90 continue

lnn = i

! Scan for digits, looking for the period also.

do i = 1, lnn
ai = a(i:i)
if (ai .eq. ' ' .and. id == 0) then
elseif (ai .eq. '.') then
if (ip >= 0) goto 210
ip = id
inz = 1
elseif (ai .eq. '+') then
if (id .ne. 0 .or. ip >= 0 .or. is .ne. 0) goto 210
is = 1
elseif (ai .eq. '-') then
if (id .ne. 0 .or. ip >= 0 .or. is .ne. 0) goto 210
is = -1
elseif (ai .eq. 'e' .or. ai .eq. 'E' .or. ai .eq. 'd' .or. ai .eq. 'D') then
goto 100
elseif (index (dig, ai) .eq. 0) then
goto 210
else
! read (ai, '(f1.0)') bi
bi = index (dig, ai) - 1
if (inz > 0 .or. bi > 0.d0) then
inz = 1
id = id + 1
call ddmuld (s1, 10.d0, s0)
f(1) = bi
f(2) = 0.d0
call dddqc (bi, f)
call ddadd (s0, f, s1)
endif
endif
enddo

100 continue
if (is .eq. -1) then
s1(1) = - s1(1)
s1(2) = - s1(2)
endif
k = i
if (ip == -1) ip = id
ie = 0
is = 0
ca = ' '

do i = k + 1, lnn
ai = a(i:i)
if (ai .eq. ' ') then
elseif (ai .eq. '+') then
if (ie .ne. 0 .or. is .ne. 0) goto 210
is = 1
elseif (ai .eq. '-') then
if (ie .ne. 0 .or. is .ne. 0) goto 210
is = -1
elseif (index (dig, ai) .eq. 0) then
goto 210
else
ie = ie + 1
if (ie .gt. 3) goto 210
ca(ie:ie) = ai
endif
enddo

! read (ca, '(i4)') ie
ie = dddigin (ca, 4)
if (is .eq. -1) ie = - ie
ie = ie + ip - id
s0(1) = 10.d0
s0(2) = 0.d0
call ddnpwr (s0, ie, s2)
call ddmul (s1, s2, b)
goto 220

210 write (6, 1)
1 format ('*** DDINPC: Syntax error in literal string.')
call ddabrt

220 continue

return
end subroutine

subroutine ddlog (a, b)



implicit none
integer k
real*8 t1, t2
real*8 a(2), al2(2), b(2), s0(2), s1(2), s2(2)
save al2
!>
! Uncomment one of the following two lines, preferably the first if it is
! accepted by the compiler.

data al2 / z'3FE62E42FEFA39EF', z'3C7ABC9E3B39803F'/
! data al2/ 6.9314718055994529D-01, 2.3190468138462996D-17/

if (a(1) .le. 0.d0) then
write (6, 1)
1 format ('*** DDLOG: Argument is less than or equal to zero.')
call ddabrt
return
endif

! Compute initial approximation of Log (A).

t1 = a(1)
t2 = log (t1)
b(1) = t2
b(2) = 0.d0

! Perform the Newton-Raphson iteration described above.

do k = 1, 3
call ddexp (b, s0)
call ddsub (a, s0, s1)
call dddiv (s1, s0, s2)
call ddadd (b, s2, s1)
b(1) = s1(1)
b(2) = s1(2)
enddo

120 continue

return
end subroutine

subroutine ddqdc (a, b)

! This converts the DD number A to DP.

implicit none
real*8 a(2), b

b = a(1)
return
end subroutine

subroutine ddqqc (a, b, c)

! This converts DD numbers A and B to DDC form in C, i.e. C = A + B i.

implicit none
real*8 a(2), b(2), c(4)

c(1) = a(1)
c(2) = a(2)
c(3) = b(1)
c(4) = b(2)
return
end subroutine

subroutine ddmul (dda, ddb, ddc)

! This routine multiplies DD numbers DDA and DDB to yield the DD product DDC.

implicit none
real*8 dda(2), ddb(2), ddc(2)
real*8 a1, a2, b1, b2, cona, conb, c11, c21, c2, e, split, t1, t2
parameter (split = 134217729.d0)



cona = dda(1) * split
conb = ddb(1) * split
a1 = cona - (cona - dda(1))
b1 = conb - (conb - ddb(1))
a2 = dda(1) - a1
b2 = ddb(1) - b1

! Multilply dda(1) * ddb(1) using Dekker's method.

c11 = dda(1) * ddb(1)
c21 = (((a1 * b1 - c11) + a1 * b2) + a2 * b1) + a2 * b2
!>
! Compute dda(1) * ddb(2) + dda(2) * ddb(1) (only high-order word is needed).

c2 = dda(1) * ddb(2) + dda(2) * ddb(1)

! Compute (c11, c21) + c2 using Knuth's trick, also adding low-order product.

t1 = c11 + c2
e = t1 - c11
t2 = ((c2 - e) + (c11 - (t1 - e))) + c21 + dda(2) * ddb(2)

! The result is t1 + t2, after normalization.

ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)

return
end subroutine

subroutine ddmuld (dda, db, ddc)

! This routine multiplies the DD number DDA by the DP number DB to yield
! the DD product DDC.

implicit none
real*8 dda(2), db, ddc(2)
real*8 a1, a2, b1, b2, cona, conb, c11, c21, c2, e, split, t1, t2
parameter (split = 134217729.d0)



cona = dda(1) * split
conb = db * split
a1 = cona - (cona - dda(1))
b1 = conb - (conb - db)
a2 = dda(1) - a1
b2 = db - b1

! Multilply dda(1) * db using Dekker's method.

c11 = dda(1) * db
c21 = (((a1 * b1 - c11) + a1 * b2) + a2 * b1) + a2 * b2
!>
! Compute dda(2) * db (only high-order word is needed).

c2 = dda(2) * db

! Compute (c11, c21) + c2 using Knuth's trick.

t1 = c11 + c2
e = t1 - c11
t2 = ((c2 - e) + (c11 - (t1 - e))) + c21

! The result is t1 + t2, after normalization.

ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine

subroutine ddmuldd (da, db, ddc)

! This subroutine computes ddc = da x db.

implicit none
real*8 a1, a2, b1, b2, cona, conb, da, db, ddc(2), split, s1, s2
parameter (split = 134217729.d0)


cona = da * split
conb = db * split
a1 = cona - (cona - da)
b1 = conb - (conb - db)
a2 = da - a1
b2 = db - b1

! Multiply da * db using Dekker's method.

s1 = da * db
s2 = (((a1 * b1 - s1) + a1 * b2) + a2 * b1) + a2 * b2
!>
ddc(1) = s1
ddc(2) = s2

return
end subroutine

subroutine ddnint (a, b)

! This sets B equal to the integer nearest to the DD number A.

implicit none
real*8 a(2), b(2), con(2), s0(2), t105, t52
parameter (t105 = 2.d0 ** 105, t52 = 2.d0 ** 52)
save con
data con / t105, t52/

! Check if A is zero.

if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
goto 120
endif

if (a(1) .ge. t105) then
write (6, 1)
1 format ('*** DDINFR: Argument is too large.')
call ddabrt
return
endif

if (a(1) .gt. 0.d0) then
call ddadd (a, con, s0)
call ddsub (s0, con, b)
else
call ddsub (a, con, s0)
call ddadd (s0, con, b)
endif

120 continue
return
end subroutine

subroutine ddnpwr (a, n, b)


implicit none
integer j, kk, kn, l1, mn, n, na1, na2, nn
real*8 cl2, t1
parameter (cl2 = 1.4426950408889633d0)
real*8 a(2), b(2), s0(2), s1(2), s2(2), s3(2)

if (a(1) .eq. 0.d0) then
if (n .ge. 0) then
s2(1) = 0.d0
s2(2) = 0.d0
goto 120
else
write (6, 1)
1 format ('*** DDCPWR: Argument is zero and N is negative or zero.')
call ddabrt
return
endif
endif

nn = abs (n)
if (nn .eq. 0) then
s2(1) = 1.d0
s2(2) = 0.d0
goto 120
elseif (nn .eq. 1) then
s2(1) = a(1)
s2(2) = a(2)
goto 110
elseif (nn .eq. 2) then
call ddmul (a, a, s2)
goto 110
endif

! Determine the least integer MN such that 2 ^ MN .GT. NN.

t1 = nn
mn = cl2 * log (t1) + 1.d0 + 1.d-14
s0(1) = a(1)
s0(2) = a(2)
s2(1) = 1.d0
s2(2) = 0.d0
kn = nn

! Compute B ^ N using the binary rule for exponentiation.

do j = 1, mn
kk = kn / 2
if (kn .ne. 2 * kk) then
call ddmul (s2, s0, s1)
s2(1) = s1(1)
s2(2) = s1(2)
endif
kn = kk
if (j .lt. mn) then
call ddmul (s0, s0, s1)
s0(1) = s1(1)
s0(2) = s1(2)
endif
enddo

! Compute reciprocal if N is negative.

110 continue

if (n .lt. 0) then
s1(1) = 1.d0
s1(2) = 0.d0
call dddiv (s1, s2, s0)
s2(1) = s0(1)
s2(2) = s0(2)
endif

120 continue

b(1) = s2(1)
b(2) = s2(2)

return
end subroutine

subroutine ddnrt (a, n, b)


implicit none
integer i, k, n
real*8 t1, t2, tn
real*8 a(2), b(2), f1(2), s0(2), s1(2)

if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
goto 140
elseif (a(1) .lt. 0.d0) then
write (6, 1)
1 format ('*** DDNRT: Argument is negative.')
call ddabrt
return
endif
if (n .le. 0) then
write (6, 2) n
2 format ('*** DDNRT: Improper value of N',i10)
call ddabrt
return
endif

! Handle cases N = 1 and 2.

if (n .eq. 1) then
b(1) = a(1)
b(2) = a(1)
goto 140
elseif (n .eq. 2) then
call ddsqrt (a, b)
goto 140
endif

f1(1) = 1.d0
f1(2) = 0.d0

! Compute the initial approximation of A ^ (-1/N).

tn = n
t1 = a(1)
t2 = exp (- log (t1) / tn)
b(1) = t2
b(2) = 0.d0

! Perform the Newton-Raphson iteration described above.

do k = 1, 3
call ddnpwr (b, n, s0)
call ddmul (a, s0, s1)
call ddsub (f1, s1, s0)
call ddmul (b, s0, s1)
call dddivd (s1, tn, s0)
call ddadd (b, s0, s1)
b(1) = s1(1)
b(2) = s1(2)
enddo

! Take the reciprocal to give final result.

call dddiv (f1, b, s1)
b(1) = s1(1)
b(2) = s1(2)

140 continue
return
end subroutine

subroutine ddout (iu, a)

! This routine writes the DD number A on logical unit iu using a standard
! E format, with lines 40 characters long.

implicit none
integer i, iu, ln
parameter (ln = 40)
character*40 cs
real*8 a(2)

call ddoutc (a, cs)
write (iu, '(a)') cs

return
end subroutine

subroutine ddoutc (a, b)


implicit none
integer i, ii, ix, ln, nx
parameter (ln = 40)
integer ib(ln)
real*8 t1
character*40 b
character*10 ca, digits
parameter (digits = '0123456789')
real*8 a(2), f(2), s0(2), s1(2)

f(1) = 10.d0
f(2) = 0.d0

do i = 1, ln
ib(i) = 0
enddo

! Determine exact power of ten for exponent.

if (a(1) .ne. 0.d0) then
t1 = log10 (abs (a(1)))
if (t1 .ge. 0.d0) then
nx = t1
else
nx = t1 - 1.d0
endif
call ddnpwr (f, nx, s0)
call dddiv (a, s0, s1)
if (s1(1) .lt. 0.d0) then
s1(1) = - s1(1)
s1(2) = - s1(2)
endif

! If we didn't quite get it exactly right, multiply or divide by 10 to fix.

i = 0

100 continue

i = i + 1
if (s1(1) .lt. 1.d0) then
nx = nx - 1
call ddmuld (s1, 10.d0, s0)
s1(1) = s0(1)
s1(2) = s0(2)
if (i <= 3) goto 100
elseif (s1(1) .ge. 10.d0) then
nx = nx + 1
call dddivd (s1, 10.d0, s0)
s1(1) = s0(1)
s1(2) = s0(2)
goto 100
endif
else
nx = 0
s1(1) = 0.d0
s1(2) = 0.d0
endif

! Compute digits.

do i = 1, ln - 8
ii = s1(1)
ib(i) = ii
f(1) = ii
call ddsub (s1, f, s0)
call ddmuld (s0, 10.d0, s1)
enddo

! Fix negative digits.

do i = ln - 8, 2, -1
if (ib(i) .lt. 0) then
ib(i) = ib(i) + 10
ib(i-1) = ib(i-1) - 1
endif
enddo

if (ib(1) .lt. 0) then
write (6, 1)
1 format ('ddoutc: negative leading digit')
call ddabrt
endif

! Round.

if (ib(ln-8) .ge. 5) then
ib(ln-9) = ib(ln-9) + 1

do i = ln - 9, 2, -1
if (ib(i) .eq. 10) then
ib(i) = 0
ib(i-1) = ib(i-1) + 1
endif
enddo

if (ib(1) .eq. 10) then
ib(1) = 1
nx = nx + 1
endif
endif

! Insert digit characters in ib.

b(1:1) = ' '
b(2:2) = ' '
if (a(1) .ge. 0.d0) then
b(3:3) = ' '
else
b(3:3) = '-'
endif
ii = ib(1)
b(4:4) = digits(ii+1:ii+1)
b(5:5) = '.'
b(ln:ln) = ' '

do i = 2, ln - 9
ii = ib(i)
b(i+4:i+4) = digits(ii+1:ii+1)
enddo

! Insert exponent.

190 continue
! write (ca, '(i4)') nx
ca = dddigout (dble (nx), 4)
b(ln-4:ln-4) = 'E'
ii = 0

do i = 1, 4
if (ca(i:i) /= ' ') then
ii = ii + 1
b(ln-4+ii:ln-4+ii) = ca(i:i)
endif
enddo

do i = ii + 1, 4
b(ln-4+i:ln-4+i) = ' '
enddo

return
end subroutine

subroutine ddpic (pi)

! This returns Pi to quad precision.

implicit none
real*8 pi(2), pic(2)
save pic
!>
! Uncomment one of the following two lines, preferably the first if it is
! accepted by the compiler.

data pic/ z'400921FB54442D18', z'3CA1A62633145C07'/
! data pic/ 3.1415926535897931D+00, 1.2246467991473532D-16/

pi(1) = pic(1)
pi(2) = pic(2)

return
end subroutine

subroutine ddpoly (n, a, x0, x)


implicit none
integer i, it, n
real*8 a(2,0:n), ad(2,0:n), t1(2), t2(2), t3(2), t4(2), t5(2), &
x(2), x0(2), dt1, eps
parameter (eps = 1.d-30)

do i = 0, n - 1
dt1 = i + 1
call ddmuld (a(1,i+1), dt1, ad(1,i))
enddo

ad(1,n) = 0.d0
ad(2,n) = 0.d0
x(1) = x0(1)
x(2) = x0(2)

do it = 1, 20
t1(1) = 0.d0
t1(2) = 0.d0
t2(1) = 0.d0
t2(2) = 0.d0
t3(1) = 1.d0
t3(2) = 0.d0

do i = 0, n
call ddmul (a(1,i), t3, t4)
call ddadd (t1, t4, t5)
t1(1) = t5(1)
t1(2) = t5(2)
call ddmul (ad(1,i), t3, t4)
call ddadd (t2, t4, t5)
t2(1) = t5(1)
t2(2) = t5(2)
call ddmul (t3, x, t4)
t3(1) = t4(1)
t3(2) = t4(2)
enddo

call dddiv (t1, t2, t3)
call ddsub (x, t3, t4)
x(1) = t4(1)
x(2) = t4(2)
if (abs (t3(1)) .le. eps) goto 110
enddo

write (6, 1)
1 format ('ddroot: failed to converge.')
call ddabrt

110 continue

return
end subroutine

subroutine ddrand (a)

! This returns a pseudo-random DD number A between 0 and 1.

implicit none
real*8 f7, r28, r30, r53, r58, s0, s1, s2, sd, t1, t2, t30
parameter (f7 = 78125.d0, s0 = 314159265.d0, r30 = 0.5d0 ** 30, &
r53 = 0.5d0 ** 53, r58 = 0.5d0 ** 58, t30 = 2.d0 ** 30)
real*8 a(2)
save sd
data sd /s0/

t1 = f7 * sd
t2 = aint (r30 * t1)
s1 = t1 - t30 * t2
t1 = f7 * s1
t2 = aint (r30 * t1)
s2 = t1 - t30 * t2
a(1) = r30 * s1 + r58 * s2
t1 = f7 * s2
t2 = aint (r30 * t1)
s1 = t1 - t30 * t2
t1 = f7 * s1
t2 = aint (r30 * t1)
s2 = t1 - t30 * t2
a(2) = r53 * a(1) * (r30 * s1 + r58 * s2)
sd = s2

return
end subroutine

subroutine ddsqrt (a, b)



implicit none
real*8 t1, t2, t3
real*8 a(2), b(2), f(2), s0(2), s1(2)

if (a(1) .eq. 0.d0) then
b(1) = 0.d0
b(2) = 0.d0
goto 100
endif
t1 = 1.d0 / sqrt (a(1))
t2 = a(1) * t1
call ddmuldd (t2, t2, s0)
call ddsub (a, s0, s1)
t3 = 0.5d0 * s1(1) * t1
s0(1) = t2
s0(2) = 0.d0
s1(1) = t3
s1(2) = 0.d0
call ddadd (s0, s1, b)

100 continue

return
end subroutine

subroutine ddsub (dda, ddb, ddc)

! This subroutine computes ddc = dda - ddb.

implicit none
real*8 dda(2), ddb(2), ddc(2)
real*8 e, t1, t2

! Compute dda + ddb using Knuth's trick.

t1 = dda(1) - ddb(1)
e = t1 - dda(1)
t2 = ((-ddb(1) - e) + (dda(1) - (t1 - e))) + dda(2) - ddb(2)

! The result is t1 + t2, after normalization.

ddc(1) = t1 + t2
ddc(2) = t2 - (ddc(1) - t1)
return
end subroutine

real*8 function dddigin (ca, n)
implicit none
real*8 d1
character*(*), ca
character*16 digits
integer i, k, n
parameter (digits = '0123456789')

d1 = 0.d0

do i = 1, n
k = index (digits, ca(i:i)) - 1
if (k < 0) then
write (6, *) 'dddigin: non-digit in character string'
elseif (k <= 9) then
d1 = 10.d0 * d1 + k
endif
enddo

dddigin = d1
end function

character*16 function dddigout (a, n)
implicit none
real*8 a, d1, d2
character*16 ca, digits
parameter (digits = '0123456789')
integer i, is, k, n

ca = ' '
is = sign (1.d0, a)
d1 = abs (a)

do i = n, 1, -1
d2 = aint (d1 / 10.d0)
k = 1.d0 + (d1 - 10.d0 * d2)
d1 = d2
ca(i:i) = digits(k:k)
if (d1 == 0.d0) goto 100
enddo

i = 0

100 continue

if (is < 0 .and. i > 1) then
ca(i-1:i-1) = '-'
elseif (i == 0 .or. is < 0 .and. i == 1) then
ca = '****************'
endif

dddigout = ca
return
end function

end module

module dddefmod



use ddfunmod
private kdb
parameter (kdb = kind (0.d0))
type dd_real
sequence
real*8 ddr(2)
end type
type dd_complex
sequence
real*8 ddc(4)
end type
type (dd_real), public:: ddeps
real*8, private:: ddt0(4), ddt1(4), ddt2(4), ddt3(4), ddt4(4)

contains

subroutine ddxzc (a, b)

! This converts the DC variable A to the DDC variable B.
! This routine is not intended to be called directly by the user.

complex (kdb) a
real*8 b(4)

b(1) = a
b(2) = 0.d0
b(3) = aimag (a)
b(4) = 0.d0
return
end subroutine

subroutine ddmzc (a, b)



real*8 a(2), b(4)

b(1) = a(1)
b(2) = a(2)
b(3) = 0.d0
b(4) = 0.d0
return
end subroutine

end module

module ddrealmod



use ddfunmod
use dddefmod
private kdb
parameter (kdb = kind (0.d0))
real*8, private:: ddt0(4), ddt1(4), ddt2(4), ddt3(4), ddt4(4)
private &
dd_eqqq, dd_eqqz, dd_eqrq, dd_eqqr, dd_eqdq, dd_eqqd, dd_eqxq, dd_eqqx, &
dd_addqq, dd_addqz, dd_adddq, dd_addqd, dd_addxq, dd_addqx, &
dd_subqq, dd_subqz, dd_subdq, dd_subqd, dd_subxq, dd_subqx, dd_negq, &
dd_mulqq, dd_mulqz, dd_muldq, dd_mulqd, dd_mulxq, dd_mulqx, &
dd_divqq, dd_divqz, dd_divdq, dd_divqd, dd_divxq, dd_divqx, &
dd_expqq, dd_expqi, dd_expdq, dd_expqd, &
dd_eqtqq, dd_eqtqz, dd_eqtdq, dd_eqtqd, dd_eqtxq, dd_eqtqx, &
dd_netqq, dd_netqz, dd_netdq, dd_netqd, dd_netxq, dd_netqx, &
dd_letqq, dd_letdq, dd_letqd, dd_getqq, dd_getdq, dd_getqd, &
dd_lttqq, dd_lttdq, dd_lttqd, dd_gttqq, dd_gttdq, dd_gttqd

! DDR operator extension interface blocks.

interface assignment (=)
module procedure dd_eqqq
module procedure dd_eqqz
module procedure dd_eqdq
module procedure dd_eqqd
module procedure dd_eqxq
module procedure dd_eqqx

module procedure dd_eqqa
end interface

interface operator (+)
module procedure dd_addqq
module procedure dd_addqz
module procedure dd_adddq
module procedure dd_addqd
module procedure dd_addxq
module procedure dd_addqx
end interface

interface operator (-)
module procedure dd_subqq
module procedure dd_subqz
module procedure dd_subdq
module procedure dd_subqd
module procedure dd_subxq
module procedure dd_subqx

module procedure dd_negq
end interface

interface operator (*)
module procedure dd_mulqq
module procedure dd_mulqz
module procedure dd_muldq
module procedure dd_mulqd
module procedure dd_mulxq
module procedure dd_mulqx
end interface

interface operator (/)
module procedure dd_divqq
module procedure dd_divqz
module procedure dd_divdq
module procedure dd_divqd
module procedure dd_divxq
module procedure dd_divqx
end interface

interface operator (**)
module procedure dd_expqq
module procedure dd_expqi
module procedure dd_expdq
module procedure dd_expqd
end interface

interface operator (.eq.)
module procedure dd_eqtqq
module procedure dd_eqtqz
module procedure dd_eqtdq
module procedure dd_eqtqd
module procedure dd_eqtxq
module procedure dd_eqtqx
end interface

interface operator (.ne.)
module procedure dd_netqq
module procedure dd_netqz
module procedure dd_netdq
module procedure dd_netqd
module procedure dd_netxq
module procedure dd_netqx
end interface

interface operator (.le.)
module procedure dd_letqq
module procedure dd_letdq
module procedure dd_letqd
end interface

interface operator (.ge.)
module procedure dd_getqq
module procedure dd_getdq
module procedure dd_getqd
end interface

interface operator (.lt.)
module procedure dd_lttqq
module procedure dd_lttdq
module procedure dd_lttqd
end interface

interface operator (.gt.)
module procedure dd_gttqq
module procedure dd_gttdq
module procedure dd_gttqd
end interface

contains

! DDR assignment routines.

subroutine dd_eqqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: qb
qa%ddr(1) = qb%ddr(1)
qa%ddr(2) = qb%ddr(2)
return
end subroutine

subroutine dd_eqqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: zb
qa%ddr(1) = zb%ddc(1)
qa%ddr(2) = zb%ddc(2)
return
end subroutine

subroutine dd_eqdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: da
intent (in):: qb
da = qb%ddr(1)
return
end subroutine

subroutine dd_eqqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: db
qa%ddr(1) = db
qa%ddr(2) = 0.d0
return
end subroutine

subroutine dd_eqxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: xa
intent (in):: qb
xa = qb%ddr(1)
return
end subroutine

subroutine dd_eqqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: qa
intent (in):: xb
qa%ddr(1) = xb
qa%ddr(2) = 0.d0
return
end subroutine

subroutine dd_eqqa (qa, ab)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
character*(*), intent (in):: ab
intent (out):: qa
character*80 t
t = ab
call ddinpc (t, qa%ddr)
return
end subroutine

! DDR add routines.

function dd_addqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_addqq
intent (in):: qa, qb
call ddadd (qa%ddr, qb%ddr, dd_addqq%ddr)
return
end function

function dd_addqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcadd (ddt1, zb%ddc, dd_addqz%ddc)
return
end function

function dd_adddq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_adddq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddadd (ddt1, qb%ddr, dd_adddq%ddr)
return
end function

function dd_addqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_addqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddadd (qa%ddr, ddt1, dd_addqd%ddr)
return
end function

function dd_addxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcadd (ddt1, ddt2, dd_addxq%ddc)
return
end function

function dd_addqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcadd (ddt1, ddt2, dd_addqx%ddc)
return
end function

! DDR subtract routines.

function dd_subqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_subqq
intent (in):: qa, qb
call ddsub (qa%ddr, qb%ddr, dd_subqq%ddr)
return
end function

function dd_subqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcsub (ddt1, zb%ddc, dd_subqz%ddc)
return
end function

function dd_subdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_subdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddsub (ddt1, qb%ddr, dd_subdq%ddr)
return
end function

function dd_subqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_subqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddsub (qa%ddr, ddt1, dd_subqd%ddr)
return
end function

function dd_subxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcsub (ddt1, ddt2, dd_subxq%ddc)
return
end function

function dd_subqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcsub (ddt1, ddt2, dd_subqx%ddc)
return
end function

! DDR negation routine.

function dd_negq (qa)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_negq
intent (in):: qa
call ddeq (qa%ddr, dd_negq%ddr)
dd_negq%ddr(1) = - qa%ddr(1)
dd_negq%ddr(2) = - qa%ddr(2)
return
end function

! DDR multiply routines.

function dd_mulqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_mulqq
intent (in):: qa, qb
call ddmul (qa%ddr, qb%ddr, dd_mulqq%ddr)
return
end function

function dd_mulqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcmul (ddt1, zb%ddc, dd_mulqz%ddc)
return
end function

function dd_muldq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_muldq
intent (in):: da, qb
call ddmuld (qb%ddr, da, dd_muldq%ddr)
return
end function

function dd_mulqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_mulqd
intent (in):: qa, db
call ddmuld (qa%ddr, db, dd_mulqd%ddr)
return
end function

function dd_mulxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcmul (ddt1, ddt2, dd_mulxq%ddc)
return
end function

function dd_mulqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcmul (ddt1, ddt2, dd_mulqx%ddc)
return
end function

! DDR divide routines.

function dd_divqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_divqq
intent (in):: qa, qb
call dddiv (qa%ddr, qb%ddr, dd_divqq%ddr)
return
end function

function dd_divqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_divqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcdiv (ddt1, zb%ddc, dd_divqz%ddc)
return
end function

function dd_divdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_divdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call dddiv (ddt1, qb%ddr, dd_divdq%ddr)
return
end function

function dd_divqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_divqd
intent (in):: qa, db
call dddivd (qa%ddr, db, dd_divqd%ddr)
return
end function

function dd_divxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_divxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcdiv (ddt1, ddt2, dd_divxq%ddc)
return
end function

function dd_divqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_divqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcdiv (ddt1, ddt2, dd_divqx%ddc)
return
end function

! DDR exponentiation routines.

function dd_expqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expqq
intent (in):: qa, qb
call ddlog (qa%ddr, ddt1)
call ddmul (ddt1, qb%ddr, ddt2)
call ddexp (ddt2, dd_expqq%ddr)
return
end function

function dd_expqi (qa, ib)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expqi
intent (in):: qa, ib
call ddnpwr (qa%ddr, ib, dd_expqi%ddr)
return
end function

function dd_expdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddlog (ddt1, ddt2)
call ddmul (ddt2, qb%ddr, ddt3)
call ddexp (ddt3, dd_expdq%ddr)
return
end function

function dd_expqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_real):: dd_expqd
intent (in):: qa, db
call ddlog (qa%ddr, ddt1)
call ddmuld (ddt1, db, ddt2)
call ddexp (ddt2, dd_expqd%ddr)
return
end function

! DDR .EQ. routines.

function dd_eqtqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .eq. 0) then
dd_eqtqq = .true.
else
dd_eqtqq = .false.
endif
return
end function

function dd_eqtqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcpr (ddt1, zb%ddc, ic1)
call ddcpr (ddt1(3), zb%ddc(3), ic2)
if (ic1 .eq. 0 .and. ic2 .eq. 0) then
dd_eqtqz = .true.
else
dd_eqtqz = .false.
endif
return
end function

function dd_eqtdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .eq. 0) then
dd_eqtdq = .true.
else
dd_eqtdq = .false.
endif
return
end function

function dd_eqtqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .eq. 0) then
dd_eqtqd = .true.
else
dd_eqtqd = .false.
endif
return
end function

function dd_eqtxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .eq. 0 .and. ic2 .eq. 0) then
dd_eqtxq = .true.
else
dd_eqtxq = .false.
endif
return
end function

function dd_eqtqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_eqtqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .eq. 0 .and. ic2 .eq. 0) then
dd_eqtqx = .true.
else
dd_eqtqx = .false.
endif
return
end function

! DDR .NE. routines.

function dd_netqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .ne. 0) then
dd_netqq = .true.
else
dd_netqq = .false.
endif
return
end function

function dd_netqz (qa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqz
intent (in):: qa, zb
call ddmzc (qa%ddr, ddt1)
call ddcpr (ddt1, zb%ddc, ic1)
call ddcpr (ddt1(3), zb%ddc(3), ic2)
if (ic1 .ne. 0 .or. ic2 .ne. 0) then
dd_netqz = .true.
else
dd_netqz = .false.
endif
return
end function

function dd_netdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .ne. 0) then
dd_netdq = .true.
else
dd_netdq = .false.
endif
return
end function

function dd_netqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .ne. 0) then
dd_netqd = .true.
else
dd_netqd = .false.
endif
return
end function

function dd_netxq (xa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netxq
intent (in):: xa, qb
call ddxzc (xa, ddt1)
call ddmzc (qb%ddr, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .ne. 0 .or. ic2 .ne. 0) then
dd_netxq = .true.
else
dd_netxq = .false.
endif
return
end function

function dd_netqx (qa, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_netqx
intent (in):: qa, xb
call ddmzc (qa%ddr, ddt1)
call ddxzc (xb, ddt2)
call ddcpr (ddt1, ddt2, ic1)
call ddcpr (ddt1(3), ddt2(3), ic2)
if (ic1 .ne. 0 .or. ic2 .ne. 0) then
dd_netqx = .true.
else
dd_netqx = .false.
endif
return
end function

! DDR .LE. routines.

function dd_letqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_letqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .le. 0) then
dd_letqq = .true.
else
dd_letqq = .false.
endif
return
end function

function dd_letdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_letdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .le. 0) then
dd_letdq = .true.
else
dd_letdq = .false.
endif
return
end function

function dd_letqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_letqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .le. 0) then
dd_letqd = .true.
else
dd_letqd = .false.
endif
return
end function

! DDR .GE. routines.

function dd_getqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_getqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .ge. 0) then
dd_getqq = .true.
else
dd_getqq = .false.
endif
return
end function

function dd_getdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_getdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .ge. 0) then
dd_getdq = .true.
else
dd_getdq = .false.
endif
return
end function

function dd_getqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_getqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .ge. 0) then
dd_getqd = .true.
else
dd_getqd = .false.
endif
return
end function

! DDR .LT. routines.

function dd_lttqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_lttqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .lt. 0) then
dd_lttqq = .true.
else
dd_lttqq = .false.
endif
return
end function

function dd_lttdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_lttdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .lt. 0) then
dd_lttdq = .true.
else
dd_lttdq = .false.
endif
return
end function

function dd_lttqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_lttqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .lt. 0) then
dd_lttqd = .true.
else
dd_lttqd = .false.
endif
return
end function

! DDR .GT. routines.

function dd_gttqq (qa, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_gttqq
intent (in):: qa, qb
call ddcpr (qa%ddr, qb%ddr, ic)
if (ic .gt. 0) then
dd_gttqq = .true.
else
dd_gttqq = .false.
endif
return
end function

function dd_gttdq (da, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_gttdq
intent (in):: da, qb
ddt1(1) = da
ddt1(2) = 0.d0
call ddcpr (ddt1, qb%ddr, ic)
if (ic .gt. 0) then
dd_gttdq = .true.
else
dd_gttdq = .false.
endif
return
end function

function dd_gttqd (qa, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
logical dd_gttqd
intent (in):: qa, db
ddt1(1) = db
ddt1(2) = 0.d0
call ddcpr (qa%ddr, ddt1, ic)
if (ic .gt. 0) then
dd_gttqd = .true.
else
dd_gttqd = .false.
endif
return
end function

end module


module ddcmpmod


use ddfunmod
use dddefmod
private kdb
parameter (kdb = kind (0.d0))
real*8, private:: ddt0(4), ddt1(4), ddt2(4), ddt3(4), ddt4(4)
private &
dd_eqzq, dd_eqzz, dd_eqdz, dd_eqzd, dd_eqxz, dd_eqzx, &
dd_addzq, dd_addzz, dd_adddz, dd_addzd, dd_addxz, dd_addzx, &
dd_subzq, dd_subzz, dd_subdz, dd_subzd, dd_subxz, dd_subzx, dd_negz, &
dd_mulzq, dd_mulzz, dd_muldz, dd_mulzd, dd_mulxz, dd_mulzx, &
dd_divzq, dd_divzz, dd_divdz, dd_divzd, dd_divxz, dd_divzx, &
dd_expzi, &
dd_eqtzq, dd_eqtzz, dd_eqtdz, dd_eqtzd, dd_eqtxz, dd_eqtzx, &
dd_netzq, dd_netzz, dd_netdz, dd_netzd, dd_netxz, dd_netzx

! DDR operator extension interface blocks.

interface assignment (=)
module procedure dd_eqzq
module procedure dd_eqzz
module procedure dd_eqdz
module procedure dd_eqzd
module procedure dd_eqxz
module procedure dd_eqzx
end interface

interface operator (+)
module procedure dd_addzq
module procedure dd_addzz
module procedure dd_adddz
module procedure dd_addzd
module procedure dd_addxz
module procedure dd_addzx
end interface

interface operator (-)
module procedure dd_subzq
module procedure dd_subzz
module procedure dd_subdz
module procedure dd_subzd
module procedure dd_subxz
module procedure dd_subzx

module procedure dd_negz
end interface

interface operator (*)
module procedure dd_mulzq
module procedure dd_mulzz
module procedure dd_muldz
module procedure dd_mulzd
module procedure dd_mulxz
module procedure dd_mulzx
end interface

interface operator (/)
module procedure dd_divzq
module procedure dd_divzz
module procedure dd_divdz
module procedure dd_divzd
module procedure dd_divxz
module procedure dd_divzx
end interface

interface operator (**)
module procedure dd_expzi
end interface

interface operator (.eq.)
module procedure dd_eqtzq
module procedure dd_eqtzz
module procedure dd_eqtdz
module procedure dd_eqtzd
module procedure dd_eqtxz
module procedure dd_eqtzx
end interface

interface operator (.ne.)
module procedure dd_netzq
module procedure dd_netzz
module procedure dd_netdz
module procedure dd_netzd
module procedure dd_netxz
module procedure dd_netzx
end interface

contains

! DDC assignment routines.

subroutine dd_eqzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: qb
call ddmzc (qb%ddr, za%ddc)
return
end subroutine

subroutine dd_eqzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: zb
call ddceq (zb%ddc, za%ddc)
return
end subroutine

subroutine dd_eqdz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: da
intent (in):: zb
da = zb%ddc(1)
return
end subroutine

subroutine dd_eqzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: db
xb = db
call ddxzc (xb, za%ddc)
return
end subroutine

subroutine dd_eqxz (xa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: xa
intent (in):: zb
db = zb%ddc(1)
dc = zb%ddc(3)
xa = cmplx (db, dc, kdb)
return
end subroutine

subroutine dd_eqzx (za, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
intent (out):: za
intent (in):: xb
call ddxzc (xb, za%ddc)
return
end subroutine

! DDC add routines.

function dd_addzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzq
intent (in):: za, qb
call ddmzc (qb%ddr, ddt1)
call ddcadd (za%ddc, ddt1, dd_addzq%ddc)
return
end function

function dd_addzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzz
intent (in):: za, zb
call ddcadd (za%ddc, zb%ddc, dd_addzz%ddc)
return
end function

function dd_adddz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_adddz
intent (in):: da, zb
xa = da
call ddxzc (xa, ddt1)
call ddcadd (ddt1, zb%ddc, dd_adddz%ddc)
return
end function

function dd_addzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzd
intent (in):: za, db
xb = db
call ddxzc (xb, ddt1)
call ddcadd (za%ddc, ddt1, dd_addzd%ddc)
return
end function

function dd_addxz (xa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addxz
intent (in):: xa, zb
call ddxzc (xa, ddt1)
call ddcadd (ddt1, zb%ddc, dd_addxz%ddc)
return
end function

function dd_addzx (za, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_addzx
intent (in):: za, xb
call ddxzc (xb, ddt2)
call ddcadd (za%ddc, ddt2, dd_addzx%ddc)
return
end function

! DDC subtract routines.

function dd_subzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzq
intent (in):: za, qb
call ddmzc (qb%ddr, ddt1)
call ddcsub (za%ddc, ddt1, dd_subzq%ddc)
return
end function

function dd_subzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzz
intent (in):: za, zb
call ddcsub (za%ddc, zb%ddc, dd_subzz%ddc)
return
end function

function dd_subdz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subdz
intent (in):: da, zb
xa = da
call ddxzc (xa, ddt1)
call ddcsub (ddt1, zb%ddc, dd_subdz%ddc)
return
end function

function dd_subzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzd
intent (in):: za, db
xb = db
call ddxzc (xb, ddt1)
call ddcsub (za%ddc, ddt1, dd_subzd%ddc)
return
end function

function dd_subxz (xa, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subxz
intent (in):: xa, zb
call ddxzc (xa, ddt1)
call ddcsub (ddt1, zb%ddc, dd_subxz%ddc)
return
end function

function dd_subzx (za, xb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_subzx
intent (in):: za, xb
call ddxzc (xb, ddt2)
call ddcsub (za%ddc, ddt2, dd_subzx%ddc)
return
end function

! DDC negation routine.

function dd_negz (za)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_negz
intent (in):: za
call ddceq (za%ddc, dd_negz%ddc)
dd_negz%ddc(1) = - za%ddc(1)
dd_negz%ddc(2) = - za%ddc(2)
dd_negz%ddc(3) = - za%ddc(3)
dd_negz%ddc(4) = - za%ddc(4)
return
end function

! DDC multiply routines.

function dd_mulzq (za, qb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulzq
intent (in):: za, qb
call ddmzc (qb%ddr, ddt1)
call ddcmul (za%ddc, ddt1, dd_mulzq%ddc)
return
end function

function dd_mulzz (za, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulzz
intent (in):: za, zb
call ddcmul (za%ddc, zb%ddc, dd_mulzz%ddc)
return
end function

function dd_muldz (da, zb)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_muldz
intent (in):: da, zb
xa = da
call ddxzc (xa, ddt1)
call ddcmul (ddt1, zb%ddc, dd_muldz%ddc)
return
end function

function dd_mulzd (za, db)
implicit real*8 (d), &
type (dd_real) (q), complex (kdb) (x), type (dd_complex) (z)
type (dd_complex):: dd_mulzd
intent (in):: za, db
xb = db
call ddxzc (xb, ddt1)
call ddcmul (za%ddc,
Back to top
View user's profile
mkcolg



Joined: 30 Jun 2004
Posts: 5871
Location: The Portland Group Inc.

PostPosted: Wed Mar 17, 2010 4:32 pm    Post subject: Reply with quote

Hi SWL_EGGBABY,

It looks like the code was too long and got cut off.

Can you please send the code via email to PGI Customer Service (trs@pgroup.com)?

Thanks,
Mat
Back to top
View user's profile
mkcolg



Joined: 30 Jun 2004
Posts: 5871
Location: The Portland Group Inc.

PostPosted: Thu Mar 18, 2010 2:22 pm    Post subject: Reply with quote

Hi SWL_EGGBABY,

Customer service sent me your code and it does look like a compiler error due to your derived types. I've created a report (TPR#16733) and sent it to our engineers for further investigation.

As a work around, please replace "type(dd_real)" with "real*8" in your matmul code.

Thanks,
Mat
Back to top
View user's profile
SWL_EGGBABY



Joined: 16 Dec 2009
Posts: 31

PostPosted: Thu Mar 18, 2010 11:06 pm    Post subject: Compute four precision on GPU but Permission denied why? Reply with quote

Hi Mat,

real*8 is double precision.
How can I compute four precision on GPU. The ddfunmod can use type(dd_real) compute four precision on CPU,the result can retain the 32-bit effective number.I need it.

Thanks,
Back to top
View user's profile
mkcolg



Joined: 30 Jun 2004
Posts: 5871
Location: The Portland Group Inc.

PostPosted: Mon Mar 29, 2010 1:24 pm    Post subject: Reply with quote

Hi SWL_EGGBABY,

GPUs don't support quad precision. You might be able to re-write the matmul kernel to simulate quad precision. If you're successful, let us know and we can find a place where we can post it since it may be useful for others as well.

- Mat
Back to top
View user's profile
Display posts from previous:   
Post new topic   Reply to topic    PGI User Forum Forum Index -> Accelerator Programming All times are GMT - 7 Hours
Goto page Previous  1, 2, 3, 4  Next
Page 3 of 4

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © phpBB Group