50 likes | 249 Views
Double precision, real and complex. To make a single precision real number, we declare REAL :: X, Y, Z These are accurate to 7 significant digits, usually stored in 32 bits. Single precision is the default If we want double precision, REAL*8 :: X, Y, Z DOUBLE PRECISION :: A, B, C
E N D
Double precision, real and complex To make a single precision real number, we declare REAL :: X, Y, Z These are accurate to 7 significant digits, usually stored in 32 bits. Single precision is the default If we want double precision, REAL*8 :: X, Y, Z DOUBLE PRECISION :: A, B, C INTEGER, PARAMETER :: Prec14 = SELECTED_REAL_KIND(14) REAL(KIND = Prec14) :: D, E, F Each of these gives 14 digit precision in 64 bits
Declaring values within equations For numbers in an equation, we use for single/double precision: 4.2d0 (4.2 to 14 digits of precision) 4.2e0 (4.2 to 7 digits of precision) Why does it matter? When mixing precisions, if something is defined as single precision, the undetermined digits can be filled with garbage numbers
Complex variables and arrays Again default is single precision COMPLEX :: psi Here, psi will be a complex number with the real and imaginary part precise to 7 digits (64 bits total!) Interestingly, COMPLEX*8 :: psi Also gives single precision
Current project: Double precision complex variables and arrays INTEGER, PARAMETER :: Prec14=SELECTED_REAL_KIND(14) INTEGER :: j,n INTEGER, PARAMETER :: jmax= ?? COMPLEX(KIND=Prec14), DIMENSION(0:jmax-1) :: a,b,c,chi COMPLEX(KIND=Prec14), DIMENSION(0:jmax-1) :: psi COMPLEX(KIND=Prec14) :: ai, ar ! Define i =sqrt(-1) and 1 as complex numbers ai, ar ai=(0.0d0,1.0d0) ar=(1.0d0,0.0d0) ! Define matrix elements of Q do j=0, jmax-1 a(j) = -ai*alpha/8.0d0 b(j) = ar*0.5d0+ai*alpha/4.0d0 c(j) = -ai*alpha/8.0d0 enddo a(0)=0.0d0 ; c(jmax-1)=0.0d0 ! Not necessary … ! Need a loop on time steps n call tridag(a,b,c,psi,chi,jmax) psi=chi-psi
In the student directory, code to invert a tridiagonal matrix SUBROUTINE TRIDAG(A,B,C,R,U,N) INTEGER :: N,J INTEGER, PARAMETER :: NMAX=100 INTEGER, PARAMETER :: DP= SELECTED_REAL_KIND(14) COMPLEX(KIND=DP), BET COMPLEX(KIND=DP), DIMENSION(NMAX) :: GAM COMPLEX(KIND=DP), DIMENSION(N) :: A,B,C,R,U IF(B(1).EQ.0.)PAUSE BET=B(1) U(1)=R(1)/BET DO J=2,N GAM(J)=C(J-1)/BET BET=B(J)-A(J)*GAM(J) IF(BET.EQ.0.)PAUSE U(J)=(R(J)-A(J)*U(J-1))/BET ENDDO DO J=N-1,1,-1 U(J)=U(J)-GAM(J+1)*U(J+1) ENDDO RETURN END