🔢 Fortran Deep Foundations

🎯 Fortran Design Philosophy & Execution Model

Fortran (Formula Translation, John Backus, IBM, 1957) is the oldest high-level programming language, designed specifically for scientific and engineering computation. Modern Fortran (Fortran 90, 95, 2003, 2008, 2018) is a compiled, statically-typed language with native array operations, parallel programming support, and mathematical excellence.

Deep Explanation

Fortran uses 1-based indexing by default (mathematical convention), not 0-based like C. Variables are case-insensitive. The language has undergone major revisions: Fortran 77 (fixed form, columns 1-72), Fortran 90 (free form, modules, arrays), Fortran 2003 (OOP, C interop), Fortran 2018 (coarrays, teams).

Why it matters: Fortran powers weather forecasting (ECMWF), climate modeling (CESM), computational fluid dynamics (OpenFOAM), nuclear simulation, and most scientific supercomputing applications. It remains the fastest language for numerical linear algebra.

📌 1-Based Indexing📌 Column-Major Order📌 Static Typing📌 HPC Legacy

! Fortran Deep Foundations
! Comments start with ! (exclamation mark)

program hello_world
    implicit none
    ! 'implicit none' is mandatory in modern Fortran
    ! Prevents implicit typing (variables starting with I-N are integer by default in old Fortran)
    
    print *, "Hello, Fortran World!"
    print *, "Modern Fortran (2018) is alive and well!"
end program hello_world

! Compilation:
! gfortran hello.f90 -o hello
! ./hello

! Basic program structure:
! program name
!   use statements
!   implicit none
!   variable declarations
!   executable statements
! end program

📊 Data Types & Declarations

🎯 Intrinsic Data Types & KIND

Fortran has INTEGER, REAL, COMPLEX, LOGICAL, CHARACTER, and derived types. KIND parameters specify precision.

Precision Control & Portability

Use selected_int_kind and selected_real_kind for portable precision. Common KIND values: real64 = selected_real_kind(15, 307) for double precision. Fortran 2008+ supports iso_fortran_env intrinsic module with real64, real32, int32.

📐 INTEGER/REAL/COMPLEX⚡ KIND Parameters📊 SELECTED_REAL_KIND🔍 iso_fortran_env

program data_types
    use iso_fortran_env, only: int32, int64, real32, real64
    implicit none
    
    ! Basic types
    integer :: i = 42
    real :: x = 3.14159
    complex :: z = (1.0, 2.0)
    logical :: flag = .true.
    character(len=20) :: name = "Fortran"
    
    ! KIND-specified types (portable precision)
    integer(kind=int64) :: big_num = 123456789012345_int64
    real(kind=real64) :: pi = 3.14159265358979323846_real64
    real(kind=real32) :: single_prec = 3.14159_real32
    
    ! Kind inquiry
    print *, "Integer kind for i:", kind(i)
    print *, "Real kind for pi:", kind(pi)
    print *, "Range of int64:", huge(big_num)
    
    ! Selected kind (portable)
    integer, parameter :: dp = selected_real_kind(15, 307)
    real(kind=dp) :: precise = 1.0_dp
    print *, "Double precision kind:", dp
    
    ! Constants with underscores for readability
    integer, parameter :: million = 1_000_000
    real, parameter :: avogadro = 6.02214076e23
end program data_types

📐 Array Programming (DOPE Vectors)

🎯 Fortran's Array Heritage

Fortran arrays are first-class citizens with whole-array operations, array sections, and dynamic allocation. Arrays are column-major (first index varies fastest).

Array Declarations & Slicing

Arrays can be explicit-shape (size known at compile time), assumed-shape (dummy arguments), deferred-shape (allocatable), or assumed-size (legacy). The DOPE vector is the internal descriptor for assumed-shape arrays.

📐 Column-Major Order📊 Array Constructors🔀 Array Sections⚡ ALLOCATABLE

