I might need a little hand-holding. How do you see this adapting to command prompts? It seems the actual variable must be encoded in the keyvars constructor, not just a string containing the variable name.
How would you construct a keyvar object at runtime that points to the correct variable based on a user-provided string?
The advantage of using user-define types in all the modules to gather the variables like a common block probably previously did is that you just add a variable to the user-defined object and the code does not have to change in the other routines, although I would imagine routines referencing it would need recompiled/reloaded. With well-chosen names the input can be even more descriptive. The disadvantage being you have to type more verbose input like color%background=“green” in some cases; but using NAMELIST groups allows you to dump the state of all the variables easily, read values from a file; etc.
So to simplify the previous example, instead of common block A and B I create user-defined variables PROP and ENV. If I want a new PROP variable I just add it where I define the type T_PROP and the other code does not change. So if I just added “prop%mass” I can just run the program and enter
prop%mass=321.56
without changing code anywhere else. In code where I use prop%mass it is more clearly some imported global value because it is not just “mass” but “prop%mass”; very easy to add new values; just a bit verbose. It is easy to edit the input to change some other character to % and to add an equal after the first word if you wanted, to
prop_mass 321.56
would work instead of native NAMELIST input; but using the NAMELIST syntax means you can use a language-defined standard syntax – just describe NAMELIST input and your documentation is done.
Just curious what was not appealing about that approach, as it technically met your desired features I thought.
! this file defines the variables in a group of variables defining properties
module a
implicit none
private
type t_prop
real :: mass=0.0
real :: volume=0.0
end type t_prop
public t_prop
end module a
! this file defines the variables in a group of variables defining the environment
module b
implicit none
private
type t_env
real :: air_temp=0.0
real :: pressure=0.0
end type t_env
public t_env
end module b
! this file aggregates the various user-defined variables into a NAMELIST and
! defines a procedure to change values.
module command_line_module
use a, only: t_prop
use b, only: t_env
implicit none
private
type(t_prop),public :: prop; namelist /globals/ prop
type(t_env),public :: env; namelist /globals/ env
public :: read_user_input, globals
contains
function read_user_input() result(status)
logical :: status
character(len=256) :: line
character(len=256) :: answer
integer :: iostat
integer :: lun
character(len=:), allocatable :: intmp
character(len=256) :: iomsg
write (*, '(a)', advance='no') 'read_mode>>'
read (*, '(a)') line
if(line=='stop')then
status=.false.
else
status=.true.
intmp = '&globals '//trim(line)//'/'
read (intmp, nml=globals, iostat=iostat, iomsg=iomsg)
if (iostat .ne. 0) then
write (*, *) 'ERROR:', trim(iomsg)
write (*, *) 'VALUES ARE:'
write (*,nml=globals )
end if
endif
end function read_user_input
end module command_line_module
program main
use command_line_module, only: read_user_input, globals, prop, env
implicit none
character(len=:), allocatable :: status
write(*,'(''>Enter "stop" to end'')')
do while (read_user_input()) ! interactively change NAMELIST group
if (status .eq. 'stop') exit
call print_values
end do
contains
subroutine print_values
write (*, '(a)')'PRINT NAMELIST:'
write (*, globals)
write (*, '(a)')'PRINT USER-DEFINED TYPES:'
write (*,*)'PROPERTIES:',prop
write (*,*)'ENVIRONMENT:',env
write (*, '(a)')'PRINT VALUES IN USER-DEFINED TYPES:'
write (*,*)Prop%MASS, Prop%VOLUME, Env%AIR_TEMP, Env%PRESSURE
end subroutine print_values
end program main
$ ./a.out
Use gfortran interactive extension “?” to print namelist group
>Enter "stop" to end
read_mode>>?
>PRINT NAMELIST:
>&GLOBALS
> PROP%MASS= 0.00000000 ,
> PROP%VOLUME= 0.00000000 ,
> ENV%AIR_TEMP= 0.00000000 ,
> ENV%PRESSURE= 0.00000000 ,
> /
>PRINT USER-DEFINED TYPES:
> PROPERTIES: 0.00000000 0.00000000
> ENVIRONMENT: 0.00000000 0.00000000
>PRINT VALUES IN USER-DEFINED TYPES:
> 0.00000000 0.00000000 0.00000000 0.00000000
Change some values
read_mode>>prop%mass=321,env%pressure=200
>PRINT NAMELIST:
>&GLOBALS
> PROP%MASS= 321.000000 ,
> PROP%VOLUME= 0.00000000 ,
> ENV%AIR_TEMP= 0.00000000 ,
> ENV%PRESSURE= 200.000000 ,
> /
>PRINT USER-DEFINED TYPES:
> PROPERTIES: 321.000000 0.00000000
> ENVIRONMENT: 0.00000000 200.000000
>PRINT VALUES IN USER-DEFINED TYPES:
> 321.000000 0.00000000 0.00000000 200.000000
exit program
read_mode>>stop
Was this question directed at me as the OP? Unclear because you posted in reply to RonShepard, but you seem to touch on some items from my original post.
You as OP
(@Machalot This is just to explore various possibilities (mainly for coding experiment), so please do not care very much about my posts…)
Hi @urbanjost , as for the reasons (or “hesitation”) why people might not prefer derived types as “parameter namespace”, it may be because (i) if existing codes are not based on derived types, a significant amount of code modification may become necessary; and (ii) to reuse old codes as directly as possible, one could use ‘associate’ to extract type components from parameter objects. But with the current functionality of associate, it may be rather tedious to extract many components in every necessary routine (for comparison, some languages seem to provide capability to extract all type components at once (e.g. with in D or Pascal (?)).
RE the tediousness of adding type prefix, it might be useful to define alternative pointer variables in parameter modules (so that the parameters could be read in via derived types, while their components could be imported directly into legacy routines via USE).
module param1_mod
type param1_t
integer :: num = 0
real :: foo = 0
endtype
type(param1_t), target :: par1
!! To allow direct access of contents without "par1".
integer, pointer :: num => par1 % num
real, pointer :: foo => par1 % foo
end
module param2_mod
type param2_t
integer :: num = 0
real :: baa = 0
endtype
type(param2_t), target :: par2
integer, pointer :: num => par2 % num
real, pointer :: baa => par2 % baa
!! Can pre-define alt names for legacy codes (if USE + renaming is tedious).
integer, pointer :: n => par2 % num
real, pointer :: b => par2 % baa
end
program main
use param1_mod, only: par1
use param2_mod, only: par2
implicit none
namelist /globals/ par1, par2
open(10, file="test.inp", status="old")
read(10, globals)
close(10)
block
use param1_mod
print *, "param1: num = ", num, " foo = ", foo
endblock
block
use param2_mod
print *, "param2: num = ", num, " baa = ", baa
print *, "param2: n = ", n, " b = ", b
endblock
end
!! test.inp
&globals
par1%num = 100
par1%foo = 1.1
par2%num = 200
par2%baa = 2.2
/
$ gfortran test.f90 && ./a.out
param1: num = 100 foo = 1.10000002
param2: num = 200 baa = 2.20000005
param2: n = 200 b = 2.20000005
That nicely minimizes or eliminates almost all the changes to the old code while retaining t the simple locality of the variable declarations to my mind.
So a modernization using no external dependencies allowing for some additional interactive commands that allows for the values to be accessed with their original names as well as an aggregated NAMELIST group and grouped via user-defined types where only one file needs changed to add a variable, applied to the originally posted code would be …
module topical_module_1
implicit none
private
type param1_t
real, public :: x1, y1, z1
end type param1_t
type(param1_t), target :: par1
real, public, pointer :: x1 => par1%x1, y1 => par1%y1, z1 => par1%z1
public par1
end module topical_module_1
module topical_module_2
implicit none
private
type param2_t
real, public :: x2, y2, z2
end type param2_t
type(param2_t), target :: par2
real, public, pointer :: x2 => par2%x2, y2 => par2%y2, z2 => par2%z2
public par2
end module topical_module_2
module command_line_module
use topical_module_1, only: par1, x1, y1, z1 ! as user-defined type and "old" names
use topical_module_2, only: par2, x2, y2, z2 ! as user-defined type and "old" names
implicit none
private
public :: read_user_input
public :: print_values
namelist /globals/ par1, par2
contains
subroutine read_user_input(status)
character(len=:), intent(out), allocatable :: status
character(len=256) :: line
character(len=256) :: answer
integer :: ios
integer :: lun
status = ''
write (*, '(a)') 'Enter "help" to list available commands'
do
write (*, '(a)', advance='no') 'read_mode>>'
read (*, '(a)') line
if (line(1:1) == '!') cycle
select case (line)
case ('help')
write (*, '(a)') [character(len=80) :: &
'READ MODE (examine and change values): ', &
' "." return to main program ', &
' NAME=VALUE change values ', &
' show show current globals values ', &
' read read globals from a NAMELIST file ', &
' write write globals to a NAMELIST file ', &
' sh start subshell. use "exit" to return.', &
' stop tell caller to stop ', &
'']
case ('.')
exit
case ('show')
write (*, *) 'SO FAR'
write (*, globals)
case ('stop', 'quit', 'adios', 'adieu', 'arrivederci', 'au revior', 'so long', 'sayonara', 'auf wiedersehen', 'cheerio')
status = 'stop'
exit
case ('sh')
call execute_command_line('bash')
case ('read')
write (*, '(a)', advance='no') 'filename:'
read (*, '(a)', iostat=ios) answer
if (ios /= 0) exit
open (file=answer, iostat=ios, newunit=lun)
if (ios /= 0) exit
read (lun, globals, iostat=ios)
close (unit=lun, iostat=ios)
case ('write')
write (*, '(a)', advance='no') 'filename:'
read (*, '(a)', iostat=ios) answer
if (ios /= 0) exit
open (file=answer, iostat=ios, newunit=lun)
if (ios /= 0) exit
write (lun, globals, iostat=ios)
close (unit=lun, iostat=ios)
case default
UPDATE: block
character(len=:), allocatable :: intmp
character(len=256) :: message
integer :: ios
intmp = '&globals '//trim(line)//'/'
read (intmp, nml=globals, iostat=ios, iomsg=message)
if (ios /= 0) then
write (*, *) 'ERROR:', trim(message)
end if
end block UPDATE
end select
end do
end subroutine read_user_input
subroutine print_values
use topical_module_1, only: x1
use topical_module_2, only: x2
implicit none
write (*, *) 'x1 = ', x1
write (*, *) 'x2 = ', x2
end subroutine print_values
end module command_line_module
program main
! main program
use command_line_module, only: read_user_input
use command_line_module, only: print_values
implicit none
character(len=:), allocatable :: status
call print_values()
do
call read_user_input(status)
if (status == 'stop') stop
call print_values()
end do
end program main
@urbanjost Mainly I haven’t had time to really digest some of the more detailed suggestions on this thread, including yours. I would say @septc has hit on my first hesitation though – I think it would take a huge amount of work to implement derived types as originally proposed. But what the two of you have worked out in the last couple of posts seems very interesting and might actually be workable.
I have no experience using Fortran pointers. Would this require re-declaring all the existing variables as pointers when they are migrated into the topical modules? Can pointers be passed as actual arguments to procedures that correspond to non-pointer dummy arguments?
I am not the ultimate owner of the software in question, I am just trying to provide the team good advice for modernization.
Yes, I am afraid so… The module definition becomes longer than that the namelist approach because we need to declare all variable information again (like types, dimensions, etc). Another potential drawback is that the extensive use of pointers might prohibit some code optimization (because the target of pointers can change at run time, though not very sure about performance).
Yes, I think this is no problem. Once pointers are passed to non-pointer dummy arguments, I think they are treated in the same way as non-pointer variables.
That said, I also feel that the definition of parameter types + additional pointers are rather redundant with the functionality of modules (if those variables are treated as global variables). So overall, if I were to modernize such codes (e.g. by changing all common variables to module variables), I may start from a rather simple approach like the following (without derived types, just manually looping over different namelists in different modules). Also, if possible, I think it would be much simpler to read each namelist group separately from the input file, as usual (rather than merging all lines as in test.inp below).
module param1_mod
integer :: num = 0
real :: foo = 0
namelist /par1/ num, foo
end
module param2_mod
integer :: num = 0
real :: arr(3) = 0
namelist /par2/ num, arr
end
module io_mod
use param1_mod, only: par1
use param2_mod, only: par2
implicit none
contains
subroutine readline(done)
logical, intent(out) :: done
character(256) :: line
character(:), allocatable :: inp_tag, inp_var, buf
integer :: pos
read(*, '(a)') line
print '(/,a)', '>> ' // trim(line)
done = merge(.true., .false., line == "quit")
if (done) return
pos = index(line, ':')
inp_tag = trim(adjustl(line(:pos-1)))
inp_var = trim(adjustl(line(pos+1:)))
select case (inp_tag)
case ('par1')
buf = '&par1 ' // inp_var // ' /'
read(buf, par1)
case ('par2')
buf = '&par2 ' // inp_var // ' /'
read(buf, par2)
case default
print *, 'no match: ignoring input'
return
end select
if (.not. done) then
print par1
print par2
endif
end
end module
program main
use io_mod
logical :: done
done = .false.
do while (.not. done)
call readline(done)
end do
end
!! test.inp
par1 : num = 100
par2 : arr(2:3) = 2.2 3.3
xyz = 777 !! invalid line
quit
$ gfortran test.f90 && ./a.out < test.inp
>> par1 : num = 100
&PAR1
NUM=100 ,
FOO= 0.00000000 ,
/
&PAR2
NUM=0 ,
ARR= 3*0.00000000 ,
/
>> par2 : arr(2:3) = 2.2 3.3
&PAR1
NUM=100 ,
FOO= 0.00000000 ,
/
&PAR2
NUM=0 ,
ARR= 0.00000000 , 2.20000005 , 3.29999995 ,
/
>> xyz = 777 !! invalid line
no match: ignoring input
>> quit
I was considering that concept too (except using a command line input in addition to input files): a separate namelist group defined in each topical module, which in turn are used in the command line module. The command line module would perform a kind of “try/catch” namelist read by checking each namelist group one by one and using iostat to determine if the variable was found and it should stop, or if it should attempt to read the next namelist group. If that works, this seems tighter than the derived-type+pointer approach and probably less work.
I am trying to modernize ODRPACK95, and I would appreciate your recommendation on some matters.
For starters, I will just focus on one question. The code has quite a number of equality comparisons of reals to 0. What is the best to correct these? Should I replace them with comparisons to tiny, or is there a better way? Thanks.
real(kind=wp) :: x
! if (x .eq. 0.0_wp) then
if (abs(x) .lt. tiny(x)) then
! do something
end if
Does x represent the difference of some quantity? Or is it a scalar input value, say set by the user?
It’s mostly some entry in a matrix, for example:
C T: The upper or lower tridiagonal system.
C ZERO: The value 0.0E0_R8.
C***First executable statement DSOLVE
C Find first nonzero diagonal entry in T
J1 = 0
DO 10 J=1,N
IF (J1.EQ.0 .AND. T(J,J).NE.ZERO) THEN ! <=
J1 = J
ELSE IF (T(J,J).EQ.ZERO) THEN ! <=
B(J) = ZERO
END IF
10 CONTINUE
IF (J1.EQ.0) RETURN
C Find last nonzero diagonal entry in T
JN = 0
DO 20 J=N,J1,-1
IF (JN.EQ.0 .AND. T(J,J).NE.ZERO) THEN ! <=
JN = J
ELSE IF (T(J,J).EQ.ZERO) THEN ! <=
B(J) = ZERO
END IF
20 CONTINUE
Context really matters here. Is it a “sentinel” value (i.e. if it’s zero I know it was set to exactly zero), or is it the result of some calculation? If it’s the former I usually do if (abs(x) <= 0), as that really is checking equality to zero, but doesn’t get yelled about when warnings are turned on (and I usually just have a helper function is_zero(x)). If it’s the latter, you need to think about what your tolerance actually is. Are you doing distance checking? At what scale (i.e. cosmological, global, mechanical, cellular)? Is it an increment of some state variable? Again, what’s the typical scale of that property (i.e. MPa, K, neutrons/cm^2, etc.)?
I use the following functions that were inspired by a pseudo code function called almostEqual given in David Kopriva’s, “Implementing Spectral Methods for Partial Differential Equations” book. His version used machine epsilon instead of spacing. I switched to SPACING based on a post on the Intel Fortran Forum back in 2016. Feel free to use and abuse as fits your code. The almostEqual functions allow for an optional user defined tolerance to be specified instead of using the SPACING tests. Both can be used for tests against a literal 0.0
Module equalTests
! Tests for equality to zero or between two real values
USE ISO_FORTRAN_ENV, ONLY: SP=>REAL32, DP=>REAL64
Implicit NONE
PRIVATE
Real(SP), Parameter :: SPACING_FACTOR32 = 2.0_SP
Real(DP), Parameter :: SPACING_FACTOR64 = 2.0_DP
Interface isZero
Module Procedure isZeroR64
Module Procedure isZeroR32
End Interface
Interface almostEqual
Module Procedure almostEqualR64tol
Module Procedure almostEqualR32tol
End Interface
PUBLIC :: isZero, almostEqual
Contains
Elemental Function isZeroR64(anum) Result(a0)
! Function to test if a double precision floating point number is
! almost zero to close to machine precision
! Code modified from post on Intel Fortran Forum 11/15/2016
! Argument variables
Real(DP), Intent(IN) :: anum
Logical :: a0
! Local variables
a0 = (ABS(anum) < SPACING(0.0_DP))
End Function isZeroR64
Elemental Function isZeroR32(anum) Result(a0)
! Function to test if a single precision floating point number is
! almost zero to close to machine precision
! Argument variables
Real(SP), Intent(IN) :: anum
Logical :: a0
! Local variables
a0 = (ABS(anum) < SPACING(0.0_SP))
End Function isZeroR32
Elemental Function almostEqualR32tol(a, b, tol) Result(aeqb)
! Function to test if two REAL32 floating point numbers
! almost equal to close to machine precision
! Code taken from post on Intel Fortran Forum 11/15/2016
! Argument variables
Real(SP), Intent(IN) :: a, b
Real(SP), Optional, Intent(IN) :: tol
Logical :: aeqb
If (PRESENT(tol)) Then
aeqb = ABS(a-b) < tol
Else
aeqb = ABS(a-b) < SPACING_FACTOR32*SPACING(MAX(ABS(a), ABS(b)))
End If
End Function almostEqualR32tol
Elemental Function almostEqualR64tol(a, b, tol) Result(aeqb)
! Function to test if two REAL64 floating point numbers
! almost equal to close to machine precision
! Code taken from post on Intel Fortran Forum 11/15/2016
! Argument variables
Real(DP), Intent(IN) :: a, b
Real(DP), Optional, Intent(IN) :: tol
Logical :: aeqb
If (PRESENT(tol)) Then
aeqb = ABS(a-b) < tol
Else
aeqb = ABS(a-b) < SPACING_FACTOR64*SPACING(MAX(ABS(a), ABS(b)))
End If
End Function almostEqualR64tol
End Module equalTests
As Brad mentioned, context matters.
The snippet shown comes from the LINPACK procedure DSOLVE. Personally I would replace it with DTRSV from BLAS Level 2, and kick it out of the ODRPACK source code entirely, and simply link against an existing version on your system.
Here are the instances you need to replace:
$ grep -n "CALL DSOLVE\(.*\)" ~/fortran/odr.f
9932: CALL DSOLVE(NQ,OMEGA,NQ,WRK1(I,1:NQ,J),4)
9933: CALL DSOLVE(NQ,OMEGA,NQ,WRK1(I,1:NQ,J),2)
9940: CALL DSOLVE(M,WRK4,M,WRK5,4)
9941: CALL DSOLVE(M,WRK4,M,WRK5,2)
9948: CALL DSOLVE(NQ,OMEGA,NQ,TFJACB(I,1:NQ,K),4)
9968: CALL DSOLVE(NQ,OMEGA,NQ,WRK2(I,1:NQ),4)
10109: CALL DSOLVE(M,WRK4,M,WRK5,4)
10110: CALL DSOLVE(M,WRK4,M,WRK5,2)
10129: CALL DSOLVE(M,WRK4,M,T(I,1:M),4)
10130: CALL DSOLVE(M,WRK4,M,T(I,1:M),2)
11621: CALL DSOLVE(M,E,LDE,WRK5,4)
The job argument (argument 5 of dsolve) is either 2 or 4, corresponding to the operations:
C 2 Solve T*X=B, T upper triangular,
C 4 Solve trans(T)*X=B, T upper triangular.
The corresponding BLAS calls would be
! Old
call dsolve(n,t,ldt,b,job)
! New
call dtrsv('U','N','N',n,t,ldt,b,1) ! JOB = 2 -> Upper-triangular, Notranspose, Nonunit-diagonal
call dtrsv('U','T','N',n,t,ldt,b,1) ! JOB = 4 -> Upper-triangular, Transpose, Nonunit-diagonal
This test looks like it is to detect denormal numbers. This would be very different than a test against a given tol value, or even a test using spacing() or x*epsilon(x) for the relative difference of a value.
Thanks, that seems a suitable solution for those particular lines of code.
However, the use of exact equality/inequality comparisons to reals is far more prevalent. There are 53 occurrences, as one may conclude by searching for “FPT - 3087” in odr.f90.
So, I suppose I will have to address each case individually, most likely also using some of the solutions proposed by @rwmsu and @everythingfunctional.