# Fortran code snippets

**URL:** <https://fortran-lang.discourse.group/t/fortran-code-snippets/2150>\
**Category:** Uncategorized\
**Created:** [October 24, 2021, 1:39pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150 "2021-10-24T13:39:39Z")\
**Posts on this page:** 20\
**Page:** 3

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [December 25, 2023, 5:14pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/41 "2023-12-25T17:14:29Z")

</div>

Using the long long int and quadruple precision real of gfortran, you can calculate the 100th Fibonacci number using the analytic [Binet’s Formula](https://artofproblemsolving.com/wiki/index.php/Binet%27s_Formula):

```auto
module m
implicit none
integer, parameter :: ikind = selected_int_kind(38) ! long long int
integer, parameter :: rkind = selected_real_kind(32) ! quadruple precision
contains
elemental function fib_binet(n) result(ifib)
! Binet's formula for nth Fibonacci number
integer, intent(in) :: n
integer(kind=ikind) :: ifib
real(kind=rkind), parameter :: sqrt_5 = sqrt(5.0_rkind)
ifib = nint((((1 + sqrt_5)/2)**n - ((1 - sqrt_5)/2)**n)/sqrt_5, kind=ikind)
end function fib_binet
!
elemental function fib_loop(n) result(ifib)
! Use loop to compute nth Fibonacci number
integer, intent(in) :: n
integer(kind=ikind) :: ifib
integer(kind=ikind) :: i, jfib(n)
if (n < 2) then
   ifib = 1
   return
end if
jfib(:2) = [1, 1]
do i=3,n
   jfib(i) = jfib(i-1) + jfib(i-2)
end do
ifib = jfib(n)
end function fib_loop
end module m

program main
use m, only: fib_binet, fib_loop
implicit none
integer, parameter :: n = 100
print "(a,*(1x, i0))", "n, Fib(n) =", n, fib_binet(n), fib_loop(n)
! n, Fib(n) = 100 354224848179261915075 354224848179261915075
end program main

```

---

<div class="post-metadata">

**Author:** ![art-rasa](https://avatars.discourse-cdn.com/v4/letter/a/90ced4/32.png) [@art-rasa](https://fortran-lang.discourse.group/u/art-rasa)\
**Post date:** [December 29, 2023, 10:19pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/42 "2023-12-29T22:19:57Z")

</div>

I made this (rather inefficient) program for solving the [Gauss circle problem](https://en.wikipedia.org/wiki/Gauss_circle_problem). It finds the number of points inside a circle of given radius.

```auto
program gauss_circle
  implicit none
  real, parameter :: pi = atan(1.0) * 4.0
  real :: radius

  print *, "hello from project gauss_circle"
  do
    write (*, '(a)') 'Enter circle radius (negative number quits):'
    read (*, *) radius
    if (radius < 0.0) then
      exit
    end if
    write (*, '(a,t25,f0.2)') 'Circle radius: ', radius
    write (*, '(a,t25,f0.2)') 'Circle area: ', area(radius)
    write (*, '(a,t25,i0)') 'Points inside circle: ', find_points(radius)
    write (*, *) ''
  end do

contains
  integer function find_points(r) result(pts)
    real, intent(in) :: r
    integer :: x, y
    if ( r < 0.0 ) then
      pts = 0
      return
    end if
    if ( r < 1.0 ) then
      pts = 1
      return
    end if
    pts = 0 ! points inside one quadrant
    do x = 1, ceiling(r)
      do y = 1, ceiling(r)
        if ( inside(r, x, y) ) then
          pts = pts + 1
        end if
      end do
    end do
    pts = 4 * pts
    pts = pts + 4 * floor(r) + 1 ! points on x, y axis + centre point
  end function

  real function area(r)
    real, intent(in) :: r
    area = pi * r ** 2.0
  end function

  logical function inside(r, x, y)
    real, intent(in) :: r
    integer, intent(in) :: x, y
    inside = sqrt(x **2.0 + y** 2.0) <= r
  end function
end program gauss_circle

```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [January 3, 2024, 5:09am UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/43 "2024-01-03T05:09:10Z")

</div>

The new year 2024 is the sum of the cubes of the integers from 2 through 9, and 2025 is the sum of cubes from 1 through 9. Given a set of integers one can write a short recursive function to determine if a specified integer can be written as the sum of a subset of the integers, without replacement. The program below decomposes integers from 1 to 2025 into sums of cubes from 1 through 9 when possible, and it runs in about 78s on my PC with `gfortran -O3 -flto -march=native`. Setting jmax to 10 instead of 9 causes the program to run much more slowly. There is likely a faster approach. John Burkardt has a [SUBSET\_SUM](https://people.math.sc.edu/Burkardt/f_src/subset_sum/subset_sum.html) library for this.

```auto
module subset_sum_mod
implicit none
private
public :: subset_sum
contains
recursive pure function subset_sum(isum, vec) result(subset)
! Return the subset of integers in vec(:) that sum to isum,
! assuming that the elements of vec(:) are positive and
! strictly ascending. Return a zero-size array if such
! a subset does not exist.
integer, intent(in) :: isum
integer, intent(in) :: vec(:)
integer, allocatable :: subset(:)
integer, allocatable :: try_subset(:)
integer :: i,n
n = size(vec)
if (n < 1 .or. isum == 0) then
   allocate (subset(0))
   return
end if
if (any(isum == vec)) then
   subset = [isum]
   return
end if
subset = [integer ::]
do i=1,n
   if (vec(i) > isum) return
   try_subset = [vec(i), subset_sum(isum - vec(i), pack(vec, vec /= vec(i)))]
   if (sum(try_subset) == isum) then
      subset = try_subset
      return
   end if
end do
end function subset_sum
end module subset_sum_mod

program main
use subset_sum_mod, only: subset_sum
implicit none
integer, allocatable :: v(:)
integer :: i, j, sumv
integer, allocatable :: subset(:)
integer, parameter :: jmax = 9
v = [(j**3, j=1,jmax)]
sumv = sum(v)
print*,"sum(v) =",sumv
do i=1,sumv
   subset = subset_sum(i, v)
   if (size(subset) > 0) print "(*(1x,i0))",i,subset
end do
end program main

```

Output

```auto
...
 2016 27 64 125 216 343 512 729
 2017 1 27 64 125 216 343 512 729
 2024 8 27 64 125 216 343 512 729
 2025 1 8 27 64 125 216 343 512 729

```

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [January 3, 2024, 12:27pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/44 "2024-01-03T12:27:41Z")

</div>

It is worth noting that your short program did not discover this solution for N = 2016:

```program
    print *, sum([2,4,6,12]**3)
    end

```

That observation leads to a suspicion that for certain values of N the algorithm used may falsely declare “No solution”!

On a related question, Lander and Parkin in '66 (probably using Fortran-66 ) ran a program on a CDC6600 to disprove an [Euler conjecture](https://en.wikipedia.org/wiki/Euler%27s_sum_of_powers_conjecture). The [journal article](https://community.intel.com/t5/image/serverpage/image-id/11600iD6EE41CC15A31C45/image-size/large/is-moderation-mode/true?v=v2&px=999) was short enough (2 sentences, 4 lines of text and 1 equation) to fit on the back of an envelope.

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [January 3, 2024, 12:36pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/45 "2024-01-03T12:36:39Z")

</div>

It is also missing the famous taxicab number 1729 😇

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [January 3, 2024, 2:58pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/46 "2024-01-03T14:58:20Z")

</div>

Poking around in the nice little program given above by @Beliavsky , I noted that the cause of some solutions having been left undiscovered is that it sets JMAX = 9. Raising JMAX to 13, changing the PRINT statement to  
`print "(*(1x,i0))",i,nint(subset**(1.0/3))`  
and changing the DO loop to range from 2000 to 2030, I find:

```auto
 2003 2 3 5 8 11
 2004 1 2 3 5 8 11
 2007 3 5 7 8 10
 2008 1 3 5 7 8 10
 2009 1 4 6 12
 2010 1 4 6 9 10
 2015 2 3 5 7 8 10
 2016 1 2 3 5 7 8 10
 2017 1 2 4 6 12
 2018 1 2 4 6 9 10
 2023 2 5 6 7 11
 2024 1 2 5 6 7 11
 2025 1 2 3 4 5 6 7 8 9

```

His original program gives this line (easy to overlook)

```auto
252 1 8 27 216

```

which, when multiplied by 8, gives the solution that I noted as missing in my earlier post:

```auto
2016 8 64 216 1728

```

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [January 3, 2024, 10:06pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/47 "2024-01-03T22:06:28Z")

</div>

Here is an alternative program that uses a packed bit representation of which cubes are included in the sum. If the sum contains 1^3, bit-0 is set; if the sum contains 2^3, bit-1 is set, and so on up to bit-12. The packed bit value is the index into an integer array S in which the sum of the selected cubes is stored. Then, for values of N within a selected range (such as 2000 to 2040), the indices of the S entries that fall within the range are noted. These entries are then sorted and output. Because every entry in array S is already known to be the sum of cubes, there is no time spent on attempting to split a test number into the sum of cubes, as is done in some other algorithms.

```auto
! Finds integers N that are cubes and are also equal to the sum of cubes of integers
! Only values of N between YLOW and YHIGH are covered.
! The count of the solutions is expected to be no more than KMAX
!
program cubem
implicit none
integer i, j, n, k, m, i0, i1, NMAX
integer, parameter :: IMAX = 13, KMAX = 100
integer, parameter :: YLOW = 2000, YHIGH = 2040
integer, parameter :: cubes(IMAX) = (/ (i*i*i,i=1,IMAX) /)
integer eul(KMAX), idx(KMAX),sol(KMAX)
logical first

k = 0
NMAX = ishft(1,IMAX)-1
do i = 1, NMAX
   n = 0
   do j = 0, iMAX-1
      if (btest(i,j)) n = n+cubes(j+1)
   end do
   if(n >= YLOW .and. n <= YHIGH)then
      k = k+1
      eul(k) = n
      idx(k) = k
      sol(k) = i
   endif
end do
print *,'Sorting ',k,' items'
call inssort(k,idx,eul)
do i = 1, k
   j = idx(i)
   i0 = sol(j)
   write(*,"(I4,2x,'= sum of cubes of {')",advance='no')eul(j)
   first = .true.
   i1 = 0
   do m = 0,IMAX-1
      if(btest(i0,m))then
         if(first)then
            first = .false.
            write(*,"(i3)",advance='no')(m+1)
         else
            write(*,"(',',i3)",advance='no')(m+1)
         endif
         i1 = i1+cubes(m+1)
      endif
   end do
   print *,'}'
   if(i1 .ne. eul(j))print *,'Error, sum of cubes does not equal ',i1
end do
end program

subroutine inssort(N,idx,A)
implicit none
integer, intent(in) :: N, A(N)
integer, intent(in out) :: idx(N)
integer :: i,j,itmp
!
do i = 2, N
   j = i-1
   do while (j > 0)
      if (A(idx(j)) <= A(idx(j+1))) exit
      itmp = idx(j+1)
      idx(j+1) = idx(j)
      idx(j) = itmp
      j = j - 1
   end do
end do
return
end subroutine inssort

```

The output of the program (run time should be less than a second) should include multiple solutions when they exist.

```auto
 Sorting 29 items
2003 = sum of cubes of { 2, 3, 5, 8, 11}
2004 = sum of cubes of { 1, 2, 3, 5, 8, 11}
2007 = sum of cubes of { 3, 5, 7, 8, 10}
2008 = sum of cubes of { 4, 6, 12}
2008 = sum of cubes of { 1, 3, 5, 7, 8, 10}
2009 = sum of cubes of { 1, 4, 6, 12}
2009 = sum of cubes of { 4, 6, 9, 10}
2010 = sum of cubes of { 1, 4, 6, 9, 10}
2015 = sum of cubes of { 5, 6, 7, 11}
2015 = sum of cubes of { 2, 3, 5, 7, 8, 10}
2016 = sum of cubes of { 2, 4, 6, 12}
2016 = sum of cubes of { 1, 5, 6, 7, 11}
2016 = sum of cubes of { 1, 2, 3, 5, 7, 8, 10}
2016 = sum of cubes of { 3, 4, 5, 6, 7, 8, 9}
2017 = sum of cubes of { 1, 2, 4, 6, 12}
2017 = sum of cubes of { 2, 4, 6, 9, 10}
2017 = sum of cubes of { 1, 3, 4, 5, 6, 7, 8, 9}
2018 = sum of cubes of { 1, 2, 4, 6, 9, 10}
2023 = sum of cubes of { 2, 5, 6, 7, 11}
2024 = sum of cubes of { 1, 2, 5, 6, 7, 11}
2024 = sum of cubes of { 2, 3, 4, 5, 6, 7, 8, 9}
2025 = sum of cubes of { 1, 2, 3, 4, 5, 6, 7, 8, 9}
2032 = sum of cubes of { 4, 5, 8, 11}
2033 = sum of cubes of { 1, 4, 5, 8, 11}
2035 = sum of cubes of { 3, 4, 6, 12}
2036 = sum of cubes of { 1, 3, 4, 6, 12}
2036 = sum of cubes of { 3, 4, 6, 9, 10}
2037 = sum of cubes of { 1, 3, 4, 6, 9, 10}
2040 = sum of cubes of { 2, 4, 5, 8, 11}

```

P.S. Thanks to @JohnCampbell for correcting my initialization of NMAX.

---

<div class="post-metadata">

**Author:** ![JohnCampbell](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@JohnCampbell](https://fortran-lang.discourse.group/u/JohnCampbell)\
**Post date:** [January 4, 2024, 3:36am UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/48 "2024-01-04T03:36:18Z")

</div>

nmax = 2\*\*imax-1 ?

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [January 4, 2024, 5:07am UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/49 "2024-01-04T05:07:27Z")

</div>

Thanks for the correction. I have made the correction for NMAX in the post in which I included the source code.

Please note that the special-purpose bit representation that I use is

N = b\_0.1^3+b\_1.2^3+b\_2.3^3+...+b\_{12}.13^3

in contrast to the familiar binary notation

N = b\_0.2^0+b\_1.2^1+b\_2.2^2+...+b\_{12}.2^{12}

where the coefficients b\_0, b\_1,...,b\_n are all bit-valued – i.e, each of them can only be 0 or 1.

The same program, adapted for 4th powers, gave the following scanty results, which may be of interest to @Beliavsky :

```auto
2002 = sum of { 3, 5, 6 } ^ 4        
2003 = sum of { 1, 3, 5, 6 } ^ 4    
2018 = sum of { 2, 3, 5, 6 } ^ 4    
2019 = sum of { 1, 2, 3, 5, 6 } ^ 4

```

No more such years in this century.

@Beliavsky 's program, modified to search for 4-th powers and the range of the sum restricted to 1900-2100, runs in 0.04 s and produces the same results.

With 5th powers, there are no results for N in the range 1900-2000.

Here is my explanation of why the subset program fails to find more than one solution for a given N, even when such solutions exist. For any given N, the main program calls the recursive function with the statement

```auto
subset = subset_sum(i, v)

```

The recursive function executes an algorithm, and either decides “failed” and returns, or decides “here is the subset that satisfies the requirements”, and program control traces back through the chain of recursive calls. To find more than one solution, the main program would have to issue another invocation of sub\_set\_sum() with the same `i` but a proper subset of `v` rather than the original `v`. As of now, the second argument `v` is always the full set of cubes that were generated in the first executable statement of the main program…

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [January 5, 2024, 2:45pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/50 "2024-01-05T14:45:08Z")

</div>

> [@art-rasa](#):
>
> `inside = sqrt(x **2.0 + y** 2.0) <= r`

These kinds of expressions show up often. In fortran, it is better to write this as

```auto
inside = (x **2 + y** 2) <= r**2
or
inside = (x*x + y*y) <= r*r

```

The reason is that floating point exponentiation is an expensive operation, and a single multiplication is less expensive than a `sqrt()`. An aggressive compiler might recognize this identity and do it this way, but the modified code is just as simple to understand, so it is usually better to hand-optimize this expression. This kind of distance test does occur rather frequently.

---

<div class="post-metadata">

**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\
**Post date:** [January 7, 2024, 1:59pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/51 "2024-01-07T13:59:41Z")

</div>

> [@JohnCampbell](#):
>
> nmax = 2\*\*imax-1 ?

Thanks, and I wish to record here that John sent me a cleaned up version of the program. This entertaining problem enables me to provide another example of the expressive power of modern Fortran, and I cannot resist the temptation to post this very short complete program, which only reports the years within a specified range that fit the stipulation, without reporting the subset of integers whose cubes sum is equal to the year.

```auto
program cubem
  implicit none
  integer :: year, i, nk
  integer :: bits
  integer, parameter :: YLOW = 2000, YHIGH = 2050
  integer, parameter :: IMAX = 12 ! sufficient to make (IMAX+1)^3 >= YHIGH
  integer, parameter :: cubes(IMAX) = [(i*i*i,i=1,IMAX)]
  
  nk = 0
  do bits = 1, ishft (1,IMAX)-1
     year = sum(cubes,mask=btest(bits,[(i-1,i=1,IMAX)]))
     if (year >= YLOW .and. year <= YHIGH) then ! add to list
        nk = nk+1
        print '(i2,i6)',nk,year
    end if
  end do
  print '(A,i4,A)','Found ',nk,' years that are sums of cubes'

end program

```

The output years are not in order, and a year will appear more than once if that year can be composed as the sum of the cubes of more than one subset integers.

By using the ASSOCIATE construct, we can flesh out this program to print the terms in the sum as well, and use the OS sort utility to restore order in the output.

```auto
! Given a range of years between 1 and 13^3-1, list years that may be
! expressed as the sum of cubes of integers with no repetitions allowed/
! Pipe the output of the program through the system SORT utility
program cubem

  implicit none
  integer, parameter :: YLOW = 2000, YHIGH = 2050
  integer, parameter :: IMAX = 12 ! smallest possible I s.t. I^3 >= YHIGH
  integer :: year, i, nk, bits
  integer, parameter :: cubes(IMAX) = [(i*i*i, i = 1, IMAX)]

  nk = 0
  do bits = 1, ishft (1, IMAX)-1
     associate (bmask => btest(bits,[(i-1, i = 1, IMAX)]))
        year = sum(cubes, mask = bmask)
        if (year < YLOW .or. year > YHIGH) cycle
        nk = nk+1
        print '(i6,4x,10i3)',year,pack([(i, i = 1, IMAX)], mask=bmask)
     end associate
  end do
  print '(/,A,i4,A)','Found ',nk,' years that are sums of cubes'

end program

```

That’s not even twenty lines of code!

---

<div class="post-metadata">

**Author:** ![art-rasa](https://avatars.discourse-cdn.com/v4/letter/a/90ced4/32.png) [@art-rasa](https://fortran-lang.discourse.group/u/art-rasa)\
**Post date:** [January 7, 2024, 9:07pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/52 "2024-01-07T21:07:19Z")

</div>

My attempt at using type-bound procedures for modelling an NTC thermistor using the [Steinhart-Hart equations](https://en.wikipedia.org/wiki/Steinhart%E2%80%93Hart_equation).

I found the **associate** construct to be very useful for abbreviating long variable names into more concise ones without introducing extra local variables.

I measured the resistance of my thermistor with a multimeter. I exposed it to water/snow “slush” mixture, room temperature air and boiling water while measuring its resistance.

| Temperature [°C] | Resistance [Ω] |
| --- | --- |
| 0 | 27 650 |
| 22.4 | 10 280 |
| 100 | 980 |

The module calculated its coefficients as:

```auto
A: 0.844666967E-3
B: 0.259054825E-3
C: 0.155988630E-6

```

I only have my multimeter’s K-type thermocouple sensor as my “reference”. I plan to make more measurements with both sensors and compare them to see if the thermistor model works as intended.

Driver program omitted.

```auto
module thermistor_data
    implicit none
    private
    public :: therm
    
    type temp_data
        real :: temp_k(3)
        real :: r_ohm(3)
    end type
    
    type therm
        real :: par_A
        real :: par_B
        real :: par_C
        type(temp_data) :: rt_data
    contains
        procedure calculate_parameters
        procedure get_temperature
    end type
    
contains
    ! Calculate Steinhart-Hart coefficients "A,B,C" from temperature and
    ! resistance data
    subroutine calculate_parameters(self)
        class(therm), intent(inout) :: self
        real :: L(3)
        real :: Y(3)
        real :: g2, g3
        real :: a, b, c
        associate (r => self%rt_data%r_ohm, t => self%rt_data%temp_k)
            L = log(r)
            Y = 1.0 / t
            g2 = (Y(2) - Y(1)) / (L(2) - L(1))
            g3 = (Y(3) - Y(1)) / (L(3) - L(1))
            c = ( (g3 - g2) / (L(3) - L(2)) ) * 1.0/(sum(L))
            b = g2 - c * ( L(1)**2 + L(1)*L(2) + L(2)**2 )
            a = Y(1) - (b + L(1)**2 * c) * L(1)
        end associate
        self%par_A = a
        self%par_B = b
        self%par_C = c
    end subroutine
    
    ! Convert a thermistor resistance value into temperature in deg. K
    real elemental function get_temperature(self, r) result(temp_k)
        class(therm), intent(in) :: self
        integer, intent(in) :: r
        real :: t_inv, ln_r
        ln_r = log(real(r))
        associate (a => self%par_A, b => self%par_B, c => self%par_C)
            t_inv = a + b * ln_r + c * ln_r**3
        end associate
        temp_k = 1.0/t_inv
    end function
end module thermistor_data

```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [January 7, 2024, 11:47pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/53 "2024-01-07T23:47:47Z")

</div>

It is simple to emulate [DataFrame.apply](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.apply.html) of Python, which applies a function along an axis of a DataFrame (matrix). The code below, for `real(kind=dp)`, implements apply for functions that map 1D arrays to scalars and 1D arrays to 1D arrays of the same size.

```auto
module kind_mod
implicit none
private
public :: dp
integer, parameter :: dp = kind(1.0d0)
end module kind_mod

module apply_mod
use kind_mod, only: dp
implicit none
private
public :: apply, mean, running_sum
interface apply
   module procedure apply_1d_0d, apply_1d_1d
end interface apply
!
interface
   pure function f_real_1d_0d(x) result(y)
   ! map 1D array to scalar
   import dp
   implicit none
   real(kind=dp), intent(in) :: x(:)
   real(kind=dp) :: y
   end function f_real_1d_0d
!
   pure function f_real_1d_1d(x) result(y)
   ! map 1D array to 1D array of same size
   import dp
   implicit none
   real(kind=dp), intent(in) :: x(:)
   real(kind=dp) :: y(size(x))
   end function f_real_1d_1d
end interface
contains
!
pure function apply_1d_0d(f,x) result(y)
! apply function that maps vector to scalar to each column of x
procedure(f_real_1d_0d) :: f
real(kind=dp), intent(in) :: x(:,:)
real(kind=dp) :: y(size(x,2))
integer :: i
do i=1,size(x,2)
   y(i) = f(x(:,i))
end do
end function apply_1d_0d
!
pure function apply_1d_1d(f,x) result(y)
! apply function that maps vector to scalar to each column of x
procedure(f_real_1d_1d) :: f
real(kind=dp), intent(in) :: x(:,:)
real(kind=dp) :: y(size(x,1),size(x,2))
integer :: i
do i=1,size(x,2)
   y(:,i) = f(x(:,i))
end do
end function apply_1d_1d
!
pure function mean(x) result(xmean)
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: xmean
xmean = sum(x)/max(1,size(x))
end function mean
!
pure function running_sum(x) result(xsum)
real(kind=dp), intent(in) :: x(:)
real(kind=dp) :: xsum(size(x))
integer :: i, n
n = size(x)
if (n < 1) return
xsum(1) = x(1)
do i=2,n
   xsum(i) = xsum(i-1) + x(i)
end do
end function running_sum
end module apply_mod

program xapply
use kind_mod, only: dp
use apply_mod, only: apply, mean, running_sum
implicit none
integer, parameter :: n1=3, n2=2
real(kind=dp) :: x(n1,n2), xsum(n1,n2)
integer :: i
character (len=*), parameter :: fmt_r="(*(f7.3))"
call random_number(x)
print "('matrix')"
do i=1,n1
   print fmt_r, x(i,:)
end do
print "(/,'column means')"
print fmt_r, apply(mean, x)
print "(/,'check column means')"
print fmt_r, (sum(x(:,i))/n1, i=1,n2)
xsum = apply(running_sum, x)
print "(/,'running sum')"
do i=1,n1
   print fmt_r, xsum(i,:)
end do
print "(/,'check running sum')"
do i=1,n1
   print fmt_r, sum(x(:i,:), dim=1)
end do
end program xapply

```

Sample output:

```auto
matrix
  0.312 0.960
  0.871 0.906
  0.523 0.173

column means
  0.569 0.680

check column means
  0.569 0.680

running sum
  0.312 0.960
  1.183 1.866
  1.706 2.039

check running sum
  0.312 0.960
  1.183 1.866
  1.706 2.039

```

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [January 8, 2024, 7:49am UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/54 "2024-01-08T07:49:54Z")

</div>

Here is a generalisation of the inner product to (compatible) multidimensional arrays:

```auto
real :: A(imax,jmax,kmax), &
        B(kmax,lmax,mmax), &
        AB(imax,jmax,lmax,mmax)
... (fill the arrays)
do j = 1,jmax
  do i = 1,imax
    AB(i,j,:,:) = &
      sum( reshape( [(A(i,j,k) * B(k,:,:), &
        k = 1,kmax), &
        [lmax,mmax,kmax]), dim = 3)
  enddo
enddo

```

It is part of my paper-to-be on the mathematics of arrays. You might even get rid of the two do-loops but that would lead to a much longer statement and of larger “nestedness”.

---

<div class="post-metadata">

**Author:** ![jkd2022](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jkd2022/32/1949_2.png) [@jkd2022](https://fortran-lang.discourse.group/u/jkd2022)\
**Post date:** [January 8, 2024, 10:54am UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/55 "2024-01-08T10:54:23Z")

</div>

Not a Fortran snippet per se, but here’s a little bash script for converting old-style logical operators like .le. to the modern style \<=.

I wrote it because LFortran makes (somewhat insistent 🙂) style suggestions when compiling code, so I converted the entire code in one go and have to say it looks better for it.

```auto
#!/bin/bash

while ["$*" != ""]
do
  echo $1
# lower case
  sed -i 's/\.lt\./ < /g' $1
  sed -i 's/\.le\./ <= /g' $1
  sed -i 's/\.gt\./ > /g' $1
  sed -i 's/\.ge\./ >= /g' $1
  sed -i 's/\.eq\./ == /g' $1
  sed -i 's/\.ne\./ \/= /g' $1
# upper case
  sed -i 's/\.LT\./ < /g' $1
  sed -i 's/\.LE\./ <= /g' $1
  sed -i 's/\.GT\./ > /g' $1
  sed -i 's/\.GE\./ >= /g' $1
  sed -i 's/\.EQ\./ == /g' $1
  sed -i 's/\.NE\./ \/= /g' $1
  shift
done

```

Usage (assuming the script is called ‘fmlop’):

Single file:

`fmlop test.f90`

or for many files:

`ls *.f90 | xargs fmlop`

---

<div class="post-metadata">

**Author:** ![msz59](https://avatars.discourse-cdn.com/v4/letter/m/3d9bf3/32.png) [@msz59](https://fortran-lang.discourse.group/u/msz59)\
**Post date:** [January 8, 2024, 12:09pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/56 "2024-01-08T12:09:22Z")

</div>

You can make it in just one invocation of `sed`:

```bash
#!/bin/bash

while ["$*" != ""]
do
  echo $1
  sed -i \
   -e 's/ *\.[lL][tT]\. */ < /g' \
   -e 's/ *\.[lL][eE]\. */ <= /g' \
   -e 's/ *\.[gG][tT]\. */ > /g' \
   -e 's/ *\.[gG][eE]\. */ >= /g' \
   -e 's/ *\.[eE][qQ]\. */ == /g' \
   -e 's/ *\.[nN][eE]\. */ \/= /g' $1
  shift
done

```

This version strips multiple spaces around the operators and also supports mixed-case like `.gE.`, however weird it might look in the original code.

---

<div class="post-metadata">

**Author:** ![jkd2022](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jkd2022/32/1949_2.png) [@jkd2022](https://fortran-lang.discourse.group/u/jkd2022)\
**Post date:** [January 8, 2024, 12:31pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/57 "2024-01-08T12:31:17Z")

</div>

> This version strips multiple spaces around the operators and also supports mixed-case like `.gE.`, however weird it might look in the original code.

Much better! Thanks.

---

<div class="post-metadata">

**Author:** ![hkvzjal](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/hkvzjal/32/3055_2.png) [@hkvzjal](https://fortran-lang.discourse.group/u/hkvzjal)\
**Post date:** [January 11, 2024, 12:09pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/58 "2024-01-11T12:09:24Z")

</div>

A snippet showing a combination of fpm libraries put together to replicate a [Lorentz system](https://en.wikipedia.org/wiki/Lorenz_system)

```Fortran
program main
  use ogpf, only: gpf
  use rklib_module
  use iso_fortran_env, only: output_unit, wp => real64

  implicit none

  integer,parameter :: n = 3 !! dimension of the system
  real(wp),parameter :: tol = 1.0e-12_wp !! integration tolerance
  real(wp),parameter :: t0 = 0.0_wp !! initial t value
  real(wp),parameter :: dt = 0.01_wp !! initial step size
  integer, parameter :: num_steps = 10000 !! number of steps
  real(wp),parameter :: x0(n) = [0.0_wp,0.1_wp,1.05_wp] !! initial x value

  integer :: incr !! increment number
  real(wp) :: time(0:num_steps)
  real(wp) :: x(n,0:num_steps)
  type(rkf45_class) :: prop !! rk integrator
  type(gpf) :: gp !! gnuplot interface
  ! ---------------------------------------------
  ! Initialize buffers
  time = [(t0+incr*dt,incr=0,num_steps)]
  x(:,0) = x0; x(:,1:) = 0._wp
  ! ---------------------------------------------
  ! Initialize method
  call prop%initialize(n=n,f=fvpol,rtol=[tol],atol=[tol])

  do incr = 1, num_steps
    call prop%integrate(time(incr-1),x(:,incr-1),dt,time(incr),x(:,incr))
  end do
  ! static 2D multiplot
  call gp%multiplot(2,2)
  call gp%plot(x(1,:), x(2,:), 'with lines lt 5 lc rgb "blue"')
  call gp%plot(x(1,:), x(3,:), 'with lines lt 5 lc rgb "green"')
  call gp%plot(x(2,:), x(3,:), 'with lines lt 5 lc rgb "red"')
  call gp%run_script()
contains

  subroutine fvpol(me,t,x,f)
    !! Right-hand side of Lorentz system

    class(rk_class),intent(inout) :: me
    real(wp),intent(in) :: t
    real(wp),intent(in) :: x(:)
    real(wp),intent(out) :: f(:)
    real(wp),parameter :: s=10._wp, r=28._wp, b=2.667_wp

    f(1) = s*(x(2) - x(1))
    f(2) = r*x(1) - x(2) - x(1)*x(3)
    f(3) = x(1)*x(2) - b*x(3)

  end subroutine fvpol
end program main

```

Which produces:

 ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/5/5e51b2183bc10d8a348630031fa8446907faea74.png)

To get here:

- create a fpm project

```auto
> fpm new lorentz

```

- Add dependencies in fpm.toml

```auto
[dependencies]
rklib = { git="https://github.com/jacobwilliams/rklib.git" }
ogpf = { git = "https://github.com/kookma/ogpf.git" }

```

(it also requires having gnuplot installed)

---

<div class="post-metadata">

**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\
**Post date:** [January 11, 2024, 2:28pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/59 "2024-01-11T14:28:35Z")

</div>

I just came up with a small program to determine an array containing only the unique elements of another array (not well-tested, but at least for the test case it works):

```auto
! unique.f90 --
! Construct an array consisting of unique elements
!
module unique_elements
    implicit none

contains
recursive function unique( array ) result(unique_array)
    character(len=*), dimension(:) :: array
    character(len=len(array)), allocatable :: unique_array(:)

    if ( size(array) > 0 ) then
        unique_array = [array(1), unique( pack(array(2:), array(2:) /= array(1)) )]
    else
        allocate( unique_array(0) )
    endif
end function unique

end module unique_elements
program test_unique
    use unique_elements
    implicit none

    character(len = 2) :: array(8) = ['AA', 'BB', 'CB', 'AA', 'BB', 'CB', 'DD', 'AA']
    character(len = len(array)), allocatable :: unique_values(:)

    unique_values = unique(array)

    write(*,'(10a3)') unique_values
end program test_unique

```

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [January 11, 2024, 4:16pm UTC](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150/60 "2024-01-11T16:16:37Z")

</div>

Thanks. In general, a library for set operations in Fortran such as those in the Python [set](https://docs.python.org/3/tutorial/datastructures.html) would be useful. John Burkardt has [SET\_THEORY - An Implementation of Set Theoretic Operations](https://people.math.sc.edu/Burkardt/f_src/set_theory/set_theory.html), and for his codes,  
`gfortran -std=legacy set_theory.f90 set_theory_test.f90` compiles and runs.

[Previous page](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150.md?page=2)

[Next page](https://fortran-lang.discourse.group/t/fortran-code-snippets/2150.md?page=4)