program array_demo
    implicit none
    
    ! Explicit-shape arrays (fixed size)
    integer, dimension(5) :: vec = [1, 2, 3, 4, 5]
    real, dimension(3,4) :: mat = 0.0
    
    ! Array constructors
    integer :: a(5) = (/1, 2, 3, 4, 5/)
    integer :: b(5) = [1, 2, 3, 4, 5]          ! Modern syntax
    
    ! Allocatable arrays (dynamic)
    integer, allocatable :: dynamic(:,:)
    integer :: n = 100
    allocate(dynamic(n, n))
    
    ! Initialize array with implied DO
    real :: values(10) = [(real(i), i=1, 10)]
    
    ! Array sections (slicing)
    integer :: arr(10) = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    integer :: first_five(5)
    first_five = arr(1:5)          ! arr(1:5:1) with stride
    integer :: even_indices(5)
    even_indices = arr(2:10:2)     ! stride of 2: [2,4,6,8,10]
    
    ! Multi-dimensional slicing
    integer :: matrix(5,5)
    integer :: row3(5)
    row3 = matrix(3, :)             ! entire row 3
    integer :: col2(5)
    col2 = matrix(:, 2)             ! entire column 2
    integer :: submat(2,2)
    submat = matrix(1:2, 4:5)       ! 2x2 submatrix
    
    ! Vector subscripting
    integer :: idx(3) = [1, 3, 5]
    integer :: selected(3)
    selected = arr(idx)              ! [1, 3, 5]
    
    deallocate(dynamic)
end program array_demo

⚡ Array Operations & Slicing

🎯 Whole-Array Arithmetic

Fortran's array operations operate on entire arrays without explicit loops — the original vectorization language.

Element-wise & Reduction

Arithmetic operations (+, -, *, /, **) on same-shape arrays are element-wise. Intrinsic functions include sum, product, matmul, dot_product, all, any, count, maxval, minval, maxloc, minloc, reshape, transpose, cshift, eoshift. The WHERE construct masks array operations.

🔢 Element-wise ops📊 SUM, PRODUCT🔀 MATMUL🎯 WHERE masks

program array_ops
    implicit none
    
    real :: a(5) = [1.0, 2.0, 3.0, 4.0, 5.0]
    real :: b(5) = [2.0, 4.0, 6.0, 8.0, 10.0]
    real :: c(5)
    
    ! Element-wise operations (no loops!)
    c = a + b           ! [3,6,9,12,15]
    c = a * b           ! [2,8,18,32,50]
    c = a ** 2          ! [1,4,9,16,25]
    c = 2.0 * a         ! [2,4,6,8,10]
    
    ! Reduction intrinsics
    real :: total = sum(a)              ! 15.0
    real :: product = product(a)        ! 120.0
    real :: avg = sum(a) / size(a)      ! average
    real :: max_val = maxval(a)         ! 5.0
    integer :: max_pos = maxloc(a, 1)   ! 5 (position)
    
    ! Matrix operations
    real :: m1(3,3), m2(3,3), m3(3,3)
    m3 = matmul(m1, m2)                 ! matrix multiplication
    real :: dot = dot_product(a(1:3), b(1:3))  ! dot product
    
    ! WHERE (conditional array operations)
    real :: x(10) = [(real(i), i=1,10)]
    where (x > 5)
        x = x * 2.0
    elsewhere
        x = x * 0.5
    end where
    
    ! FORALL (parallel assignment - Fortran 95)
    integer :: matrix(10,10)
    forall (i=1:10, j=1:10, i==j) matrix(i,j) = 1
    forall (i=1:10, j=1:10, i/=j) matrix(i,j) = 0
    
end program array_ops

🔀 Control Structures (IF, SELECT, DO)

🎯 Modern Control Flow

Fortran has IF-THEN-ELSE, SELECT CASE, DO loops (with and without indices), and DO WHILE.

Loop Control & EXIT/CYCLE

DO CONCURRENT (Fortran 2008) indicates loop is safe for parallel execution. EXIT exits a loop, CYCLE skips to next iteration. Named loops allow nested control. SELECT CASE works with integer, character, or logical expressions.

🔀 IF-THEN-ELSE📊 SELECT CASE🔄 DO loops⚡ DO CONCURRENT

program control_flow
    implicit none
    integer :: i, x, y
    
    ! IF-THEN-ELSE
    if (x > 0) then
        print *, "Positive"
    else if (x == 0) then
        print *, "Zero"
    else
        print *, "Negative"
    end if
    
    ! SELECT CASE
    select case (x)
    case (1)
        print *, "One"
    case (2, 3, 4)
        print *, "Two, three, or four"
    case (5:10)
        print *, "Five through ten"
    case default
        print *, "Other"
    end select
    
    ! Simple DO loop
    do i = 1, 10
        print *, i
    end do
    
    ! DO with step
    do i = 10, 1, -1
        print *, i
    end do
    
    ! DO WHILE
    i = 1
    do while (i <= 10)
        print *, i
        i = i + 1
    end do
    
    ! Infinite DO with EXIT
    i = 1
    do
        if (i > 10) exit
        if (mod(i, 2) == 0) then
            i = i + 1
            cycle           ! skip to next iteration
        end if
        print *, i
        i = i + 1
    end do
    
    ! Named loops (for nested control)
    outer: do i = 1, 10
        inner: do j = 1, 10
            if (i == j) cycle outer
            print *, i, j
        end do inner
    end do outer
    
    ! DO CONCURRENT (parallel hint)
    integer :: A(1000)
    do concurrent (i = 1:1000)
        A(i) = i * i
    end do
    
