Quasi Engineer

100 Days of Code Log R1D15

2021.10.12 2021.10.12 ⌨️ Coding/➿ 100DaysofCode 385 2 mins

Randomly Generated plot in Fortran

I owe this to one of the participants of matlab contest - Matlab Mini Hack - where this simple but powerful method was shared. Here is my implementation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
module charchal
    use, intrinsic :: iso_fortran_env, only : wp=>real64
    implicit none
    private
    real, parameter :: pi = 4.*atan(1.)
    public :: plot_charchal

contains
subroutine plot_charchal
    integer(wp), parameter :: m = 2e5
    real(wp) :: a(2,2), b(2), c(2,2), d(2), e(2,2), f(2), x(2), n, data(2, m)
    integer(wp) :: i, clr(m), fn
    call randinit(a)
    call randinit(c)
    call randinit(e)

    b = [0.0, 2*rand()]
    d = [0.0, rand()]
    f = [0.0, rand()]
    x = 0
    clr = 0

    do i = 1, m
        call random_number(n)
        if (n <= 0.85) then
            x = matmul(a, x) + b
            clr(i) = randint(3)
        else if (n <= 0.94) then
            x = matmul(c, x) + d
            clr(i) = randint(3) + 1
        else if (n <= 1) then
            x = matmul(e, x) + f
            clr(i) = 2 +randint(4)
        end if
        data(:, i) = x
!         clr(i) = randint(3)
    end do

    open(newunit=fn, file='charchal.dat')
    do i= 1, m
        write(fn, *) data(:,i), clr(i)
    end do
    close(fn)

    call execute_command_line('gnuplot -p charchal.plt')
end subroutine plot_charchal

subroutine randinit(a)
    real(wp), intent(inout) :: a(2,2)
    call random_number(a)
    a = 2*a - 1
end subroutine randinit

function randint(seed) result(val)
    integer, intent(in) :: seed
    integer :: val
    val = floor(seed*(1+rand()))
end function randint

end module charchal

program main

    use charchal, only : plot_charchal
    implicit none
    call plot_charchal()

end program main

And the following goes into the charchal.plt file.

100 Days of Code Log R1D14

2021.10.10 2021.10.10 ⌨️ Coding/➿ 100DaysofCode 184 1 min

Sieve of Eratosthenes

Here’s the fortran implementation of Sieve of Eratosthenes, an ancient method for finding prime upto n numbers (and one of the efficient method for finding small primes).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
program sieve_of_eratosthenes
    use, intrinsic :: iso_fortran_env, only : error_unit
    implicit none

    call prime_sieve_eratosthenes(170)

    contains
    subroutine prime_sieve_eratosthenes(n)
        integer, intent(in) :: n
        logical :: is_prime(2:n) ! Boolean for knowing whether number is prime or not.
        integer :: p, i

        if(n < 1) then
            write(error_unit, *) 'n must be > 1'
            stop 1
        end if

        is_prime = .true.
        p = 2

        do while (p**2 <= n)
            if (is_prime(p)) then
                forall (i = p**2:n:p)
                    is_prime(i) = .false.
                end forall
            end if
            p = p + 1
        end do

        write (*, '("The primes upto ",i0," are: ")', advance='no') n
        do p = 2,n
            if ( is_prime(p) ) write(*, '(i0, ", ")', advance='no') p
        end do
        write (*,*)
    end subroutine prime_sieve_eratosthenes
end program sieve_of_eratosthenes

100 Days of Code Log R1D13

2021.10.10 2021.10.10 ⌨️ Coding/➿ 100DaysofCode 20 1 min

Today, I’ve forked All Algorithms archive for Fortran and have implemented few algorithms. Today’s highlight has been gauss elimination method!

Gauss Elimination in Fortran

2021.10.10 2021.10.10 🗃️ Tech/🧰 Engg 479 3 mins

Gauss Elimination is a well known method for solving system of linear equations. It has time complexity of order $O(n^3)$ and hence, it fares better than cramer’s rule or Gauss-Jordan method.

The following is my implementation of Gauss elimination method in fortran. I’ve made use of scaled pivoting in order to reduce round-off errors.

