Nugroho's blog.: fortran
Showing posts with label fortran. Show all posts
Showing posts with label fortran. Show all posts

Monday, December 22, 2014

Compare Native Loop Time in Python with "homemade" Fortran Module

This code print d and e as result of two matrix addition, e's using python native code, d's using fortran module compiled with F2PY

The code
import numpy as np
import aravir as ar
import time

n = 1000

u = np.ones((n,n))
v = np.ones((n,n))
e = np.ones((n,n))

t = time.clock()
d = ar.add3(u,v)
tfortran= time.clock()-t

t = time.clock()
for i in range (n):
for j in range (n):
e[i,j] = u[i,j]+v[i,j]
tnative = time.clock()-t

print 'fortran ', d
print 'native', e
print 'tfortran = ', tfortran, ', tnative = ', tnative


The fortran module I imported to python
        subroutine add3(a, b, c, n)
double precision a(n,n)
double precision b(n,n)
double precision c(n,n)

integer n
cf2py intent(in) :: a,b
cf2py intent(out) :: c,d
cf2py intent(hide) :: n
do 1700 i=1, n
do 1600 j=1, n
c(i,j) = a(i,j)
$ +b(i,j)

1600 continue
1700 continue
end

save it as aravir.f and compile using
$ f2py -c aravir.f -m aravir

And here the result
$ python cobamodul.py 
fortran [[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]
...,
[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]]
native [[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]
...,
[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]
[ 2. 2. 2. ..., 2. 2. 2.]]
tfortran = 0.069974 , tnative = 1.202547

The Desktop, :)


Friday, December 19, 2014

Using 'Home-Made' Fortran Binary as Python module

Python is easy to use, but it's slow, especially for loop computation. So I compute it using fortran like this

        subroutine subs(a, b, n)
double precision a(n)
double precision b(n)
integer n
cf2py intent(in) :: a
cf2py intent(out) :: b
cf2py intent(hide) :: n
! b(1) = a(1)
do 100 i=2, n
b(i) = a(i)-1
100 continue
end


save it as aravir.py and do the following command
$ f2py -c aravir.f -m aravir

To use the module on the python I use the code below
import numpy as np
import aravir as ar

a = np.linspace(0,1,100)

b = ar.subs(a)

print a
print b

:)


Thursday, November 7, 2013

Laggy Mavericks and QuickLook that Neither Quick nor Look on Anything

Don't know what happened, but my Maverick is getting slower and slower.

QuickLook can't preview .f,  .py and .tex files (and other files type too as I tried to "quick-looking" other type but fortran, python and latex  is what I "quick-looking" most).

Thursday, May 2, 2013

Monte Carlo PI

 
c calculating pi by throwing stones
c
PROGRAM stones
IMPLICIT none
c
c declarations
REAL*8 area, x, y, DRAND48
INTEGER i, max, pi, seed
c
c set parameters (number of stones, seed for generator)
max = 1000
seed = 11168
c
c open file, set initial value, seed generator
OPEN(6, FILE='pif.dat')
pi=0
call seed48(seed)
c
c execute
DO 10 i=1, max
x = DRAND48()*2-1
y = DRAND48()*2-1
IF ((x*x + y*y) .LT. 1) THEN
pi = pi+1
ENDIF
area = 4.0 * pi/REAL(i)
WRITE(6,*) i, area
10 CONTINUE
STOP
END

PI on MPI

still trying...
 program main
include "mpif.h"
double precision PI25DT
parameter (PI25DT = 3.141592653589793238462643d0)
double precision mypi, pi, h, sum, x, f, a
integer n, myid, numprocs, i, ierr
c function to integrate
f(a) = 4.d0 / (1.d0 + a*a)

call MPI_INIT(ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)


10 if ( myid .eq. 0 ) then
print *, 'Enter the number of intervals: (0 quits) '
read(*,*) n
endif
c broadcast n
call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
c check for quit signal
if ( n .le. 0 ) goto 30
c calculate the interval size
h = 1.0d0/n
sum = 0.0d0
do 20 i = myid+1, n, numprocs
x = h * (dble(i) - 0.5d0)
sum = sum + f(x)
20 continue
mypi = h * sum
c collect all the partial sums
call MPI_REDUCE(mypi,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,
& MPI_COMM_WORLD,ierr)
c node 0 prints the answer.
if (myid .eq. 0) then
print *, 'pi is ', pi, ' Error is', abs(pi - PI25DT)
endif
goto 10
30 call MPI_FINALIZE(ierr)
stop
end

Wednesday, November 19, 2008

First on Fortran

My first program written in fortran is below


!================================================
PROGRAM Grafik