end program control_flow

📦 Subroutines & Functions

🎯 SUBROUTINE vs FUNCTION

Subroutines (call) have no return value, functions return a scalar or array. Both support intent(in), intent(out), intent(inout).

Argument Intent & Purity

INTENT(IN) = read-only, INTENT(OUT) = write-only (undefined on entry), INTENT(INOUT) = read/write. PURE functions have no side effects, ELEMENTAL functions work on arrays element-wise. OPTIONAL arguments can be omitted.

📦 SUBROUTINE📊 FUNCTION🔧 INTENT⚡ ELEMENTAL

program functions_demo
    implicit none
    
    ! Function call
    real :: x = 5.0, y
    y = square(x)
    print *, "Square of", x, "is", y
    
    ! Subroutine call
    call swap(x, y)
    print *, "After swap:", x, y
    
    ! Elemental function (works on arrays)
    real :: arr(5) = [1.0, 2.0, 3.0, 4.0, 5.0]
    arr = half(arr)          ! works element-wise
    print *, arr
    
contains

    ! FUNCTION (returns a value)
    function square(x) result(res)
        real, intent(in) :: x
        real :: res
        res = x * x
    end function square
    
    ! Alternative simple function syntax
    real function cube(x)
        real, intent(in) :: x
        cube = x * x * x
    end function cube
    
    ! SUBROUTINE (no return value)
    subroutine swap(a, b)
        real, intent(inout) :: a, b
        real :: temp
        temp = a
        a = b
        b = temp
    end subroutine swap
    
    ! ELEMENTAL function (works on array elements)
    elemental real function half(x)
        real, intent(in) :: x
        half = x / 2.0
    end function half
    
    ! PURE function (no side effects)
    pure real function add(a, b)
        real, intent(in) :: a, b
        add = a + b
    end function add
    
    ! OPTIONAL argument
    subroutine greet(name, greeting)
        character(*), intent(in) :: name
        character(*), optional, intent(in) :: greeting
        character(20) :: msg
        
        if (present(greeting)) then
            msg = greeting
        else
            msg = "Hello"
        end if
        print *, msg, name
    end subroutine greet
    
end program functions_demo

📁 Modules & Interfaces

🎯 Code Organization & Encapsulation

Modules are Fortran's namespace and encapsulation mechanism — containing variables, derived types, procedures, and interfaces.

Public/Private Control

PUBLIC and PRIVATE attributes control visibility. USE module, ONLY: symbol imports specific items. Interface blocks provide explicit interfaces — mandatory for passing procedures as arguments.

📁 MODULE structure🔧 PUBLIC/PRIVATE📊 USE with ONLY🔍 INTERFACE

! Module definition
module math_constants
    implicit none
    private                     ! everything private by default
    public :: pi, e, factorial, gcd
    
    real, parameter, public :: pi = 3.141592653589793
    real, parameter, public :: e = 2.718281828459045
    
    interface
        ! Interface for external procedure
        subroutine external_thing(x)
            real, intent(inout) :: x
        end subroutine external_thing
    end interface
    
contains

    recursive function factorial(n) result(res)
        integer, intent(in) :: n
        integer :: res
        if (n <= 1) then
            res = 1
        else
            res = n * factorial(n-1)
        end if
    end function factorial
    
    function gcd(a, b) result(res)
        integer, intent(in) :: a, b
        integer :: res
        integer :: ta, tb
        ta = a
        tb = b
        do while (tb /= 0)
            res = mod(ta, tb)
            ta = tb
            tb = res
        end do
        res = ta
    end function gcd

end module math_constants

! Submodule (Fortran 2008+ for splitting modules)
submodule (math_constants) math_impl
contains
    
    module procedure factorial
        ! Implementation goes here (optional)
    end procedure
    
end submodule math_impl

