Patran Neutral Mesh File Writer updated

Sphere

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. 🙂

 

Few Python Regex for FORTRAN Code

 

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) 

 

Simple Python code to Extract Fortran Routines and its Argument

fortranRoutinesWithPython
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

Lights Out Game In Fortran

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.

LightsOut_game_in_fortran

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


Fortran from python – Python Binding for Fortran Code

fortran from python
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 Simple Fortran code to Simulate Stochastic Stock Price Movement

Stochastic_Stock_fortran_ProgramA 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.

Knight in Fortran

Knight Move Game in Fortran 90
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..

Continue reading

The development of FORTRAN I

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

Calling Fortran From C Sharp

Fortran DLL and ExcelIf 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++.

Few Toy Programs to Mesh Excel and Fortran together

Fortran DLL and Excel
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.

Continue reading

Patran Neutral Mesh File Writer

Patran Neutral Mesh File Writer
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.

Continue reading

New Fortran STL Binary File Reader

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=&quot;/Users/Sukhbinder/Desktop/projects/projects/work/stlreader/zDino/Dino.stl&quot;
       end if

       iunit=13
       print *, filname
       open(unit=iunit,file=filname,status='old',access='stream', &amp;
     &amp; 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=&quot;/Users/Sukhbinder/Desktop/projects/projects/work/stlreader/zDino/Dino.stl&quot;
      end if

      iunit=13
     
      open(unit=iunit,file=filename,status='old',form='unformatted', &amp;
    &amp; access='direct',recl=80)

      read(iunit,rec=1)title
      close(iunit,status='keep')

      open(unit=iunit,file=filename,status='old',form='unformatted', &amp;
    &amp; access='direct',recl=4)
     
      read(iunit,rec=21)ntri
      close(iunit,status='keep')
           
      open(unit=iunit,file=filename,status='old',form='unformatted', &amp;
    &amp; 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(*,*) &quot;normals: &quot;,n
         write(*,*) &quot;v1 : &quot;,x1
         write(*,*) &quot;v2 : &quot;,x2
         write(*,*) &quot;v3 : &quot;,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

Installing gfortran on Mac OS X Mavericks

installing gfortran on Maverick
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!

Interface Block in Fortran

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

Software Slave: Fortran Program to Create ‘Find the Word’ Puzzle

20131104-213329.jpg

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;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


Few Strategies for Refactoring old Fortran code to new Fortran

Fortran 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:

  • Common blocks
  • Implicit variables
  • Two or three DO loops with shared CONTINUE statements
  • GOTO’s in place of DO loops
  • Arithmetic IF statements
  • Computed GOTO’s
  • Equivalence REAL/INTEGER/other in some common block

Strategies for solving these involve:

  1. Get Spag / plusFORT, worth the money, it solves a lot of them automatically and Bug Free(tm)
  2. Move to Fortran 90 if at all possible, if not move to free-format Fortran 77
  3. Add IMPLICIT NONE to each subroutine and then fix every compile error, time consuming but ultimately necessary, some programs can do this for you automatically (or you can script it)
  4. Moving all COMMON blocks to MODULEs, low hanging fruit, worth it
  5. Convert arithmetic IF statements to IF..ELSEIF..ELSE blocks
  6. Convert computed GOTOs to SELECT CASE blocks
  7. Convert all DO loops to the newer named F90 syntax
  8. Convert equivalenced common block members to either ALLOCATABLE memory allocated in a module, or to their true character routines if it is Hollerith being stored in a REAL

Check out other fortran related posts here

Slowing Down is the Key…..

slow down

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.