100 Days of Code Log R1D12

2021.10.9 2021.10.9 ⌨️ Coding/➿ 100DaysofCode 431 3 mins

It’s been a great experience today, solving the tower of hanoi. Though a simple solution, the implementation was hard especially since I made some mistake in variable which took more than an hour (and resorting to gdb) to debug. So, the victory feels sweet. Here’s the code. Will try to make it more efficient and make a blog post out of it tomorrow.

100 Days of Code Log R1D11

2021.10.7 2021.10.9 ⌨️ Coding/➿ 100DaysofCode 251 2 mins

Today’s venture was rather short. I reimplemented the logistic map with elemental function. Here’s the code.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
program logisticmap
    use, intrinsic :: iso_fortran_env, only : dp=>real64, error_unit
    implicit integer(i-k)
    integer, parameter :: l=2001, b = 2001
    real(dp) :: r(l), x0(b), x(l, b)
    integer :: fileno, istat, n
    character(len=1024) :: msg
    print *, "No of iterations: "
    read *, n
    x0 = linspace(0._dp, 1._dp, b)
    r = linspace(3.5_dp, 4._dp, l)

    forall (i=1:l)
        x(i,:) = lm(x0, r(i), n)
    end forall
    print *, "Calculations done"

    open(newunit=fileno, file='logisticmap.dat', iostat=istat, iomsg=msg)
    if (istat /= 0) then
        write(error_unit, '(*(A))') trim(msg)
        stop 1
    end if
    do i=1,l
        write(fileno, *) r(i), x(i, :)
    end do
    close(fileno)
    call execute_command_line('gnuplot -p logisticmap.plt')

    contains
    elemental function lm(x0, r, n) result(xn)
        integer, intent(in) :: n
        real(dp), intent(in) :: x0, r
        real(dp) :: xn
        integer :: i

        xn = x0
        do i = 1, n
            xn = r*xn*(1-xn)
        end do
    end function lm

    function linspace(x0, xn, n) result(x)
        integer, intent(in) :: n
        real(dp), intent(in) :: x0, xn
        real(dp) :: x(n), step
        step = (xn-x0)/(n-1)
        do concurrent(i=1:n)
            x(i) = x0 + step*(i-1)
        end do
    end function linspace
end program logisticmap

And here’s the plot generated.

The World as I See It - Albert Einstein

2021.10.7 2021.10.9 💭 Musings/📝 Write-ups 1276 6 mins

Albert Einstein is the scientific genius that the world worships whether they understand his works or not. The change in perspective he brought is so mind boggling that it’s difficult to believe that he is actually a human. But, he is truly a human and the book The World as I see it is a proof of that. A collection of various letters and speeches made by Albert Einstein, this book adds a whole new dimension to whom we know as a frizzly haired genius!

100 Days of Code Log R1D10

2021.10.6 2021.10.7 ⌨️ Coding/➿ 100DaysofCode 177 1 min

Today’s learning can be split mainly into the following two sections.

I/O Formatting

The formatting string in fortran is passed as a string to print or write statements.

1
2
3
4
character(len=10) :: fmt_str
fmt_str = '(*(i3))'
print fmt_str, 900
write (*, fmt_str) 900

The following formats are available

100 Days of Code Log R1D9

2021.10.5 2021.10.7 ⌨️ Coding/➿ 100DaysofCode 772 4 mins

Continuing with the lessons on fortran, the following stuff were covered today.

Procedures in detail

Intents:

A fortran function’s parameters may have three possible intents

in
out
inout

1
2
3
4
5
6
function one(x) result (y)
	real, intent(inout) :: x(:)
	logical :: y
	x = 1.
	y = .true.
end function one

Note that the result variable doesn’t (and cannot) have intent.

100 Days of Code Log R1D8

2021.10.4 2021.10.4 ⌨️ Coding/➿ 100DaysofCode 415 2 mins

Today has been the beginning of Week 2 in Fortran for scientific Computing course. This had the introduction to arrays, mainly, along with elemental procedures. A word or two were mentioned about “pure functions” but those remain a mystery for now. But, here are some highlights.