! Main program using the module
program use_math
    use math_constants, only: pi, e, factorial, gcd
    implicit none
    
    print *, "pi =", pi
    print *, "e =", e
    print *, "10! =", factorial(10)
    print *, "GCD(48, 18) =", gcd(48, 18)
    
end program use_math

🏗️ Derived Types & OOP

🎯 TYPE & CLASS (Fortran 2003+)

Derived types are Fortran's struct-like constructs. Fortran 2003 added OOP: inheritance, polymorphism, type-bound procedures, and abstract types.

Type-Bound Procedures & Polymorphism

TYPE :: ... END TYPE defines components. CONTAINS inside type binds procedures. CLASS declares polymorphic variables. EXTENDS enables inheritance. ABSTRACT defines abstract types.

📦 TYPE definition🔧 CONTAINS (methods)📊 EXTENDS (inheritance)⚡ CLASS (polymorphic)

program oop_demo
    implicit none
    
    ! Simple derived type
    type :: point
        real :: x, y
    end type point
    
    ! Derived type with type-bound procedures
    type :: vector
        real :: x, y, z
    contains
        procedure :: length => vector_length
        procedure :: normalize => vector_normalize
        procedure :: dot => vector_dot
        generic :: operator(+) => vector_add
        generic :: operator(*) => scalar_mult
    end type vector
    
    ! Inheritance
    type, extends(vector) :: colored_vector
        character(20) :: color
    contains
        procedure :: describe => colored_describe
    end type colored_vector
    
    ! Abstract type
    type, abstract :: shape
    contains
        procedure(area_interface), deferred :: area
    end type shape
    
    abstract interface
        real function area_interface(this)
            import :: shape
            class(shape), intent(in) :: this
        end function area_interface
    end interface
    
    ! Type extending abstract
    type, extends(shape) :: circle
        real :: radius
    contains
        procedure :: area => circle_area
    end type circle
    
    type(point) :: p = point(3.0, 4.0)
    type(vector) :: v1, v2, v3
    type(colored_vector) :: cv
    
    v1 = vector(1.0, 2.0, 3.0)
    v2 = vector(4.0, 5.0, 6.0)
    v3 = v1 + v2      ! uses defined operator
    
    print *, "Vector length:", v1%length()
    print *, "Dot product:", v1%dot(v2)
    
    cv = colored_vector(1.0, 0.0, 0.0, "red")
    call cv%describe()
    
    ! Polymorphism
    class(shape), allocatable :: s
    allocate(circle :: s)
    select type(s)
    type is (circle)
        s%radius = 5.0
        print *, "Circle area:", s%area()
    end select
    
contains

    real function vector_length(this)
        class(vector), intent(in) :: this
        vector_length = sqrt(this%x**2 + this%y**2 + this%z**2)
    end function
    
    subroutine vector_normalize(this)
        class(vector), intent(inout) :: this
        real :: len
        len = this%length()
        if (len > 0) then
            this%x = this%x / len
            this%y = this%y / len
            this%z = this%z / len
        end if
    end subroutine
    
    real function vector_dot(this, other)
        class(vector), intent(in) :: this
        type(vector), intent(in) :: other
        vector_dot = this%x*other%x + this%y*other%y + this%z*other%z
    end function
    
    type(vector) function vector_add(u, v)
        type(vector), intent(in) :: u, v
        vector_add%x = u%x + v%x
        vector_add%y = u%y + v%y
        vector_add%z = u%z + v%z
    end function
    
    type(vector) function scalar_mult(a, v)
        real, intent(in) :: a
        type(vector), intent(in) :: v
        scalar_mult%x = a * v%x
        scalar_mult%y = a * v%y
        scalar_mult%z = a * v%z
    end function
    
    subroutine colored_describe(this)
        class(colored_vector), intent(in) :: this
        print *, "Vector with color", this%color
        print *, "Components:", this%x, this%y, this%z
    end subroutine
    
    real function circle_area(this)
        class(circle), intent(in) :: this
        circle_area = 3.141592653589793 * this%radius**2
    end function

end program oop_demo

💾 I/O & Formatting (READ, WRITE)

🎯 READ, WRITE, PRINT, OPEN, CLOSE

Fortran has powerful formatted I/O with format descriptors, internal I/O (reading/writing from strings), and non-advancing I/O.

Format Specifiers & File Operations

Format descriptors: I (integer), F (real), E (scientific), A (character), L (logical), X (space), / (newline). Internal file I/O uses character variables. NAMED and STREAM access available.