IMPLICIT NONE
!declare variable
REAL :: gaya,ygaya,k,ky,a,ay,m,v,vy,t,x,y,z
REAL :: sbx,sby,sbz,g,fx,fy,fz
REAL :: r,dt,dx,dy,dz,dv,dvy,theta,phi
INTEGER :: i,j,n
!initial value
r =100
g =10
n =10000
k =10000
ky =1000
t=0
dt=0.001
!dx=0.01
m=10
v=0
x=10
y=10
vy=0
!create new file named data.dat (you may change the extension. I used that as I think .dat is cool extension )
OPEN (80, FILE = 'data.dat', STATUS = 'NEW', FORM = 'FORMATTED')
WRITE(80,*) ' i t x y'
CLOSE(80)


OPEN(80, ACCESS = 'APPEND', FILE = 'data.dat', FORM = 'FORMATTED')
!the algorithm
DO i=1,n
gaya=-k*x
a=gaya/m
dv=a*dt
v=v+dv
dx=v*dt
x=x+dx

ygaya=-ky*y
ay=ygaya/m
dvy=ay*dt
vy=vy+dvy
dy=vy*dt
y=y+dy
t=t+dt

!write the result of calclation to data.dat
WRITE(80,81) i, t, x, y
81 FORMAT(i7,' 'e14.7,' ',e14.7,' ',e14.7)
END DO

CLOSE(80)


END PROGRAM Grafik
!===========================================

save to grafik.f90 then compile using this syntax
$gfortran grafik.f90

that'll produce file named a.out

run using
$./a.out

and data.dat will appear
to plot the graphic using value in data.dat I use kst

$kst -x 2 -y 3 data dat

mean kst will use 2nd column as x-axis n 3rd as y-axis


Yeah, I know that the steps look boring, so I used sh script below to generate it automatically.
using kwrite, create text file like below:
rm *.dat
gfortran grafik.f90
./a.out
kst -x 2 -y 3 data.dat

save as run.sh, make it executable then run:
$./run.sh




!I also used script below to generate x, y, z in text file and plot it using kst
!
!DO i=1 ,100
! dt=i
! dt=dt/100
! sbx=r*sin(dt)
! sby=r*cos(dt)
! sbz=dt
! WRITE(80,81) sbx, sby, sbz
! 81 FORMAT(F10.2,' ',e14.7,' ',e14.7)
!END DO

!DO i=1, n
!DO j=1, n
! theta=i
! phi=j
! sbx=r*sin(theta)*cos(phi)
! sby=r*sin(theta)*sin(phi)
! sbz=r*cos(theta)
! !sbz=2
! WRITE(80,81) i, j, sbx, sby, sbz
! 81 FORMAT(i7,' ',i7,' 'e14.7,' ',e14.7,' ',e14.7)
! PRINT*, 'theta: ',i,' -------------'
! PRINT*, 'phi: ',j,' -------------'
! PRINT*, 'sb-x :', sbx
! PRINT*, 'sb-y :', sby
! PRINT*, 'sb-z :', sbz
!END DO
!END DO



There is something strange in two script above. The result is x, y, z (THREE variable) but why the plot is just in TWO axis? Err..., that the weakness of kst.

At prior time :), I used gnuplot to generate "3D" picture--of course, it just projection of z-axis (or x-axis if you using left-handed coordinate)--but I have'nt yet found to automating it using .sh script.

If you like to using gnuplot:
$gnuplot
>plot "data.dat" using 2:3 (oops, that still 2D. Plotting data.dat using 2nd and 3rd column)



>splot "data.dat" using 2:3:4 (that is.)


I've skimmed "man gnuplot" but did'nt found how to do it automatically

anybody can help me automating it using gnuplot?


oh yeah, almost forget, I use Linux.
323f (5) amp (1) android (12) apple (7) arduino (18) art (1) assembler (21) astina (4) ATTiny (23) blackberry (4) camera (3) canon (2) cerita (2) computer (106) crazyness (11) debian (1) delphi (39) diary (286) flash (8) fortran (6) freebsd (6) google apps script (8) guitar (2) HTML5 (10) IFTTT (7) Instagram (7) internet (12) iOS (5) iPad (6) iPhone (5) java (1) javascript (1) keynote (2) LaTeX (6) lazarus (1) linux (29) lion (15) mac (28) macbook air (8) macbook pro (3) macOS (1) Math (3) mathematica (1) maverick (6) mazda (4) microcontroler (35) mountain lion (2) music (37) netbook (1) nugnux (6) os x (36) php (1) Physicist (29) Picture (3) programming (189) Python (109) S2 (13) software (7) Soliloquy (125) Ubuntu (5) unix (4) Video (8) wayang (3) yosemite (3)