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

 

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.

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

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

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

How to Create PGPLOT Fortran Programs With Intel Fortran Compiler ? An animated Journey…

How to create PGPLOT P\programs with Intel Compiler?

How to create PGPLOT programs with Intel Compiler?

How to create PGPLOT Fortran programs with Intel compiler is one of the popular post on this blog and every few days I get a request to share the powerpoint, so here’s the animated version of the powerpoint for quick reference..

Using PGPLOT, then you should visit and see these pages. Want the make file shown in the above animation, visit  How to use pgplot in windows with intel fortran. to download it!

Plants Become Weeds when they…..

FortranDLL in visual studio

How to Create Fortran DLL in Visual Studio? – The simple way


Plants become weeds when they obstruct our plans or our tidy maps of the world. Fortran is also perceived in the same way. Although still a formidable programming language, many now treat it as weeds.

They encounter it and they go back to fortran 66 mindset. Fortran’s new capabilities are discounted and people look past the new features of the language to talk about the horrible “go to” mess!! They might want to look at this Myths Of Fortran post!

There’s so much legacy code in fortran, but no one wants to touch it. The easiest solution commonly used is to convert it into a dll or some library and use with new languages. This explains the popularity of the post how to create fortran dll in visual studio on this blog.

This brings me to today’s post. Here’s an even more convenient animated gif to explain the process of converting a fortran program to a dll in visual studio with intel fortran compiler!

Program to extract results from Substructure Analysis!

ANSYS substructure analysis
Substructuring is a technique in finite element analysis for solving mechanical engineering problems involving large structures which are too enormous for existing computer hardware and software to handle.

This technique involves breaking the system into parts, analyzing the parts, and then re-assembling the total system using selected degrees of freedom from the parts.

Learn how to do sub structure analysis in ansys with this excellent PowerPoint titled The Art Of Matrix Reduction

Well the fun starts only after this, when you need to access the mass and stiffness matrices. I wrote a small code last year just to do that.

Take a look.
Program to extract Mass and Stiffness Matrix from Ansys Sub Structure analysis result file

Automatic Differentiation and Dual Numbers

dual numbers fortran In my book, its a great week if you learn something new on the first day of the week. A few weeks back I stumbled upon automatic differentiation and dual numbers.

Fascinated I decided to explore these two topics and actually coded this Fortran module to do some basic operations with dual numbers.

What are dual numbers?
They are like complex numbers with a twist.

Read more about them here at this following links

  1. Dual Numbers SImple Math.
  2. Autodiff in Perl

Meanwhile here’s the Fortran module for using dual numbers. Still in early development. I hope to take this up latter when I have some more time.

Fortran Module For Dual Numbers

Meshing Fortran and Matlab together with Mex

matlab fortran A few days back colleague of mine needed a help in calling Fortran from matlab. At that time I knew about the call system command.

But as we tried solving his problem, I learnt about Mex. I knew Mex was matlab executable but had never explored them or used or created them.

My colleague’s problem was solved once we got matlab and the c compiler. That day and few days after, I continued exploring the dark alleys of Mex.

Learnt a lot.  Here’s the list of links that were great at explaining the meshing of fortran and c with matlab with Mex.

Perhaps someday I will use them in some project, but till then let them be kept here.

Meet fable and fortwarp

I met fable and fortwarp a few weeks back. I am beginning to know them and before I get too comfy with them, I thought of introducing them to the readers of this blog.

If you work with c,c++ and fortran, you will definitely like their company!

Fable
fable converts fixed-format Fortran sources to C++. The generated C++ code is designed to be human-readable and suitable for further development. In many cases it can be compiled and run without manual intervention.

fable comes with the C++ fem Fortran EMulation library.

fable is written in Python (version 2.x)

The name “fable” is short for “Fortran ABLEitung”. “Ableitung” is a german word which can mean both “derivative” and “branching off”.

FortWarp
FortWrap is a python script that parses Fortran 90/95/200X source files and generates wrapper code for interfacing with the original Fortran code from C++.

FortWrap is compatible with both g95
and gfortran and is intended to be used with Fortran code that takes an object oriented approach and makes use of Fortran derived types.

The resulting wrapper code provides a C++ interface that wraps the Fortran derived types with C++ proxy classes.

Top 10 Algorithms of 20th century!

As i was reading this old PDF mentioning the top 10 algorithms of 20th century, was pleasantly surprised to see creation of fortran compiler listed as one of them!

The creation of Fortran may rank as the single most important event in the history of computer programming: Finally, scientists (and others) could tell the computer what they wanted it to do, without having to descend into the netherworld of machine code. Although modest by modern compiler standards—Fortran I consisted of a mere 23,500 assembly-language instructions—the early compiler was nonetheless capable of surprisingly sophisticated computations. As Backus himself recalls in a recent history of Fortran I, II, and III, published in 1998 in the IEEE Annals of the History of Computing, the compiler “produced code of such efficiency that its output would startle the programmers who studied it.”

Interestingly being in aeronautical/aerospace industry, I have encountered most of the algorithms listed there. Do read the PDF report and see how many of them you have encountered?