📁 OPEN/CLOSE📊 FORMAT specifiers🔧 Internal I/O⚡ NAMELIST

program io_demo
    implicit none
    integer :: i, ios
    real :: x, y
    character(100) :: line
    character(20) :: name
    integer :: unit = 10
    
    ! Simple output
    print *, "Hello World"
    write(*, *) "Hello World"      ! same as print
    
    ! Formatted output
    print '(A, I5, F10.3)', "Value:", 42, 3.14159
    write(*, '(A, I0, A, F0.3)') "Result = ", 42, " pi = ", 3.14159
    
    ! Format specifiers:
    ! I5   - integer with width 5
    ! F10.3 - real with total width 10, 3 decimal places
    ! E12.5 - scientific notation
    ! A10   - character with width 10
    ! L1    - logical (T/F)
    ! X     - space
    ! /     - newline
    
    ! Reading input
    read(*, *) i, x, name
    
    ! Formatted read
    read(*, '(I5, F10.3, A20)') i, x, name
    
    ! File operations
    open(unit=unit, file="output.txt", status="replace", action="write", iostat=ios)
    if (ios == 0) then
        write(unit, *) "Writing to file"
        write(unit, '(A, I0)') "Count = ", 100
        close(unit)
    else
        print *, "Error opening file"
    end if
    
    ! Reading from file
    open(10, file="data.txt", status="old", action="read")
    do
        read(10, *, iostat=ios) i, x
        if (ios /= 0) exit
        print *, i, x
    end do
    close(10)
    
    ! Internal I/O (read from / write to string)
    line = "42 3.14"
    read(line, *) i, x
    write(line, '(I5, F10.3)') i, x
    print *, trim(line)
    
    ! NAMELIST (grouped I/O)
    namelist /params/ i, x, name
    i = 10
    x = 2.718
    name = "Fortran"
    open(20, file="params.nml", status="replace")
    write(20, nml=params)
    close(20)
    
    ! Stream I/O (unformatted, direct access)
    open(30, file="binary.dat", form="unformatted", access="stream")
    write(30) 42, 3.14159, "Hello"
    close(30)
    
end program io_demo

📍 Pointers & Allocatables

🎯 ALLOCATABLE vs POINTER

Allocatables (Fortran 90) are automatically deallocated when out of scope. Pointers require manual deallocation and can point to targets.

Dynamic Memory Management

ALLOCATABLE arrays are simpler and safer — compiler handles cleanup. POINTER attributes allow aliasing, linked structures, and pointer arithmetic (limited). TARGET attribute marks variables as pointer-targetable. CONTIGUOUS attribute ensures array contiguity.

📍 ALLOCATE/DEALLOCATE🔗 POINTER/TARGET📊 ALLOCATABLE⚡ NULLIFY

program pointers_demo
    implicit none
    
    ! Allocatable arrays (automatic cleanup)
    integer, allocatable :: arr(:), matrix(:,:)
    integer :: n = 100
    
    allocate(arr(n))
    arr = [(i, i=1, n)]
    print *, "Allocated arr of size", size(arr)
    deallocate(arr)          ! explicit deallocation
    ! arr automatically deallocated when out of scope if not done
    
    ! Reallocate on assignment (Fortran 2003+)
    arr = [1, 2, 3, 4, 5]   ! automatically allocates if needed
    
    ! Pointers
    integer, pointer :: p(:), q(:)
    integer, target :: t(5) = [1, 2, 3, 4, 5]
    integer, pointer :: r => null()   ! initialize to null
    
    ! Pointer assignment
    p => t
    print *, p(3)            ! 3
    
    ! Pointer to array section
    q => t(2:4)
    print *, q               ! [2, 3, 4]
    
    ! Allocatable pointer
    allocate(p(10))
    p = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    
    ! Pointer association checking
    if (associated(p)) then
        print *, "p is associated"
    end if
    
    ! Linked list using derived types with pointers
    type :: node
        integer :: value
        type(node), pointer :: next => null()
    end type node
    
    type(node), pointer :: head, current
    integer :: i
    
    ! Create linked list
    nullify(head)
    do i = 1, 5
        allocate(current)
        current%value = i
        current%next => head
        head => current
    end do
    
    ! Traverse linked list
    current => head
    do while (associated(current))
        print *, current%value
        current => current%next
    end do
    
    ! Clean up linked list
    do while (associated(head))
        current => head
        head => head%next
        deallocate(current)
    end do
    
    ! Clean up pointer
    if (associated(p)) deallocate(p)
    nullify(p)
    
