PGI User Forum
 Search   Memberlist       Profile    Log in

 Compute four precision on GPU but Permission denied why? Goto page Previous  1, 2, 3, 4  Next
Author Message
SWL_EGGBABY

Joined: 16 Dec 2009
Posts: 31

 Posted: Wed Mar 17, 2010 12:24 am    Post subject: Compute four precision on GPU but Permission denied why? 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,
mkcolg

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

 Posted: Wed Mar 17, 2010 4:32 pm    Post subject: 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
mkcolg

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

 Posted: Thu Mar 18, 2010 2:22 pm    Post subject: 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
SWL_EGGBABY

Joined: 16 Dec 2009
Posts: 31

 Posted: Thu Mar 18, 2010 11:06 pm    Post subject: Compute four precision on GPU but Permission denied why? 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,
mkcolg

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

 Posted: Mon Mar 29, 2010 1:24 pm    Post subject: 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
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT - 7 HoursGoto page Previous  1, 2, 3, 4  Next Page 3 of 4

 Jump to: Select a forum General Information----------------New Release Announcements User Forums----------------Programming and CompilingAccelerator ProgrammingDebugging and ProfilingLicenses and Installation
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