A recent comment from a blog reader made me relook at the patran neutral mesh file writer routine that I had on this blog.
Since the code was pasted on this blog, it lacked proper update and ability to have a test case or example file.
So spent this weekend, creating an example file and have uploaded the code and the example on github.
The entire writer is now availabe at https://github.com/sukhbinder/PatranMeshWriter
Hope to add and update the code more often. 🙂
Â
Most of my day job involves Fortran code and finding quick details is a regular need. Here are four regex in python that I have been using for some time now. Sharing it here.
To Find Subroutines and Functions
import re procedule_name=re.compile(r'\s*(RECURSIVE)?\s+(SUBROUTINE|FUNCTION)\s+\S+\(.*',re.IGNORECASE)
Â
To Find All Call Statements
callname=re.compile(r'^\s*call\s*\S*[(]?\S*[)]?',re.IGNORECASE)
   Â
Â
To Find Variables (Integer, Double Precision, Logical, Real)
variregex=re.compile(r'^\s+(logical|integer|double.+precision|real(\*\s*\d)?)\s+\S*',re.IGNORECASE)
Â
To Find Character Variables
characteregex=re.compile(r'^\s+(character+([\*\s*\d+]|\(\S*)*)\s+\S*',re.IGNORECASE)
Â
Suppose you got a long fortran file full of code, and want to parse it and get the subroutine/functions defined in the file along with its argument.
Here’s a simple python code to extract all this info in one go
import re
from collections import defaultdict
def processsubroutine(line):
if ')' in line:
return line.split("(")[1].split(')')[0].split(',')
def collectsubroutines(fname):
procnam=defaultdict()
fort_proc_def = re.compile(r'\s*(RECURSIVE)?\s*(SUBROUTINE|FUNCTION)\s+\S+\(*', re.IGNORECASE)
fp = open(fname,"r")
datalines=fp.readlines()
fp.close()
for i,line in enumerate(datalines):
if fort_proc_def.match(line):
procnam[(line.split()[1].split('(')[0])]=processsubroutine(line)
return procnam
def main():
for key,item in procnam.items():
print key,item
Here’s an implementation of the Lights out game in fortran. Wrote this in 2011 and found it yesterday. Works as it is and can be improved. Checkout the screenshot.
For those who don’t know, here’s a brief description from wikipedia
The game consists of a 5 by 5 grid of lights. When the game starts, a random number or a stored pattern of these lights is switched on. Pressing any of the lights will toggle it and the four adjacent lights. The goal of the puzzle is to switch all the lights off, preferably in as few button presses as possible.
module lightsout
implicit none
integer :: elements(0:25)
integer :: iswon,moves
logical,parameter :: debug =.false.
character*30 :: sdata(4)
contains
subroutine printtable1()
integer :: i,j,ii,jj
sdata(1)=' _\|/_ '
sdata(2)=' (o o) '
sdata(3)='+----oOO-{_}-OOo----+'
do i=1,3
write(*,'(a30)') sdata(i)
end do
write(*,'(6(2x,a))') ' ','a','b','c','d','e'
ii=0
jj=4
do j=1,5
! write(*,'(6i3)') j,(elements(i),i=ii,jj)
write(*,'(i3,5a3)') j,(achar(elements(i)+46),i=ii,jj)
ii=jj+1
jj=ii+4
end do
! if(debug) then
! print *,"Degub mode:"
! write(*,'(5i8)')(elements(i),i=0,24)
! end if
end subroutine
subroutine printtable()
integer :: i,j,ii,jj
sdata(1)=' _\|/_ '
sdata(2)=' (o o) '
sdata(3)='+----oOO-{_}-OOo----+'
do i=1,3
write(*,'(a30)') sdata(i)
end do
write(*,'(6(1x,a,2x))') ' ','a','b','c','d','e'
ii=1
write(*,'(7(5a))')' ','___','.','___','.','___','.','___','.','___','.'
do i=1,5
do j=1,5
if(elements(ii) .eq. 0) then
write(*,'(4a)',advance='no') achar(32),achar(32),achar(32),'|'
else
write(*,'(4a)',advance='no')achar(219),achar(219),achar(219),'|'
end if
ii=ii+1
end do
write(*,*)
write(*,'(i3,5(5a))')i,'.','___','|','___','|','___','|','___','|','___','|'
end do
! if(debug) then
! print *,"Degub mode:"
! write(*,'(5i8)')(elements(i),i=0,24)
! end if
end subroutine
subroutine fillall()
integer :: i
do i =0 , 24
elements(i)=1
end do
end subroutine
subroutine clearall()
integer :: i
do i = 0, 24
elements(i)=0
end do
end subroutine
subroutine level1()
iswon=0
moves=0
call clearall()
elements(10)=1
elements(12)=1
elements(14)=1
end subroutine
subroutine level2()
iswon=0
moves=0
call clearall();
elements (12)=1;
elements (16)=1;
elements (17)=1;
elements (18)=1;
elements (20)=1;
elements (21)=1;
elements (22)=1;
elements (23)=1;
elements (24)=1;
end subroutine
subroutine level3()
iswon=0
moves=0
call fillall();
elements (4)=0;
elements (6)=0;
elements (7)=0;
elements (8)=0;
elements (11)=0;
elements (12)=0;
elements (13)=0;
elements (16)=0;
elements (17)=0;
elements (18)=0;
elements (24)=0;
end subroutine
subroutine level4()
iswon=0
moves=0
call clearall();
elements (2)=1;
elements (6)=1;
elements (8)=1;
elements (10)=1;
elements (12)=1;
elements (14)=1;
elements (16)=1;
elements (18)=1;
elements (22)=1;
end subroutine
subroutine level5()
iswon=0
moves=0
call clearall();
elements (11)=1;
elements (16)=1;
elements (21)=1;
end subroutine
subroutine checkall()
integer :: win,i
win=1;
do i = 0 ,24
if (elements (i)==0) then
win=0;
end if
end do
if (win==1) then
call printtable()
write(*,*) "You Won!!"
write(*,*) "You did it in ",moves, "moves"
iswon=1
end if
if(debug) print *,'7',iswon
end subroutine
subroutine check(v)
integer :: q,w,row,intv,a,b,c,d,kc,kd,kv,is,v
q=5;
w=1;
moves=moves+1
row=Int(v/q)+w;
intv =Int(v);
a=intv+q;
b=intv-q;
c=intv+w;
d=intv-w;
if(debug) print *, 'check',iswon
if(debug) print *, 'row a b c d', row,a,b,c,d
if (a<0 .or. a>24) then
a=25;
end if
if (b<0 .or. b>24) then
b=25;
end if
if (c<0 .or. c>24) then
c=25;
end if
if (d<0 .or. d>24) then
d=25;
end if
kc = (Int(c/q)+w)
kd = (Int(d/q)+w)
kv = row;
if (kc /= kv) then
c=25;
end if
if (kd /= kv) then
d=25;
end if
if (v==5) d=25;
if(debug) print *,'a b c d second', a,b,c,d
if ( elements (a) == 1) then
elements (a)=0;
else
elements (a)=1;
end if
if ( elements (b) == 1) then
elements (b)=0;
else
elements (b)=1;
end if
if ( elements (c) == 1) then
elements (c)=0;
else
elements (c)=1;
end if
if ( elements (d) == 1) then
elements (d)=0;
else
elements (d)=1;
end if
if ( elements (v) == 1) then
elements (v)=0;
else
elements (v)=1;
end if
if(debug) print *,'6',iswon
call checkall()
if(debug) print *, 'check',iswon
end subroutine
end module
program test
use lightsout
character*2 ::ss
integer :: ival,ilev
if(debug) call fillall()
if(debug) call printtable()
if(debug) print *, "*******"
if(debug) call clearall()
if(debug) call printtable()
if(debug) print *, "*******"
if(debug) call level1()
if(debug) call printtable()
if(debug) print *, "*******"
if(debug) call level2()
if(debug) call printtable()
if(debug) print *, "*******"
if(debug) call level3()
if(debug) call printtable()
if(debug) print *, "*******"
if(debug) call level4()
if(debug) call printtable()
if(debug) print *, "*******"
if(debug) call level5()
if(debug) call printtable()
ival=0
ilev=1
call level1()
do while(.true.)
if(debug) print *,iswon,ilev
if(iswon .ne. 0) then
ilev=ilev+1
print *, "Entering level :" ,ilev
if (ilev == 2) call level2
if (ilev == 3) call level3
if (ilev == 4) call level4
if (ilev == 5) call level5
end if
call printtable()
if(debug) print *,'2',iswon,ilev
write(*,'(a1)',advance="no") ":"
read(*,'(2a)') ss
call givevalue(ss,ival)
if(debug) print *, ival
if (ival .eq. 999) exit
if(debug) print *,'3',iswon,ilev
if(ival .ne. -1) call check(ival)
if(debug) print *,'4',iswon,ilev
end do
contains
subroutine givevalue(sstext,value)
character*2 :: sstext
integer :: value,i,k
if (sstext(1:1) .eq. 'q') then
value=999
end if
select case(sstext(1:1))
case('a')
i=1
case ('b')
i=2
case ('c')
i=3
case ('d')
i=4
case ('e')
i=5
case default
i=-1
value=-1
end select
read(sstext(2:2),'(i1)',err=100) k
if (k <7 .and. k .ne. 0) then
value=(i-1)+(k-1)*5
else
value=-1
k=-1
end if
if(debug) print *, k,i,value
return
100 value =999
end subroutine
end program
As an aeronautical engineer, I received fortran as a legacy. From NAL to working for Rolls-Royce, work life revolves around fortran. This  is the language that I have been programming for 12+ years now. That explains this many fortran related posts here and here.
But from couple of years, I am falling in love with Python. Here in this post, I will walk through the steps to use fortran and python together.
There are many times, when we don’t want to use the methods included in Python and need to use the trusted, optimised, time tested fortran library. This post will show how to create a python binding of the fortran code and use it in python.
I will present the complete code first and then explain all the parts of it.
I will use this simple NACA code library.
naca4_sym(c,t,n) evaluates y(x) for a NACA symmetric 4-digit airfoil. C the chord length. T the maximum relative thickness. N the number of sample points.
file name: naca.f
This is the actual fortran code to be used from python.
subroutine naca4_symmetric ( t, c, n, x, y )
!*****************************************************************************80
!
!! NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 21 May 2014
!
! Author:
!
! John Burkardt
!
! Reference:
!
! Eastman Jacobs, Kenneth Ward, Robert Pinkerton,
! "The characteristics of 78 related airfoil sections from tests in
! the variable-density wind tunnel",
! NACA Report 460, 1933.
!
! Parameters:
!
! Input, real ( kind = 8 ) T, the maximum relative thickness.
!
! Input, real ( kind = 8 ) C, the chord length.
!
! Input, integer ( kind = 4 ) N, the number of sample points.
!
! Input, real ( kind = 8 ) X(N), points along the chord length.
! 0.0 <= X(*) <= C.
!
! Output, real ( kind = 8 ) Y(N), for each value of X, the corresponding
! value of Y so that (X,Y) is on the upper wing surface, and (X,-Y) is on the
! lower wing surface.
!
implicit none
integer ( kind = 4 ) n
real ( kind = 8 ) c
real ( kind = 8 ) t
real ( kind = 8 ) x(n)
real ( kind = 8 ) y(n)
y(1:n) = 5.0D+00 * t * c * ( &
0.2969D+00 * sqrt ( x(1:n) / c ) &
+ (((( &
– 0.1015D+00 ) * ( x(1:n) / c ) &
+ 0.2843D+00 ) * ( x(1:n) / c ) &
– 0.3516D+00 ) * ( x(1:n) / c ) &
– 0.1260D+00 ) * ( x(1:n) / c ) )
return
end
file name: aefoilnaca.py
This is the python wrapper using the nacalib created from fortran code
def naca4_sym(c,t,n): """ NACA4_SYMMETRIC evaluates y(x) for a NACA symmetric 4-digit airfoil. C the chord length. T the maximum relative thickness. N the number of sample points. python wrapper for the code by John Burkardt by Sukhbinder Singh """ import numpy as np import nacalib x = np.linspace(0.0,c,n) y=np.zeros(n) y=nacalib.naca4_symmetric(t,c,x,n) xy=np.zeros((2,2*n)) xy[0,0:n]=x xy[0,n:]=x[n::-1] xy[1,0:n]=-y xy[1,n:]=y[n::-1] return xy
file name: test.py
This file tests the python wrapper module.
""" Testing aefoilnaca """ import aefoilnaca import matplotlib.pyplot as plt c = 10.0 # Chord length t=0.15 # the maximum relative thickness n=80 # number of point m=0.02 # the maximum camber p=0.4 # the location of maximum camber y = aefoilnaca.naca4_sym(c,t,n) plt.plot(y[0,:],y[1,:]) plt.show()
The steps:
First, we create a signature file from naca.f by running
f2py naca.f -m nacalib -h naca.pyf
The signature file is saved to naca.pyf
Second step edit naca.pyf file.
Next we’ll open the naca.pyf and tell F2PY that the argument t,c,x,n are input arguments (use intent(in) attribute) and that the result, i.e. the contents of y after calling Fortran function naca4_symmetric should be returned to Python (use intent(out) attribute). In addition, an array a should be created dynamically using the size given by the input argument n (use depend(n) attribute to indicate dependence relation).
This is the content of the modified naca.pyf
file name: naca.pyf
! -*- f90 -*-
! Note: the context of this file is case sensitive.
python module nacalib
interface
subroutine naca4_symmetric(t,c,n,x,y)
real(kind=8),intent(in) :: t
real(kind=8),intent(in) :: c
integer(kind=4), optional,check(len(x)>=n),depend(x) :: n=len(x)
real(kind=8) dimension(n),intent(in) :: x
real(kind=8) dimension(n),depend(n),intent(out) :: y
end subroutine naca4_symmetric
end interface
end python module nacalib
! This file was auto-generated with f2py (version:2).
! See http://cens.ioc.ee/projects/f2py2e/
And finally, we build the extension module by running
f2py -c naca.pyf naca.f
This creates the naca.so file and it can be used as shown in the aefoilnaca.py
import nacalib x = np.linspace(0.0,c,n) y=np.zeros(n) y=nacalib.naca4_symmetric(t,c,x,n)
Hope this is clear and simple, if you have any questions, feel free to comment.
All this code is available at my github account
For further reading refer to this two excellent reference
1. Python as a glue
2. Using f2py – The smart way
A few months ago I posted this pgplot powered fortran animation program on calculating the approximate value of PI by Monte Carlo method using random numbers.
And the other found myself reading about random or stochastic phenomenons. After reading about stochastic processes, I was invariably reminded of the stock market. The most stochastic place one can imagine.
After a bit of googling and a bit of tinkering, created the following program.
Given the annual drift, volatility, initial price and number of trials, the program simulates the movement of a stock price.
program MonteCarloStocks
double precision a,b,c
character(len=100) buffer
integer i,n,noarg,IARGC
noarg=IARGC()
if(noarg == 4) then
CALL GETARG (1, buffer)
read(buffer,'(f6.2)') a
CALL GETARG (2, buffer)
read(buffer,'(f6.2)') b
CALL GETARG (3, buffer)
read(buffer,'(f6.2)') c
CALL GETARG (4, buffer)
read(buffer,'(i7)') i
print *, a,b,c,i
else
write(*,'(a)',advance='no') "Please enter annual drift :"
read(*,'(f8.2)')a
write(*,'(a)',advance='no') "Please enter annual Volatility :"
read(*,'(f8.2)')b
write(*,'(a)',advance='no') "Please enter Initial Price :"
read(*,'(f8.2)')c
write(*,'(a)',advance='no') "Please enter no of iteration :"
read(*,'(i7)')i
end if
call monte(a,b,c,i)
call system("stock.xls")
end program
subroutine monte(driftyearly,volatility,initialprice,not)
double precision driftyearly,volatility,initialprice
double precision driftdaily,voldaily,p,pri(252),PPND16
integer,parameter :: nod=252
integer not
driftyearly = driftyearly/100
volatility = volatility/100
driftdaily=driftyearly/nod
voldaily=volatility/(nod**0.5)
open(200,file="fort.200")
call random_seed
if(not <= 0) not=1
if(not >=25000) not =25000
do j=1,not
pri(1)=initialprice
do i=2,nod
call random_number(p)
p=PPND16 (p, IFAULT)
pri(i)=pri(i-1)*EXP(driftdaily+voldaily*p)
end do
write(200,'(252f8.2)')pri
end do
close(200)
end subroutine
DOUBLE PRECISION FUNCTION PPND16 (P, IFAULT)
!
! ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
!
! Produces the normal deviate Z corresponding to a given lower
! tail area of P; Z is accurate to about 1 part in 10**16.
!
! The hash sums below are the sums of the mantissas of the
! coefficients. They are included for use in checking
! transcription.
!
DOUBLE PRECISION ZERO, ONE, HALF, SPLIT1, SPLIT2, CONST1, &
CONST2, A0, A1, A2, A3, A4, A5, A6, A7, B1, B2, B3, &
B4, B5, B6, B7, &
C0, C1, C2, C3, C4, C5, C6, C7, D1, D2, D3, D4, D5, &
D6, D7, E0, E1, E2, E3, E4, E5, E6, E7, F1, F2, F3, &
F4, F5, F6, F7, P, Q, R
PARAMETER (ZERO = 0.D0, ONE = 1.D0, HALF = 0.5D0, &
SPLIT1 = 0.425D0, SPLIT2 = 5.D0, &
CONST1 = 0.180625D0, CONST2 = 1.6D0)
!
! Coefficients for P close to 0.5
!
PARAMETER (A0 = 3.3871328727963666080D0, &
A1 = 1.3314166789178437745D+2,&
A2=1.9715909503065514427D+3,&
A3=1.3731693765509461125D+4,&
A4=4.5921953931549871457D+4,&
A5=6.7265770927008700853D+4,&
A6=3.3430575583588128105D+4,&
A7=2.5090809287301226727D+3,&
B1=4.2313330701600911252D+1,&
B2=6.8718700749205790830D+2,&
B3=5.3941960214247511077D+3,&
B4=2.1213794301586595867D+4,&
B5=3.9307895800092710610D+4,&
B6=2.8729085735721942674D+4,&
B7=5.2264952788528545610D+3)
! HASH SUM AB 55.88319 28806 14901 4439
!
! Coefficients for P not close to 0, 0.5 or 1.
!
PARAMETER (C0 = 1.42343711074968357734D0,&
C1=4.63033784615654529590D0,&
C2=5.76949722146069140550D0,&
C3=3.64784832476320460504D0,&
C4=1.27045825245236838258D0,&
C5=2.41780725177450611770D-1,&
C6=2.27238449892691845833D-2,&
C7=7.74545014278341407640D-4,&
D1=2.05319162663775882187D0,&
D2=1.67638483018380384940D0,&
D3=6.89767334985100004550D-1,&
D4=1.48103976427480074590D-1,&
D5=1.51986665636164571966D-2,&
D6=5.47593808499534494600D-4,&
D7=1.05075007164441684324D-9)
! HASHSUMCD49.33206503301610289036
!
! CoefficientsforPnear0or1.
!
PARAMETER(E0=6.65790464350110377720D0,&
E1=5.46378491116411436990D0,&
E2=1.78482653991729133580D0,&
E3=2.96560571828504891230D-1,&
E4=2.65321895265761230930D-2,&
E5=1.24266094738807843860D-3,&
E6=2.71155556874348757815D-5,&
E7=2.01033439929228813265D-7,&
F1=5.99832206555887937690D-1,&
F2=1.36929880922735805310D-1,&
F3=1.48753612908506148525D-2,&
F4=7.86869131145613259100D-4,&
F5=1.84631831751005468180D-5,&
F6=1.42151175831644588870D-7,&
F7=2.04426310338993978564D-15)
! HASH SUM EF 47.52583 31754 92896 71629
!
IFAULT = 0
Q = P - HALF
IF (ABS(Q) .LE. SPLIT1) THEN
R = CONST1 - Q * Q
PPND16 = Q * (((((((A7 * R + A6) * R + A5) * R + A4) * R + A3) &
* R + A2) * R + A1) * R + A0) / &
(((((((B7 * R + B6) * R + B5) * R + B4) * R + B3) &
* R + B2) * R + B1) * R + ONE)
RETURN
ELSE
IF (Q .LT. ZERO) THEN
R = P
ELSE
R = ONE - P
END IF
IF (R .LE. ZERO) THEN
IFAULT = 1
PPND16 = ZERO
RETURN
END IF
R = SQRT(-LOG(R))
IF (R .LE. SPLIT2) THEN
R = R - CONST2
PPND16 = (((((((C7 * R + C6) * R + C5) * R + C4) * R + C3) &
* R + C2) * R + C1) * R + C0) / &
(((((((D7 * R + D6) * R + D5) * R + D4) * R + D3) &
* R + D2) * R + D1) * R + ONE)
ELSE
R = R - SPLIT2
PPND16 = (((((((E7 * R + E6) * R + E5) * R + E4) * R + E3) &
* R + E2) * R + E1) * R + E0) / &
(((((((F7 * R + F6) * R + F5) * R + F4) * R + F3) &
* R + F2) * R + F1) * R + ONE)
END IF
IF (Q .LT. ZERO) PPND16 = - PPND16
RETURN
END IF
END
Cool I would say.
An example of Monte Carlo method in action.
This is a brain game that you might enjoy. In the chess board you can move the knight as if you’re playing a regular game of chess, but you need to cover the whole board with the least possible moves, a good score would be 65 moves or less. This is a very good training for strategy games.
Written in fortran. Just for fun!! As always work in progress..
Stumbled upon this link from an old bookmark.
The first FORTRAN compiler was a milestone in the history of computing,
at that time computers had very small memories (on the order of 15KB,
it was common then to count memory capacities in bits), they were slow
and had very primitive operating systems (if they had them at all).
At those days it seemed that the only practical way is to program in
assembly language.The pioneers of FORTRAN didn’t invent the idea of writing programs in a
High Level Language (HLL) and compiling the source code to object code
with an optimizing compiler, but they produced the first successful HLL.
They designed an HLL that is still widely used, and an optimizing compiler
that produced very efficient code, in fact the FORTRAN I compiler held
the record for optimizing code for 20 years!This wonderful first FORTRAN compiler was designed and written from
scratch in 1954-57 by an IBM team lead by John W. Backus and staffed with
super-programmers like Sheldon F. Best, Harlan Herrick, Peter Sheridan,
Roy Nutt, Robert Nelson, Irving Ziller, Richard Goldberg, Lois Haibt
and David Sayre. By the way, Backus was also system co-designer of the
computer that run the first compiler, the IBM 704.The new invention caught quickly, no wonder, programs computing nuclear
power reactor parameters took now hours instead of weeks to write, and
required much less programming skill. Another great advantage of the new
invention was that programs now became portable. Fortran won the battle
against Assembly language, the first in a series of battles to come,
and was adopted by the scientific and military communities and used
extensively in the Space Program and military projects.The phenomenal success of the FORTRAN I team, can be attributed in part
to the friendly non-authoritative group climate. Another factor may be
that IBM management had the sense to shelter and protect the group,
even though the project took much more time than was first anticipated.
Interesting, read the complete article here
And while you are at it, don’t miss reading the FORTRAN SAGA
If you are involved in a project thats calls legacy Fortran code from c sharp, then this is one post (Accessing Fortran Legacy DLL in C#) by Ranjeet Sharma that you shouldn’t miss.
These are few important considerations from the post
Important Considerations:
The FORTRAN DLL cannot be accessed from C# using Add Reference as it is not a managed DLL. It can be accessed only through P/Invoke syntax.
When doing mixed-language programming with FORTRAN and C# we must be concerned with data types, because each language handles them differently.
Along with the Data type differences there is a marked difference between C# and FORTRAN in the way the Memory is allocated to the Arrays.C# used a Row Major memory allocation i.e. data of a row is stored in Contiguous memory locations however Fortran uses a Column Major memory allocation i.e. data of a column is stored in contiguous memory locations. So when sending arrays from C# to FORTRAN we need to transpose the array so that FORTRAN understands the data send to it.
In C# we distinguish between parameters as being passed either by reference or by value. In FORTRAN all parameters are passed by reference only. Hence every scalar parameter to a FORTRAN routine has to be qualified by the keyword ref. This does not apply to arrays as arrays are passed by reference in C#.
Complex type as such does not exist in C#. A complex number Z may be represented as an ordered pair of real numbers consisting of the real part of Z and the imaginary part of Z. In FORTRAN this is exactly how a complex type is stored .Specifically we may use a C# double array to hold the real and imaginary parts (in that order) of complex numbers and pass this into the FORTRAN DLL to represent a complex array or complex number.
If FORTRAN DLL subroutine uses a character Array then to call that subroutine in C# we need to pass two parameters for that First parameters will actually hold the Character content and second parameters will hold the length of the string builder.
For the picture on data type information and sample code visit the original post.
Most of considerations listed here are applicable for mixing Fortran with any other language, be it Java, Python, VBA or C++.
Fortran 90 with Excel – gfortran example is one of the post that constantly gets lot of attention as measured by the top posts sidebar widget on this blog.
So thought about revisiting it and this post is to show some toy programs that shows how to link fortran program as a dll with excel.
I am using gfortran and the following simple fortran programs to create a DLL which is called from excel.
A few years back, while learning to code in a new in-house Finite element software, I had to display some results in the GUI. And the GUI didn’t allow any GUI display unless a mesh or model was loaded.
With no FE model available, I hacked it by creating a on the fly one 8 noded, single element brick mesh and tricked the software to get my results displayed.
What I did was write a simple patran neutral mesh file. Over the years I have built on that idea and here’s a fortran module that can be used to generate a patran neural mesh file given a set of elements and node coordinates.
Was cleaning my home laptop and found this fortran transformation module that I wrote long time ago.
Has all the transformation that are generally used. Just posting it in case, I need it again.
The old one still wins
Well Dan won’t be happy with it. I am not very happy with it. But While thinking about optimising the old fortran binary stl reader for speed, I knew the speed can be reduced significantly if I can reduce the number of I/O calls in the routine. For that need a mechanism to get similar functionality like seek and tell in fortran.
After researching it, found this good notes on the new stream I/O feature of modern fortran. And using Dr. Clive Page’s some notes on Stream I/O., built this new version of binary stl file reader. The present code is compact, straight-forward and easy to understand but isn’t speedier.
so the Quest continues.
program stlread
character header*80,filname*80
real*4 n(3),x1(3),x2(3),x3(3)
integer*2 padding
integer*4 ntri,iunit
real*4, allocatable :: normals(:,:),triangles(:,:)
integer iargc
if (iargc() .ge. 1) then
call getarg(1,filname)
else
filname="/Users/Sukhbinder/Desktop/projects/projects/work/stlreader/zDino/Dino.stl"
end if
iunit=13
print *, filname
open(unit=iunit,file=filname,status='old',access='stream', &
& form='unformatted')
read(iunit)header
read(iunit) ntri
allocate(normals(3,ntri))
allocate(triangles(3,ntri*3))
k=1
do i = 1,ntri
read(iunit) normals(1,i),normals(2,i),normals(3,i)
read(iunit) triangles(1,k),triangles(2,k),triangles(3,k)
read(iunit) triangles(1,k+1),triangles(2,k+1),triangles(3,k+1)
read(iunit) triangles(1,k+2),triangles(2,k+2),triangles(3,k+2)
read(iunit) padding
k=k+3
end do
close(iunit)
print *, header
print *, ntri
print *, normals(:,1)
print *, triangles(:,1)
print *, triangles(:,2)
print *, triangles(:,3)
end
And here’s the old one.
program stlread
! Program to Read Binary STL file.
! programmed by Sukhbinder Singh
! stores normals and the triangles in the normals and triangles arrays.
!
! 23rd May 2011
real*4 n(3),x1(3),x2(3),x3(3)
integer*2 padding
integer*4 ntri,iunit,irc
character(len=80) title, filename
character(len=52) sentry
double precision, allocatable :: normals(:,:),triangles(:,:)
integer iargc
if (iargc() .ge. 1) then
call getarg(1,filename)
else
filename="/Users/Sukhbinder/Desktop/projects/projects/work/stlreader/zDino/Dino.stl"
end if
iunit=13
open(unit=iunit,file=filename,status='old',form='unformatted', &
& access='direct',recl=80)
read(iunit,rec=1)title
close(iunit,status='keep')
open(unit=iunit,file=filename,status='old',form='unformatted', &
& access='direct',recl=4)
read(iunit,rec=21)ntri
close(iunit,status='keep')
open(unit=iunit,file=filename,status='old',form='unformatted', &
& access='direct',recl=2)
irc=(80+4)/2+1
!call fseek(iunit,84,0,ierr)
allocate(normals(3,ntri))
allocate(triangles(3,ntri*3))
!read(iunit,rec=1)sentry
!print *, sentry
k=1
do i=1,2 !ntri
call readbin(iunit,n(1),irc)
call readbin(iunit,n(2),irc)
call readbin(iunit,n(3),irc)
call readbin(iunit,x1(1),irc)
call readbin(iunit,x1(2),irc)
call readbin(iunit,x1(2),irc)
call readbin(iunit,x2(1),irc)
call readbin(iunit,x2(2),irc)
call readbin(iunit,x2(3),irc)
call readbin(iunit,x3(1),irc)
call readbin(iunit,x3(2),irc)
call readbin(iunit,x3(3),irc)
normals(1,i)=n(1)
normals(2,i)=n(2)
normals(3,i)=n(3)
triangles(1,k)=x1(1)
triangles(2,k)=x1(2)
triangles(3,k)=x1(3)
triangles(1,k+1)=x2(1)
triangles(2,k+1)=x2(2)
triangles(3,k+1)=x2(3)
triangles(1,k+2)=x3(1)
triangles(2,k+2)=x3(2)
triangles(3,k+2)=x3(3)
if (i .eq. 1) then
write(*,*) "normals: ",n
write(*,*) "v1 : ",x1
write(*,*) "v2 : ",x2
write(*,*) "v3 : ",x3
end if
k=k+3
read(iunit,rec=irc)padding
irc=irc+1
end do
write(*,*) trim(filename),' has this title ',trim(title),' and has ',ntri, ' triangles'
end program
subroutine readbin(iunit,a,irc)
real*4 a,b
integer*2 ib(2)
equivalence(b,ib)
read(iunit,rec=irc)ib(1)
irc=irc+1
read(iunit,rec=irc)ib(2)
irc=irc+1
a=b
end
Until now, didn’t actually need fortran on my mac, having it in office and the home desktop served all needs.
But yesterday really had to install it on mac, so here’s a simple step for installing gfortran on Mac os Mavericks.
There are many ways to install gfortran, but I used the following way and this information is provided to aid my memory and in the hope that others will find it a useful resource.
I choose to use the relevant tar-file provided at the High Performance Computing web page. For convenience and simplicity I choose to install the “gfortran only” distributions.
The files are gzip’d tar-files which should be unpacked by using the Terminal application, as follows:
sudo tar xvfz [tar-file] -C /
That’s it, gfortran is installed at /usr/local/bin location
Enjoy!
The interface block is perhaps the least understood part of fortran syntax. It was introduced in Fortran 90 and serves some very important functions.
Coding in modern fortran, then its worth understanding interface Statement.
INTERFACE
type FUNCTION name(arg-1, arg-2, ..., arg-n)
type, INTENT(IN) :: arg-1
type, INTENT(IN) :: arg-2
..........
type, INTENT(IN) :: arg-n
END FUNCTION name
....... other functions .......
END INTERFACE
The INTERFACE statement is the first statement in an interface block. The interface block is a powerful structure that was introduced in FORTRAN 90. When used, it gives a calling procedure the full knowledge of the types and characteristics of the dummy arguments that are used inside of the procedure that it references. This can be a very good thing as it provides a way to execute some safety checks when compiling the program. Because the main program knows what argument types should be sent to the referenced procedure, it can check to see whether or not this is the case. If not, the compiler will return an error message when you attempt to compile the program. Also with just a slight variation, it can then be used to create generic functions for use inside of a program unit.
Another additional function that the interface statement can perform is to provide a way to extend the performance of the assignment operator ( the equals sign ). Our general conception of this operator is that it just stores some value into a memory location associated with some variable. However when you think about it, the assignment operator will actually perform operations more complicated than that. For instance, when a real number is truncated to be stored as an integer or when spaces are added on to a character string to completely fill out the its allotted memory space. The interface statement can extend such abilities to instances that involve your own derived types.
The last function that the interface statement can preform is the extension of FORTRAN’s intrinsic operators ( +, -, .lt., .eq., etc. ) or the ability to create your own operators. For example, you could create your own operator called .pi. that you could use to transform angles input in degrees to their corresponding value in radians. Often times, things like that are just as easily carried out using a call to a function or a subroutine so don’t be concerned about the details. Just remember that the power exists to do such things in FORTRAN.
via this
More info on interface in fortran at the following two links
My daughter’s books are littered with word search puzzles as shown above. I don’t know if she likes them but they are everywhere. From environment studies to English to gk, every book has it share of word search exercises.
While helping her one day on one such puzzle, decided to code a program to create these puzzles. So here’s my attempt in fortran.
Given a list of words, the fortran program creates ‘find the word’ puzzle with all the given words. The output can be printed or exported to excel and solved.
Why not give it a try? Another software slave in action!
MODULE MODWORDS
INTEGER, PARAMETER :: COL=8
INTEGER, PARAMETER :: iwrd=12
Character(len=1) SPACE
character(len=14) :: swords(iwrd),TEMP
CHARACTER(LEN=1) :: BOARD(COL,COL)
CHARACTER(LEN=1) :: USED(26)
CHARACTER(LEN=26) :: CHARS
INTEGER :: IUSED
LOGICAL :: LBOARD(COL,COL)
data space/" "/
data swords/"bat","anmol","cat","song","carpenter","school","book","milk"/
END MODULE
program words
USE MODWORDS
implicit none
integer i
CALL INIT
CALL SHOW
PRINT *, USED(1:IUSED-1)
end program
SUBROUTINE SHOW
USE MODWORDS
DO I=1,COL
DO J = 1,COL
WRITE(*,'(A2)', ADVANCE='NO') BOARD(J,I)
END DO
WRITE(*,*)
END DO
WRITE(*,'(4(A14))') (SWORDS)
!CALL SHOW2
END
SUBROUTINE SHOW2
USE MODWORDS
DO I=1,COL
DO J = 1,COL
WRITE(*,'(L2)', ADVANCE='NO') LBOARD(J,I)
END DO
WRITE(*,*)
END DO
END
SUBROUTINE INIT
USE MODWORDS
CHARS='ABCDEFGHIJKLMNOPQRSTUVWXYZ'
USED=SPACE
IUSED=1
LBOARD =.FALSE.
! WLEN=LEN_TRIM(SWORDS(1))
BOARD=SPACE
!call show
CALL FILLBOARD
END
FUNCTION ISUSEDORNOT(CHAR) RESULT(STER)
USE MODWORDS
LOGICAL STER
CHARACTER(LEN=1) CHAR
STER=.FALSE.
DO I=1,IUSED
IF(USED(I) == CHAR) THEN
STER =.TRUE.
EXIT
END IF
END DO
END FUNCTION
SUBROUTINE FILLBOARD
USE MODWORDS
REAL A,B(2)
INTEGER WLEN,X,Y
LOGICAL CHECKOK,SS,ISUSEDORNOT
!XXXXXXXXXX
!Y 4 3 2
!Y 5 % 1
!Y 6 7 8
!
CALL random_seed
DO I=1,iwrd
WLEN=LEN_TRIM(SWORDS(I))
do2: do
do1: DO
call random_number(B)
X=INT(B(1)*COL)+1
Y=INT(B(2)*COL)+1
call random_number(A)
NUM= INT(A*8)+1
! print*, x,y,num,wlen
IF(.NOT. LBOARD(X,Y)) EXIT do1
END DO do1
! print *,x,y,wlen,num
SS= CHECKOK(X,Y,WLEN,NUM,I)
IF(SS) THEN
LBOARD(X,Y) =.TRUE.
TEMP=SWORDS(I)
BOARD(X,Y)=TEMP(1:1)
IF(.NOT. ISUSEDORNOT(TEMP(1:1))) THEN
USED(IUSED)=TEMP(1:1)
IUSED=IUSED+1
END IF
! print *, temp
CALL FILLB(X,Y,WLEN,I,NUM)
! call show
! print *, swords(i),i,num
! read *,a
exit do2
END IF
end do do2
END DO
DO I=1,COL
DO J=1,COL
IF(BOARD(I,J) == SPACE) THEN
call random_number(A)
NUM=INT(A*(IUSED-1))+1
BOARD(I,J)=USED(NUM)
LBOARD(I,J)=.TRUE.
END IF
END DO
END DO
END
SUBROUTINE FILLB(X,Y,WLEN,II,IDIR)
USE MODWORDS
INTEGER X,Y,WLEN,II,IDIR
LOGICAL ISUSEDORNOT
SELECT CASE(IDIR)
CASE (1)
XLIM=1
YLIM=0
CASE (2)
XLIM=1
YLIM=-1
CASE (3)
XLIM=0
YLIM=-1
CASE (4)
XLIM=-1
YLIM=-1
CASE (5)
XLIM=-1
YLIM=0
CASE (6)
XLIM=-1
YLIM=1
CASE (7)
XLIM=0
YLIM=1
CASE (8)
XLIM=1
YLIM=1
END SELECT
DO I=2,WLEN
X=X+XLIM
Y=Y+YLIM
BOARD(X,Y)=TEMP(I:I)
LBOARD(X,Y) =.TRUE.
IF(.NOT. ISUSEDORNOT(TEMP(I:I))) THEN
USED(IUSED)=TEMP(I:I)
IUSED=IUSED+1
END IF
END DO
END
FUNCTION CHECKOK(X,Y,WLEN,IDIR,IWORD) RESULT(SRET)
USE MODWORDS
LOGICAL SRET
INTEGER X,Y,WLEN,IDIR,IWORD
INTEGER XLIM,YLIM,X1,Y1,XLIM1,YLIM1,I
! 1(1,0) 2(1,-1) 3(0,-1) 4(-1,-1)
! 5(-1,0) 6(-1,1) 7(0,1) 8(1,1)
SRET=.TRUE.
SELECT CASE(IDIR)
CASE (1)
XLIM=1
YLIM=0
CASE (2)
XLIM=1
YLIM=-1
CASE (3)
XLIM=0
YLIM=-1
CASE (4)
XLIM=-1
YLIM=-1
CASE (5)
XLIM=-1
YLIM=0
CASE (6)
XLIM=-1
YLIM=1
CASE (7)
XLIM=0
YLIM=1
CASE (8)
XLIM=1
YLIM=1
END SELECT
XLIM1=XLIM
YLIM1=YLIM
XLIM=X+XLIM*WLEN
YLIM=Y+YLIM*WLEN
! print *, 'xlim ',xlim,ylim,sret
IF(((XLIM .GT. COL) .or. (XLIM .LE. 0)) .or. &amp;
((YLIM .GT. COL) .or. (YLIM .LE. 0))) THEN
SRET=.FALSE.
RETURN
END IF
IF(.NOT. SRET) RETURN
TEMP=SWORDS(IWORD)
X1=X
Y1=Y
! print *,'x y ',x,y,wlen,idir,xlim1,ylim1,sret
DO I=2,WLEN
X1=X1+XLIM1
Y1=Y1+YLIM1
! IF((TEMP(I:I) == BOARD(X1,Y1)) .OR. LBOARD(X1,Y1)) THEN
IF(LBOARD(X1,Y1)) THEN
SRET=.FALSE.
! print *, 'heelo ',temp,x1,y1
RETURN
END IF
END DO
END FUNCTION
Someone at Stackoverflow had this following question on Fortran. I’ve recently come to maintain a large amount of scientific calculation-intensive FORTRAN code. I’m having difficulties getting a handle on all of the, say, nuances, of a forty year old language, despite google & two introductory level books. The code is rife with “performance enhancing improvements”. Does anyone have any guides or practical advice for de-optimizing FORTRAN into CS 101 levels? I have spent two years, optimising a mammoth fortran IV code to modern fortran 90 and so my advice is similar to this reply. Common FORTRAN-isms I deal with, that hurt readability are:
Strategies for solving these involve:
Check out other fortran related posts here
Slowing down is the key to increased speed.
Past couple of months I was dabbing with fortran GUI and trying pgplot graphics library. I have produced gui’s in c, vb and then integrated them with fortran, but creating GUIs from fortran was new to me.
As the exploration began I took the fire aim adjust approach!! Dived deep into the tutorials and anything that I could lay my hands on.
Quickly from tutorials I graduated to actually creating my own little programs. This went on for a couple of months.
In the beginning I was sprinting as hard as possible. Learning, doing, getting stuck, reading and then doing again. The pace was fast.
But as I become comfortable, my approach shifted. I slowed.
I wrote a program and pondered how and what am I actually doing. This slowing down and pondering doubled my learning. It felt like I was learning at greater pace with this slowdown.
So the technique I want to advocate to anyone learning a new programming language, a new analysis tool or cad software, is to sprint in the first few weeks. Race and learn as much as you are able to handle. Dive deep and continue the pace as long as you are able to.
When exhaustion, sense of acheivement begins to creep in, slow down. Become deliberate in what you do? Question why and what you are doing?
I hope applying this method will help you as much as it has helped me.
What are your views, do let it out in the comments.
SukhbinderSingh.com 