end program pointers_demo

🔧 Generic & Pure Procedures

🎯 INTERFACE for Generics

Generic interfaces allow multiple procedures with the same name but different argument types. PURE procedures guarantee no side effects.

Overloading & Side-Effect Free

INTERFACE blocks create generic names. PURE procedures are guaranteed side-effect free (enables compiler optimizations). ELEMENTAL procedures are PURE by default and work on arrays element-wise. IMPURE ELEMENTAL allows side effects.

📐 Generic INTERFACE⚡ PURE attribute🔁 ELEMENTAL📊 IMPURE ELEMENTAL

program generic_demo
    implicit none
    
    interface swap
        module procedure swap_int, swap_real, swap_char
    end interface swap
    
    interface to_string
        module procedure int_to_string, real_to_string
    end interface to_string
    
    integer :: a = 5, b = 10
    real :: x = 3.14, y = 2.718
    character(10) :: s1 = "Hello", s2 = "World"
    
    call swap(a, b)      ! integer swap
    call swap(x, y)      ! real swap
    call swap(s1, s2)    ! character swap
    
    print *, "A=", a, "B=", b
    print *, to_string(42)
    print *, to_string(3.14159)
    
contains

    ! Pure function (no side effects)
    pure integer function square(n)
        integer, intent(in) :: n
        square = n * n
    end function square
    
    ! Elemental function (works on array elements automatically)
    elemental real function cube(x)
        real, intent(in) :: x
        cube = x * x * x
    end function cube
    
    ! Generic implementations
    subroutine swap_int(a, b)
        integer, intent(inout) :: a, b
        integer :: temp
        temp = a; a = b; b = temp
    end subroutine
    
    subroutine swap_real(a, b)
        real, intent(inout) :: a, b
        real :: temp
        temp = a; a = b; b = temp
    end subroutine
    
    subroutine swap_char(a, b)
        character(*), intent(inout) :: a, b
        character(len(a)) :: temp
        temp = a; a = b; b = temp
    end subroutine
    
    function int_to_string(n) result(str)
        integer, intent(in) :: n
        character(20) :: str
        write(str, '(I0)') n
    end function
    
    function real_to_string(x) result(str)
        real, intent(in) :: x
        character(20) :: str
        write(str, '(F0.6)') x
    end function
    
end program generic_demo

🌍 Coarrays (Parallel Programming)

🎯 Fortran 2008/2018 Parallelism

Coarrays provide PGAS (Partitioned Global Address Space) parallel programming — each image has its own data with remote access via square bracket syntax.

Image-Based Parallelism

Program launches multiple images (processes). this_image() returns current image number, num_images() returns total. Access remote data: arr[j]@image. SYNC ALL synchronizes images. LOCK/UNLOCK for mutual exclusion.

🌍 this_image()📊 num_images()🔀 SYNC ALL🔒 LOCK/UNLOCK

! Compile with: gfortran -fcoarray=single hello_coarray.f90
! Or: gfortran -fcoarray=lib hello_coarray.f90 -lcaf_mpi

program coarray_demo
    implicit none
    
    ! Coarray variable (accessible from all images)
    real, codimension[*] :: shared_value
    integer :: me, nimages, i
    
    me = this_image()
    nimages = num_images()
    
    ! Each image initializes its own copy
    shared_value = real(me)
    
    ! Synchronize all images
    sync all
    
    ! Image 1 collects values from all images
    if (me == 1) then
        print *, "Values from all images:"
        do i = 1, nimages
            print *, "Image", i, ":", shared_value[i]
        end do
    end if
    
    ! SYNC ALL ensures all images have consistent state
    sync all
    
    ! Collective reduction (sum across images)
    real :: sum_val
    if (me == 1) then
        sum_val = shared_value[1]
        do i = 2, nimages
            sum_val = sum_val + shared_value[i]
        end do
        print *, "Sum across images:", sum_val
    end if
    
end program coarray_demo

! More complex coarray example: parallel dot product
program coarray_dot
    use iso_fortran_env, only: real64
    implicit none
    
    integer, parameter :: n = 1000000
    real(real64), allocatable :: a(:), b(:)
    real(real64), codimension[*] :: partial_sum
    real(real64) :: total_sum
    integer :: me, nimages, i, chunk, start_idx, end_idx
    
    me = this_image()
    nimages = num_images()
    
    ! Distribute work across images
    chunk = n / nimages
    start_idx = (me-1)*chunk + 1
    if (me == nimages) then
        end_idx = n
    else
        end_idx = me*chunk
    end if
    
    ! Each image works on its chunk
    allocate(a(start_idx:end_idx), b(start_idx:end_idx))
    a = [(real(i, real64), i=start_idx, end_idx)]
    b = 1.0_real64
    
    ! Compute partial dot product
    partial_sum = sum(a * b)
    
    ! Synchronize before reduction
    sync all
    
    ! Image 1 collects and sums
    if (me == 1) then
        total_sum = partial_sum[1]
        do i = 2, nimages
            total_sum = total_sum + partial_sum[i]
        end do
        print *, "Total dot product:", total_sum
    end if
    
    deallocate(a, b)
    
end program coarray_dot

! Lock example (mutual exclusion)
program coarray_lock
    implicit none
    integer, codimension[*] :: counter
    integer(8), codimension[*] :: lock_var
    integer :: me
    
    me = this_image()
    counter = 0
    
    sync all
    
    ! Critical section using lock
    lock(lock_var)
    counter = counter + 1
    unlock(lock_var)
    
    sync all
    
    if (me == 1) then
        print *, "Final counter:", counter[1]
    end if
    
end program coarray_lock

🔄 C Interoperability (ISO_C_BINDING)

🎯 ISO_C_BINDING Module

Fortran 2003 introduced C interoperability, allowing direct calling between Fortran and C with standardized type mappings.

BIND(C) & C Types

USE, INTRINSIC :: ISO_C_BINDING provides C-compatible types: C_INT, C_DOUBLE, C_CHAR, C_PTR, etc. BIND(C) attribute ensures C calling convention. C_LOC gets C pointer to data.

🔗 ISO_C_BINDING📊 BIND(C)🔧 C_LOC, C_F_POINTER⚡ VALUE attribute

! Fortran code calling C
module fortran_c_interop
    use, intrinsic :: iso_c_binding
    implicit none
    
    interface
        ! Call C function: double c_sqrt(double x)
        real(c_double) function c_sqrt(x) bind(C, name="sqrt")
            import :: c_double
            real(c_double), value :: x
        end function c_sqrt
        
        ! Call C function: void c_print(const char* str)
        subroutine c_print(str) bind(C, name="print_string")
            import :: c_char
            character(kind=c_char), dimension(*), intent(in) :: str
        end subroutine c_print
        
        ! Call C function with array
        subroutine c_sum_array(arr, n, result) bind(C, name="sum_array")
            import :: c_int, c_double
            real(c_double), dimension(*), intent(in) :: arr
            integer(c_int), value :: n
            real(c_double), intent(out) :: result
        end subroutine c_sum_array
    end interface
    
contains

    ! Fortran function that can be called from C
    real(c_double) function fortran_add(a, b) bind(C, name="fortran_add")
        real(c_double), value :: a, b
        fortran_add = a + b
    end function fortran_add
    
    ! Working with C pointers
    subroutine use_c_pointer(c_ptr_array, n)
        type(c_ptr), value :: c_ptr_array
        integer(c_int), value :: n
        real(c_double), pointer :: fortran_array(:)
        
        ! Convert C pointer to Fortran pointer
        call c_f_pointer(c_ptr_array, fortran_array, [n])
        
        ! Use as normal Fortran array
        print *, "Array values:", fortran_array
    end subroutine use_c_pointer

end module fortran_c_interop

! Main program using C interop
program c_interop_demo
    use, intrinsic :: iso_c_binding
    use fortran_c_interop
    implicit none
    
    real(c_double) :: x = 2.0_c_double
    real(c_double) :: result
    character(kind=c_char), dimension(100) :: c_string
    
    ! Call C sqrt function
    result = c_sqrt(x)
    print *, "sqrt(", x, ") =", result
    
    ! Prepare C string (null-terminated)
    c_string = transfer("Hello from Fortran!" // c_null_char, c_string)
    call c_print(c_string)
    
    ! Call Fortran from C (through header)
    print *, "2.0 + 3.0 =", fortran_add(2.0_c_double, 3.0_c_double)
    
end program c_interop_demo

! Example C function implementations (in separate .c file):
/*
#include 
#include 

double c_sqrt(double x) {
    return sqrt(x);
}

void print_string(const char* str) {
    printf("%s\n", str);
}

void sum_array(double* arr, int n, double* result) {
    *result = 0.0;
    for (int i = 0; i < n; i++) {
        *result += arr[i];
    }
}
*/

🚀 HPC & Optimization (SIMD, OPENMP)

🎯 High-Performance Computing

Fortran is the king of HPC — with native optimizations, SIMD directives, and seamless OpenMP integration.

Compiler Optimizations & Directives

!$OMP directives for shared-memory parallelism. !DIR$ IVDEP and !DIR$ SIMD for vectorization hints. CONTIGUOUS attribute for pointer/assumed-shape arrays. DO CONCURRENT for safe parallelism. assume_aligned for alignment hints.

🔧 OpenMP directives⚡ SIMD vectorization📊 DO CONCURRENT🔀 CONTIGUOUS

program hpc_optimization
    implicit none
    integer, parameter :: n = 10000000
    real, allocatable :: a(:), b(:), c(:)
    integer :: i
    real :: start_time, end_time
    
    allocate(a(n), b(n), c(n))
    a = 1.0
    b = 2.0
    
    call cpu_time(start_time)
    
    ! Basic loop (autovectorized by compiler with -O3)
    do i = 1, n
        c(i) = a(i) + b(i)
    end do
    
    call cpu_time(end_time)
    print *, "Serial time:", end_time - start_time
    
    ! DO CONCURRENT (parallel hint)
    call cpu_time(start_time)
    do concurrent (i = 1:n)
        c(i) = a(i) * b(i)
    end do
    call cpu_time(end_time)
    print *, "DO CONCURRENT time:", end_time - start_time
    
    ! Whole-array operation (fastest)
    call cpu_time(start_time)
    c = a - b
    call cpu_time(end_time)
    print *, "Array op time:", end_time - start_time
    
    deallocate(a, b, c)
    
end program hpc_optimization

! OpenMP parallel example
! Compile with: gfortran -fopenmp omp_demo.f90

program omp_demo
    use omp_lib
    implicit none
    integer, parameter :: n = 100000000
    real, allocatable :: a(:), b(:), c(:)
    integer :: i, thread_num
    
    allocate(a(n), b(n), c(n))
    a = 1.0
    b = 2.0
    
    !$omp parallel do private(i) shared(a,b,c)
    do i = 1, n
        c(i) = a(i) + b(i)
    end do
    !$omp end parallel do
    
    !$omp parallel private(thread_num)
    thread_num = omp_get_thread_num()
    !$omp single
    print *, "Number of threads:", omp_get_num_threads()
    !$omp end single
    !$omp end parallel
    
    ! Reduction example
    real :: sum_val = 0.0
    
    !$omp parallel do reduction(+:sum_val)
    do i = 1, n
        sum_val = sum_val + a(i)
    end do
    !$omp end parallel do
    
    print *, "Sum:", sum_val
    
    deallocate(a, b, c)
    
end program omp_demo

! SIMD directives
program simd_demo
    implicit none
    integer, parameter :: n = 1000000
    real :: a(n), b(n), c(n)
    integer :: i
    
    !DIR$ SIMD
    do i = 1, n
        c(i) = a(i) * b(i)
    end do
    
    !DIR$ IVDEP (ignore vector dependencies)
    !DIR$ SIMD
    do i = 2, n
        a(i) = a(i-1) + b(i)
    end do
    
end program simd_demo

! Matrix multiplication optimization
program matmul_optimized
    implicit none
    integer, parameter :: n = 1024
    real :: a(n,n), b(n,n), c(n,n)
    integer :: i, j, k
    
    ! Initialize
    a = 1.0
    b = 2.0
    c = 0.0
    
    ! Optimized loop order (j-k-i for cache efficiency)
    ! Note: Fortran is column-major, so iterating i first is better
    do j = 1, n
        do k = 1, n
            do i = 1, n
                c(i,j) = c(i,j) + a(i,k) * b(k,j)
            end do
        end do
    end do
    
    ! Better: Use intrinsic matmul (highly optimized BLAS)
    c = matmul(a, b)
    
    ! Or call BLAS directly
    ! call sgemm('N', 'N', n, n, n, 1.0, a, n, b, n, 0.0, c, n)
    
end program matmul_optimized

! Compiler optimization flags for gfortran:
! -O3 -march=native -ffast-math -funroll-loops
! -fopenmp -fcoarray=single (or lib)
! -flto (link-time optimization)
! -ftree-vectorize -ftree-loop-vectorize