Numerical experiments, Tips, Tricks and Gotchas
All classical numerical algorithms were first implemented in FORTRAN and ALGOL [1-6] and later translated to many other languages. Numerous robust, efficient and reliable algorithms in FORTRAN still are freely available from various repositories [7-13].
The algorithms usually are well documented and include testing programs and data. Sometimes the old FORTRAN libraries are the only source of certain algorithms. If necessary, I borrowed and translated some algorithms to Pascal, C++, and C#.
There are many numerical packages in Python. However, the ability to convert the old good FORTRAN algorithm is still of interest (at least for me).
The developed conversion algorithm is focused on Fortran 90. This version is closer to modern languages, which makes it easier to convert using simple string manipulation. On the other hand, a lot of code was already converted to this (and later) versions. The earlier versions will be addressed in a separate project.
The algorithm includes the folliwing steps.
Fortran | Python | Remark |
---|---|---|
'!' | '#' | Comment |
'\t' | ' ' | Tab → 2 spaces |
'.d0' | '.0' | |
'd0' | '0' | |
'.eq.' | ' == ' | |
'.ne.' | ' != ' | |
'/=' | ' != ' | |
'.lt.' | ' < ' | |
'.gt.' | ' > ' | |
'.le.' | ' <= ' | |
'.ge.' | ' >= ' | |
'then' | ':' | |
'else if' | 'elif' | |
'else' | 'else:' | |
'&' | '\\' | Continuation |
'.and.' | ' and ' | |
'.or.' | ' or ' | |
'.true.' | ' True ' | |
'.false.' | ' False ' | |
'end if' | '' | |
'do while' | 'while' | |
'program' | 'def' | |
'function' | 'def' | |
'subroutine' | 'def' |
The replacements are stored as a list of pairs.
map_list=[('!','#'), ('\t',' '), ('.d0','.0'), ('d0','0'), ('.eq.',' == '), . . . . . . . . .
The below case-insensitive replace function is used. It skips comments (from '#' to the rest of a line).
def ireplace(old, new, text): ''' Case Insensitive Replace excluding comments based on http://stackoverflow.com/questions/919056/python-case-insensitive-replace ''' idx = 0 lim = text.find('#') if lim < 0: lim = len(text) while idx < lim: index_l = text.lower().find(old.lower(), idx) if index_l == -1: return text text = text[:index_l] + new + text[index_l + len(old):] idx = index_l + len(old) return text
The list of double precision functions (incomplete).
Fortran | Python | Remark |
---|---|---|
'dabs' | 'abs' | |
'dsqrt' | 'math.sqrt' | |
'dlog' | 'math.log' | |
'dexp' | 'math.exp' |
math_list=[ ('dabs','abs'), ('dsqrt','math.sqrt'), ('dlog','math.log'), ('dexp','math.exp'), ]
Fortran | Python | Remark |
---|---|---|
'end do' | '' | Can be '# end do' |
'end program' | '# end program' | |
'end function' | '# end function' | |
'end subroutine' | '# end subroutine' | |
'end' | '# end' | |
'integer' | '#integer' | |
'real' | '#real' | |
'allocate' | '#allocate' |
There are two options:
Te latter can be helpful for manual formatting and adjustment.
This function is called after statement and math function substitutions:
def process_do(line): ''' replaice Fortran do to Python for: do i=1, N => for i in range(0, N): do j=1, M, 2 => for j in range(0, M, 2): do iq=ip+1, N => for iq in range(ip+1, N): ''' if line.startswith('#'): return line if line.strip().startswith( 'do ' ): if '#' in line: do_var_lims,comment = line.split('#') comment = ' #'+''.join(comment) else: do_var_lims,comment = line, '' do_var,lims = do_var_lims.split('=') for_var = do_var.replace('do','for') + ' in range(' lims=lims.split(',') lims=[i.strip() for i in lims] if str(lims[0])=='1': # heuristic correction for 0-based arrays lims[0]='0' # warning: this can be unnecessery sometimes. result = for_var + ', '.join(lims) + '):' if comment: result += comment return result else: return line
The loops in Fortran are usually used for processing arrays.
The arrays in Fortran are 1-bsed and in Python are 0-based.
Therefore, in addition to conversion from Fortran 'do' loops to Python 'for' loops, the indices are also shifted:
do i=1, N → for i in range(0, N)
.
The indices are not shifted if a lower limit is a variable: do i=j, N → for i in range(j, N)
.
Note that this heurictic correction can be unnecessary sometimes.
do i=1, 50 sm=0.d0 do ip=1, N-1 !sum off-diagonal elements do iq=ip+1, N sm=sm+DABS(A(ip,iq)) end do end do . . . . .
The above Fortran code is converted to
for i in range(0, 50): sm=0.0 for ip in range(0, N-1): #sum off-diagonal elements for iq in range(ip+1, N): sm=sm+abs(A[ip,iq]) . . . . . . .
This function replaces parentheses with square brackets.
def adjust_array(line, arr): ''' Conversion of Fortran array access to Python: line = 'A(i,j) = A(m,A(k,l))' arr = 'A' result = 'A[i,j] = A[m,A[k,l]]' ''' if line.startswith('#'): return line i = line.find(arr+'(') while i >= 0: j = line.find('(',i) line = line[:j] + '[' + line[j+1:] c = 1 for k in range(j+1,len(line)): if line[k] == '(': c += 1 if line[k] == ')': c -= 1 if c == 0: line = line[:k] + ']' + line[k+1:] # i = line.find(arr+'(', k+1) # does not process nested arrays i = line.find(arr+'(') break return line
This function transforms Fortran function/subroutine declarations to Python ones.
def adjust_functions(content): """ Adds ':' after ')' """ for n,line in enumerate(content): count = 0 if line.strip().startswith('def'): i = line.find('(') if i >= 0: count = 1 for k in range(i,len(line)): #print(k, line[k]) if line[k] == '#': break if line[k] == ')': count = 0 content[n] = line[:k+1] + ':' + line[k+1:] break else: count = 0 i = line.find('#') if i >= 0: content[n] = line[:i] + ':' + line[i:] else: content[n] = line + ':' else: if count > 0: i = 0 for k in range(i,len(line)): if line[k] == '#': break if line[k] == ')': count = 0 content[n] = line[:k+1] + ':' + line[k+1:] break return content
Subroutine Jacobi(A,N,D,V,NROT) . . . . .
def Jacobi(A,N,D,V,NROT): . . . . .
The algorithm was implemented as IPython HTML notebook. The notebook can be viewed or downloaded.
The project started when I needed the Jacobi diagonalization algorithm. I found this Fortran 90 implementation ujacobi.f90.
The initial file with syntax highlighting: ujacobi.f90.
The resulting Python file: ujacobi.py.
Obviously, the generated Python file requires some additional manual adjustment and formatting. It also requires some refactoring.
In particular, the loop limits at lines 35, 73, 79 must be adjusted. The line 107 must be commented out. The allocations should be replaced with array creation.The working version and the test (based on [15]) are available at Downloads below.
Download project files (f90_to_py.ipynb, ujacobi.f90, ujacobi.py, jacobi.py, test_jacobi.py, f90_to_py.py).
I would like to thank Prof. Jianjun Hu, University of South Carolina, for improving the code in section 2.4.
© Nikolai Shokhirev, 2012-2024
email: nikolai(dot)shokhirev(at)gmail(dot)com