Last active
September 1, 2021 19:22
-
-
Save weslleyspereira/ac5babcc0bc370cb238012e3048440e3 to your computer and use it in GitHub Desktop.
Test double complex division in Fortran
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| *> \brief zdiv tests the robustness and precision of the double complex division | |
| *> \author Weslley S. Pereira, University of Colorado Denver, U.S. | |
| * | |
| *> \verbatim | |
| *> | |
| *> Real values for test: | |
| *> (1) x = 2**m, where m = MINEXPONENT-DIGITS, ..., MINEXPONENT-1. Stop on the first success. | |
| *> Mind that not all platforms might implement subnormal numbers. | |
| *> (2) x = 2**m, where m = MINEXPONENT, ..., 0. Stop on the first success. | |
| *> (3) x = OV, where OV is the overflow threshold. OV^2 overflows but the norm is OV. | |
| *> (4) x = 2**m, where m = MAXEXPONENT-1, ..., 1. Stop on the first success. | |
| *> | |
| *> Tests: | |
| *> (a) y = x + 0 * I, y/y = 1 | |
| *> (b) y = 0 + x * I, y/y = 1 | |
| *> (c) y = x + x * I, y/y = 1 | |
| *> (d) y1 = 0 + x * I, y2 = x + 0 * I, y1/y2 = I | |
| *> (e) y1 = 0 + x * I, y2 = x + 0 * I, y2/y1 = -I | |
| *> (f) y = x + x * I, y/conj(y) = I | |
| *> | |
| *> Special cases: | |
| *> | |
| *> (i) Inf inputs: | |
| *> (1) y = ( Inf + 0 * I) | |
| *> (2) y = ( 0 + Inf * I) | |
| *> (3) y = (-Inf + 0 * I) | |
| *> (4) y = ( 0 - Inf * I) | |
| *> (5) y = ( Inf + Inf * I) | |
| *> Tests: | |
| *> (a) 0 / y is either 0 or NaN. | |
| *> (b) 1 / y is either 0 or NaN. | |
| *> (c) y / y is NaN. | |
| *> | |
| *> (n) NaN inputs: | |
| *> (1) y = (NaN + 0 * I) | |
| *> (2) y = (0 + NaN * I) | |
| *> (3) y = (NaN + NaN * I) | |
| *> Tests: | |
| *> (a) 0 / y is NaN. | |
| *> (b) 1 / y is NaN. | |
| *> (c) y / y is NaN. | |
| *> | |
| *> \endverbatim | |
| * | |
| program zdiv | |
| integer N, i, nNaN, nInf, min, Max, m | |
| parameter ( N = 4, nNaN = 3, nInf = 5 ) | |
| double precision X( N ), threeFourth, fiveFourth, aInf, aNaN, b, | |
| $ eps, blueMin, blueMax, OV, Xj, stepX(N), limX(N) | |
| parameter ( threeFourth = 3.0d0 / 4, | |
| $ fiveFourth = 5.0d0 / 4 ) | |
| double complex Y, Y2, R, cInf( nInf ), cNaN( nNaN ), czero, | |
| $ cone | |
| parameter ( czero = DCMPLX( 0.0d0, 0.0d0 ), | |
| $ cone = DCMPLX( 1.0d0, 0.0d0 ) ) | |
| * | |
| intrinsic DCONJG, DBLE, RADIX, CEILING, TINY, DIGITS, | |
| $ MAXEXPONENT, MINEXPONENT, FLOOR, HUGE, DCMPLX, | |
| $ EPSILON | |
| * | |
| min = MINEXPONENT(0.0d0) | |
| Max = MAXEXPONENT(0.0d0) | |
| m = DIGITS(0.0d0) | |
| b = DBLE(RADIX(0.0d0)) | |
| eps = EPSILON(0.0d0) | |
| blueMin = b**CEILING( (min - 1) * 0.5d0 ) | |
| blueMax = b**FLOOR( (Max - m + 1) * 0.5d0 ) | |
| OV = HUGE(0.0d0) | |
| * | |
| X(1) = TINY(0.0d0) * b**( DBLE(1-m) ) | |
| X(2) = TINY(0.0d0) | |
| X(3) = OV | |
| X(4) = b**( DBLE(Max-1) ) | |
| * | |
| stepX(1) = 2.0 | |
| stepX(2) = 2.0 | |
| stepX(3) = 0.0 | |
| stepX(4) = 0.5 | |
| * | |
| limX(1) = X(2) | |
| limX(2) = 1.0 | |
| limX(3) = 0.0 | |
| limX(4) = 2.0 | |
| * | |
| print *, '# X :=', X | |
| print *, '# Blue min constant :=', blueMin | |
| print *, '# Blue max constant :=', blueMax | |
| * | |
| Xj = X(1) | |
| if( Xj .eq. 0.0d0 ) then | |
| print *, "# Subnormal numbers treated as 0" | |
| else | |
| do 100 i = 1, N | |
| Xj = X(i) | |
| if( Xj .eq. 0.0d0 ) print *, | |
| $ "# Subnormal numbers may be treated as 0" | |
| 100 continue | |
| endif | |
| * | |
| aInf = OV * 2 | |
| cInf(1) = DCMPLX( aInf, 0.0d0 ) | |
| cInf(2) = DCMPLX(-aInf, 0.0d0 ) | |
| cInf(3) = DCMPLX( 0.0d0, aInf ) | |
| cInf(4) = DCMPLX( 0.0d0,-aInf ) | |
| cInf(5) = DCMPLX( aInf, aInf ) | |
| * | |
| aNaN = aInf / aInf | |
| cNaN(1) = DCMPLX( aNaN, 0.0d0 ) | |
| cNaN(2) = DCMPLX( 0.0d0, aNaN ) | |
| cNaN(3) = DCMPLX( aNaN, aNaN ) | |
| * | |
| * Test (a) y = x + 0 * I, y/y = 1 | |
| do 10 i = 1, N | |
| Xj = X(i) | |
| if( Xj .eq. 0.0d0 ) then | |
| print *, "# [a] Subnormal numbers may be treated as 0" | |
| else | |
| do while( Xj .ne. limX(i) ) | |
| Y = DCMPLX( Xj, 0.0d0 ) | |
| R = Y / Y | |
| if( R .ne. 1.0D0 ) | |
| $ WRITE( *, FMT = 9999 ) 'a',i, Y, Y, R, 1.0D0 | |
| Xj = Xj * stepX(i) | |
| end do | |
| endif | |
| 10 continue | |
| * | |
| * Test (b) y = 0 + x * I, y/y = 1 | |
| do 20 i = 1, N | |
| Xj = X(i) | |
| if( Xj .eq. 0.0d0 ) then | |
| print *, "# [b] Subnormal numbers may be treated as 0" | |
| else | |
| do while( Xj .ne. limX(i) ) | |
| Y = DCMPLX( 0.0d0, Xj ) | |
| R = Y / Y | |
| if( R .ne. 1.0D0 ) | |
| $ WRITE( *, FMT = 9999 ) 'b',i, Y, Y, R, 1.0D0 | |
| Xj = Xj * stepX(i) | |
| end do | |
| endif | |
| 20 continue | |
| * | |
| * Test (c) y = x + x * I, y/y = 1 | |
| do 30 i = 1, N | |
| Xj = X(i) | |
| if( Xj .eq. 0.0d0 ) then | |
| print *, "# [c] Subnormal numbers may be treated as 0" | |
| else | |
| do while( Xj .ne. limX(i) ) | |
| Y = DCMPLX( Xj, Xj ) | |
| R = Y / Y | |
| if( R .ne. 1.0D0 ) | |
| $ WRITE( *, FMT = 9999 ) 'c',i, Y, Y, R, 1.0D0 | |
| Xj = Xj * stepX(i) | |
| end do | |
| endif | |
| 30 continue | |
| * | |
| * Test (d) y1 = 0 + x * I, y2 = x + 0 * I, y1/y2 = I | |
| do 40 i = 1, N | |
| Xj = X(i) | |
| if( Xj .eq. 0.0d0 ) then | |
| print *, "# [d] Subnormal numbers may be treated as 0" | |
| else | |
| do while( Xj .ne. limX(i) ) | |
| Y = DCMPLX( 0.0d0, Xj ) | |
| Y2 = DCMPLX( Xj, 0.0d0 ) | |
| R = Y / Y2 | |
| if( R .ne. DCMPLX(0.0D0,1.0D0) ) | |
| $ WRITE( *, FMT = 9999 ) 'd',i, Y, Y2, R, | |
| $ DCMPLX(0.0D0,1.0D0) | |
| Xj = Xj * stepX(i) | |
| end do | |
| endif | |
| 40 continue | |
| * | |
| * Test (e) y1 = 0 + x * I, y2 = x + 0 * I, y2/y1 = -I | |
| do 50 i = 1, N | |
| Xj = X(i) | |
| if( Xj .eq. 0.0d0 ) then | |
| print *, "# [e] Subnormal numbers may be treated as 0" | |
| else | |
| do while( Xj .ne. limX(i) ) | |
| Y = DCMPLX( 0.0d0, Xj ) | |
| Y2 = DCMPLX( Xj, 0.0d0 ) | |
| R = Y2 / Y | |
| if( R .ne. DCMPLX(0.0D0,-1.0D0) ) | |
| $ WRITE( *, FMT = 9999 ) 'e',i, Y2, Y, R, | |
| $ DCMPLX(0.0D0,-1.0D0) | |
| Xj = Xj * stepX(i) | |
| end do | |
| endif | |
| 50 continue | |
| * | |
| * Test (f) y = x + x * I, y/conj(y) = I | |
| do 60 i = 1, N | |
| Xj = X(i) | |
| if( Xj .eq. 0.0d0 ) then | |
| print *, "# [f] Subnormal numbers may be treated as 0" | |
| else | |
| do while( Xj .ne. limX(i) ) | |
| Y = DCMPLX( Xj, Xj ) | |
| R = Y / DCONJG( Y ) | |
| if( R .ne. DCMPLX(0.0D0,1.0D0) ) | |
| $ WRITE( *, FMT = 9999 ) 'f',i, Y, DCONJG( Y ), R, | |
| $ DCMPLX(0.0D0,1.0D0) | |
| Xj = Xj * stepX(i) | |
| end do | |
| endif | |
| 60 continue | |
| * | |
| * Test (g) Infs | |
| do 70 i = 1, nInf | |
| Y = cInf(i) | |
| R = czero / Y | |
| if( (R .ne. czero) .and. (R .eq. R) ) then | |
| WRITE( *, FMT = 9998 ) 'ia',i, czero, Y, R, 'NaN and 0' | |
| endif | |
| R = cone / Y | |
| if( (R .ne. czero) .and. (R .eq. R) ) then | |
| WRITE( *, FMT = 9998 ) 'ib',i, cone, Y, R, 'NaN and 0' | |
| endif | |
| R = Y / Y | |
| if( R .eq. R ) then | |
| WRITE( *, FMT = 9998 ) 'ic',i, Y, Y, R, 'NaN' | |
| endif | |
| 70 continue | |
| * | |
| * Test (h) NaNs | |
| do 80 i = 1, nNaN | |
| Y = cNaN(i) | |
| R = czero / Y | |
| if( R .eq. R ) then | |
| WRITE( *, FMT = 9998 ) 'na',i, czero, Y, R, 'NaN' | |
| endif | |
| R = cone / Y | |
| if( R .eq. R ) then | |
| WRITE( *, FMT = 9998 ) 'nb',i, cone, Y, R, 'NaN' | |
| endif | |
| R = Y / Y | |
| if( R .eq. R ) then | |
| WRITE( *, FMT = 9998 ) 'nc',i, Y, Y, R, 'NaN' | |
| endif | |
| 80 continue | |
| * | |
| 9998 FORMAT( '[',A2,I1, '] (', (ES10.3,SP,ES10.3,"*I"), ' ) / ( ', | |
| $ (ES10.3,SP,ES10.3,"*I"), ' ) = ', | |
| $ (ES10.3,SP,ES10.3,"*I"), ' differs from ', A10 ) | |
| * | |
| 9999 FORMAT( '[',A2,I1, '] (', (ES10.3,SP,ES10.3,"*I"), ' ) / ( ', | |
| $ (ES10.3,SP,ES10.3,"*I"), ' ) = ', | |
| $ (ES10.3,SP,ES10.3,"*I"), ' differs from ', | |
| $ (ES10.3,SP,ES10.3,"*I") ) | |
| * | |
| * End of zdiv | |
| * | |
| END |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment