Fortran stdlib

program kinds_version_optval
  use stdlib_kinds, only: sp, dp, qp, xdp, int32, int64
  use stdlib_version, only: get_stdlib_version, stdlib_version_string, stdlib_version_compact
  use stdlib_optval, only: optval
  use stdlib_constants, only: pi => PI_dp, stefan_boltzmann
  implicit none
  integer :: major, minor, patch
  character(len=:), allocatable :: version

  call get_stdlib_version(major, minor, patch, version)
  print '(a, i0, a, i0, a, i0, a, a)', 'stdlib ', major, '.', minor, '.', patch, ' = ', version
  print '(a, a, a, i0)', 'compile-time string ', stdlib_version_string, ', compact ', stdlib_version_compact

  print '(a, 4(i0, 1x))', 'sp dp qp xdp = ', sp, dp, qp, xdp
  print '(a, 2(i0, 1x))', 'int32 int64  = ', int32, int64

  print '(a, f8.6)', 'diffusion length, default scale  ', diffusion_length(0.01_dp, 4.0_dp)
  print '(a, f8.6)', 'diffusion length, scale = 2      ', diffusion_length(0.01_dp, 4.0_dp, 2.0_dp)
  print '(a, es12.5)', 'Stefan-Boltzmann constant        ', stefan_boltzmann
  print '(a, f0.10)', 'PI_dp                            ', pi

contains

  pure function diffusion_length(alpha, t, scale) result(l)
    real(dp), intent(in) :: alpha, t
    real(dp), intent(in), optional :: scale
    real(dp) :: l
    l = optval(scale, 1.0_dp) * sqrt(4.0_dp * alpha * t)
  end function diffusion_length

end program kinds_version_optval

stdlib 0.8.1 = 0.8.1
compile-time string 0.8.1, compact 801
sp dp qp xdp = 4 8 -1 -1
int32 int64  = 4 8
diffusion length, default scale  0.400000
diffusion length, scale = 2      0.800000
Stefan-Boltzmann constant         5.67037E-08
PI_dp                            3.1415926536

program stdlib_io_demo
  use stdlib_kinds, only: dp
  use stdlib_io, only: loadtxt, savetxt, open, get_line, get_file, FMT_REAL_DP
  use stdlib_io_npy, only: load_npy, save_npy
  implicit none
  integer, parameter :: nx = 5, ny = 3
  real(dp) :: u(nx, ny)
  real(dp), allocatable :: v(:, :), w(:, :)
  character(len=:), allocatable :: line, text
  integer :: i, j, unit, stat

  ! a small temperature field: u(i, j) = 10*(i-1) + (j-1)
  do j = 1, ny
    do i = 1, nx
      u(i, j) = 10.0_dp * real(i - 1, dp) + real(j - 1, dp)
    end do
  end do

  ! text round trip: one array row per line, columns separated by ';'
  call savetxt('field.csv', u, delimiter=';')
  call loadtxt('field.csv', v, delimiter=';')
  print '(a, 2(1x, i0), a, es9.2)', 'loaded shape', shape(v), ', max |u - v| = ', maxval(abs(u - v))

  ! skiprows and max_rows select a block of rows
  call loadtxt('field.csv', w, skiprows=1, max_rows=2, delimiter=';')
  print '(a, 2(1x, i0))', 'rows 2..3 shape', shape(w)

  ! open() with a mode string returns a unit; get_line reads a whole line
  unit = open('field.csv', 'r', iostat=stat)
  print '(a, i0)', 'open iostat = ', stat
  call get_line(unit, line, iostat=stat)
  print '(a, a)', 'first line: ', line
  close(unit)

  ! get_file reads the entire file into one allocatable string
  call get_file('field.csv', text)
  print '(a, i0, a)', 'file holds ', len(text), ' characters'

  ! .npy: binary and self-describing, any rank, readable by NumPy
  call save_npy('field.npy', u)
  call load_npy('field.npy', v)
  print '(a, es9.2)', 'npy round trip, max |u - v| = ', maxval(abs(u - v))

  ! format constant that prints a real(dp) with full precision
  print FMT_REAL_DP, 1.0_dp / 3.0_dp

  ! missing file with iostat, instead of the default error stop
  unit = open('does-not-exist.csv', 'r', iostat=stat)
  print '(a, l1)', 'missing file reported through iostat: ', stat /= 0
end program stdlib_io_demo

loaded shape 5 3, max |u - v| =  0.00E+00
rows 2..3 shape 2 3
open iostat = 0
first line:  0.0000000000000000E+000; 1.0000000000000000E+000; 2.0000000000000000E+000
file holds 375 characters
npy round trip, max |u - v| =  0.00E+00
 3.3333333333333331E-001
missing file reported through iostat: T

 0.0000000000000000E+000; 1.0000000000000000E+000; 2.0000000000000000E+000

program stdlib_linalg_demo
  use stdlib_kinds, only: dp
  use stdlib_linalg, only: solve, det, inv, eye, diag, trace, norm, eigvalsh, svdvals, cholesky, lstsq
  use stdlib_linalg_state, only: linalg_state_type
  implicit none
  integer, parameter :: n = 5
  real(dp), parameter :: r = 0.4_dp          ! alpha * dt / dx**2
  real(dp) :: a(n, n), u(n), l(n, n), c(n, n)
  real(dp) :: lambda(n), sigma(n), xs(6), ys(6), design(6, 2)
  real(dp), allocatable :: x(:), fit(:)
  type(linalg_state_type) :: err
  integer :: i

  ! One implicit (backward Euler) step of 1D heat flow: (I + r*K) u_new = u_old,
  ! where K = tridiag(-1, 2, -1) is the discrete second difference.
  a = (2.0_dp * r + 1.0_dp) * eye(n)
  do i = 1, n - 1
    a(i, i + 1) = -r
    a(i + 1, i) = -r
  end do
  u = [0.0_dp, 0.0_dp, 100.0_dp, 0.0_dp, 0.0_dp]

  x = solve(a, u)                               ! LAPACK *GESV underneath
  print '(a, 5f9.4)', 'u_new         =', x
  print '(a, es10.2)', '|A u_new - u|_2 = ', norm(matmul(a, x) - u, 2)
  print '(a, f9.4)', 'heat kept     =', sum(x)

  print '(a, f12.6)', 'det(A)        = ', det(a)
  print '(a, f12.6)', 'trace(A)      = ', trace(a)
  print '(a, es10.2)', '|A * inv(A) - I|_max = ', maxval(abs(matmul(a, inv(a)) - eye(n)))

  lambda = eigvalsh(a)                          ! symmetric: ascending eigenvalues
  print '(a, 5f8.4)', 'eigenvalues   =', lambda
  sigma = svdvals(a)                            ! descending singular values
  print '(a, 5f8.4)', 'cond_2(A)     =', sigma(1) / sigma(n)

  call cholesky(a, c, lower=.true.)             ! A = L L^T
  l = c
  print '(a, es10.2)', '|L L^T - A|_max = ', maxval(abs(matmul(l, transpose(l)) - a))
  print '(a, 5f8.4)', 'diag(L)       =', diag(l)

  ! least squares through the same module: fit T(x) = p + q x to noisy samples
  xs = [0.0_dp, 1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp]
  ys = [1.1_dp, 2.9_dp, 5.2_dp, 7.1_dp, 8.8_dp, 11.2_dp]
  design(:, 1) = 1.0_dp
  design(:, 2) = xs
  fit = lstsq(design, ys)
  print '(a, 2f8.4)', 'fit (p, q)    =', fit

  ! expert interface: a singular matrix reports an error instead of error stop
  a = 0.0_dp
  x = solve(a, u, err=err)
  print '(a, l1)', 'singular matrix flagged: ', err%error()
  print '(a)', err%print()
end program stdlib_linalg_demo

u_new         =   3.2206  14.4928  61.9968  14.4928   3.2206
|A u_new - u|_2 =   3.09E-15
heat kept     =  97.4235
det(A)        =    15.301440
trace(A)      =     9.000000
|A * inv(A) - I|_max =   2.22E-16
eigenvalues   =  1.1072  1.4000  1.8000  2.2000  2.4928
cond_2(A)     =  2.2515
|L L^T - A|_max =   2.22E-16
diag(L)       =  1.3416  1.3081  1.3063  1.3062  1.3062
fit (p, q)    =  1.0429  2.0029
singular matrix flagged: T
[solve] returned Algebra Error: singular matrix

program stdlib_logger_demo
  use stdlib_logger, only: logger_type, global_logger, debug_level, warning_level
  implicit none
  type(logger_type) :: run_log
  integer :: unit, stat, step

  ! The global logger writes to output_unit while its log_units list is empty.
  call global_logger%configure(time_stamp=.false.)
  call global_logger%log_information('heatflow started', module='demo', procedure='main')

  ! A local logger with its own file; messages below the level are ignored.
  call run_log%configure(time_stamp=.false., level=warning_level)
  call run_log%add_log_file('run.log', unit, stat=stat)
  print '(a, i0, a, i0)', 'log files opened: ', run_log%log_units_assigned(), ', stat = ', stat

  do step = 1, 3
    call run_log%log_debug('step detail')                 ! dropped: below warning_level
    call run_log%log_information('step finished')         ! dropped: below warning_level
  end do
  call run_log%log_warning('residual still large', module='solver', procedure='run')
  call run_log%log_error('did not converge in 3 steps', module='solver', procedure='run', stat=1, errmsg='max_steps reached')

  ! lower the threshold: everything is logged from now on
  call run_log%configure(level=debug_level)
  call run_log%log_debug('now visible')
  call run_log%log_message('a bare message without a label')

  call run_log%remove_log_unit(unit, close_unit=.true., stat=stat)
  print '(a, i0)', 'units left: ', run_log%log_units_assigned()

  ! echo the file content to show what was recorded
  call echo('run.log')

contains

  subroutine echo(file)
    character(len=*), intent(in) :: file
    character(len=200) :: line
    integer :: u, ios
    open(newunit=u, file=file, status='old', action='read')
    do
      read(u, '(a)', iostat=ios) line
      if (ios /= 0) exit
      print '(2a)', '| ', trim(line)
    end do
    close(u)
  end subroutine echo

end program stdlib_logger_demo

demo % main: INFO: heatflow started
log files opened: 1, stat = 0
units left: 0
| solver % run: WARN: residual still large
| solver % run: ERROR: did not converge in 3 steps
| With stat = 1
| With errmsg = "max_steps reached"
| DEBUG: now visible
| a bare message without a label

program stdlib_strings_demo
  use stdlib_kinds, only: dp
  use stdlib_strings, only: strip, chomp, starts_with, ends_with, slice, find, replace_all, &
                            count, join, to_string
  use stdlib_string_type, only: string_type, assignment(=), operator(//), operator(==), &
                                len, to_upper, write(formatted)
  implicit none
  character(len=*), parameter :: line = '   output = "heat_0042.csv"   '
  character(len=:), allocatable :: name
  type(string_type) :: s, t
  integer :: step

  ! stdlib_strings: functions over default character and string_type
  name = strip(line)
  print '(3a)', '[', name, ']'
  print '(3a)', '[', chomp('field.csv', '.csv'), ']'
  print '(l1, 1x, l1)', starts_with(name, 'output'), ends_with(name, '.csv"')
  print '(a)', slice('abcdefghij', 2, 6, 2)
  print '(i0)', find('heat_0042_0042.csv', '0042', 2)
  print '(i0)', count('heat_0042_0042.csv', '_')
  print '(a)', replace_all('heat_0042.csv', 'heat', 'flux')
  print '(a)', join(['nx', 'ny', 'lx'], ', ')

  step = 42
  print '(a)', 'file heat_' // to_string(step, '(i0.6)') // '.csv'
  print '(a)', to_string(3.14159_dp, '(f8.3)')
  print '(a)', to_string(.true.)

  ! string_type: a deferred-length string with the intrinsic string procedures overloaded
  s = 'Heat'
  t = s // ' flow'
  print '(dt)', t
  print '(i0)', len(t)
  print '(l1)', t == 'Heat flow'
  t = to_upper(t)
  write (*, '(dt)') t
end program stdlib_strings_demo

[output = "heat_0042.csv"]
[field]
T T
bdf
11
2
flux_0042.csv
nx, ny, lx
file heat_000042.csv
   3.142
T
Heat flow
9
T
HEAT FLOW

program stdlib_padding_demo
  use stdlib_strings, only: padl, padr, zfill
  implicit none
  print '(3a)', '[', padl('7', 4, '0'), ']'
  print '(3a)', '[', padr('ny', 6), ']'
  print '(3a)', '[', zfill('42', 6), ']'
end program stdlib_padding_demo

05b_padding.f90:2:6:

    2 |   use stdlib_strings, only: padl, padr, zfill
      |      1
Warning: Type specified for intrinsic function 'len' at (1) is ignored [-Wsurprising]

[0007]
[ny    ]
[000042]

program stdlib_text_demo
  use stdlib_kinds, only: sp, dp, int32
  use stdlib_ascii, only: is_alpha, is_digit, is_white, is_punctuation, to_upper, to_lower, to_title, &
                          reverse, TAB, digits
  use stdlib_str2num, only: to_num
  use stdlib_stringlist_type, only: stringlist_type, fidx, bidx
  use stdlib_string_type, only: string_type, write(formatted)
  implicit none
  character(len=*), parameter :: record = 'nx=64'//TAB//'ny=32'//TAB//'alpha=1.0e-2'
  type(stringlist_type) :: keys
  type(string_type) :: first
  integer :: i, nd, start
  integer :: n(2)
  real(dp) :: alpha

  ! character classification, one character at a time
  nd = 0
  do i = 1, len(record)
    if (is_digit(record(i:i))) nd = nd + 1
  end do
  print '(a, i0, a, l1, l1, l1)', 'digits in record: ', nd, '; classes of "a", TAB, ",": ', &
    is_alpha('a'), is_white(TAB), is_punctuation(',')
  print '(2a)', 'digit set: ', digits

  ! case conversion on default characters
  print '(a)', to_upper('heatflow')
  print '(a)', to_lower('HEATFLOW')
  print '(a)', to_title('heat flow solver')
  print '(a)', reverse('abc')

  ! to_num converts text to the type and kind of the mold, without an internal read
  n(1) = to_num('64', 1_int32)
  n(2) = to_num('32', 1_int32)
  alpha = to_num('1.0e-2', 1.0_dp)
  print '(a, 2i4, a, es10.3)', 'parsed n =', n, ', alpha =', alpha
  print '(a, f8.3)', 'single precision:', to_num('0.125', 1.0_sp)

  ! split a tab separated record by hand and keep the fields in a stringlist
  start = 1
  do i = 1, len(record) + 1
    if (i > len(record)) then
      call keys%insert_at(bidx(1), record(start:i - 1))
    else if (record(i:i) == TAB) then
      call keys%insert_at(bidx(1), record(start:i - 1))
      start = i + 1
    end if
  end do
  print '(a, i0)', 'fields: ', keys%len()
  first = keys%get(fidx(1))
  print '(dt)', first
  first = keys%get(bidx(1))
  print '(dt)', first
end program stdlib_text_demo

digits in record: 7; classes of "a", TAB, ",": TTT
digit set: 0123456789
HEATFLOW
heatflow
Heat Flow Solver
cba
parsed n =  64  32, alpha = 1.000E-02
single precision:   0.125
fields: 3
nx=64
alpha=1.0e-2

program stdlib_sorting_demo
  use stdlib_kinds, only: dp, int32
  use stdlib_sorting, only: sort, ord_sort, sort_index, radix_sort, sort_adjoint, int_index
  use stdlib_selection, only: select, arg_select
  implicit none
  real(dp) :: t(8), w(8), median
  integer(int_index) :: order(8)
  integer(int32) :: id(8)
  integer :: k, kth
  integer :: idx(8)

  t = [310.5_dp, 295.0_dp, 301.2_dp, 295.0_dp, 330.9_dp, 288.4_dp, 301.2_dp, 305.0_dp]
  id = [1_int32, 2_int32, 3_int32, 4_int32, 5_int32, 6_int32, 7_int32, 8_int32]

  w = t
  call sort(w)                          ! introsort, not stable
  print '(a, 8f7.1)', 'sort            ', w
  w = t
  call sort(w, reverse=.true.)
  print '(a, 8f7.1)', 'sort, reverse   ', w

  w = t
  call ord_sort(w)                      ! stable merge sort, fast on partly ordered data
  print '(a, 8f7.1)', 'ord_sort        ', w

  w = t
  call sort_index(w, order)             ! sorted values plus the permutation
  print '(a, 8f7.1)', 'sort_index      ', w
  print '(a, 8i7)', 'permutation     ', order

  w = t
  call sort_adjoint(w, id)              ! carry a second array along
  print '(a, 8i7)', 'ids sorted by T ', id

  call radix_sort(w)                    ! integers and reals only
  print '(a, l1)', 'radix_sort keeps sorted data sorted: ', all(w(1:7) <= w(2:8))

  ! selection: the 4th smallest value without sorting everything
  w = t
  call select(w, 4, median)             ! select reorders its argument
  print '(a, f6.1)', '4th smallest    ', median
  idx = [(k, k = 1, 8)]
  call arg_select(t, idx, 4, kth)       ! arg_select leaves the data alone
  print '(a, i0, a, f6.1)', 'its index       ', kth, ', value ', t(kth)
end program stdlib_sorting_demo

sort              288.4  295.0  295.0  301.2  301.2  305.0  310.5  330.9
sort, reverse     330.9  310.5  305.0  301.2  301.2  295.0  295.0  288.4
ord_sort          288.4  295.0  295.0  301.2  301.2  305.0  310.5  330.9
sort_index        288.4  295.0  295.0  301.2  301.2  305.0  310.5  330.9
permutation           6      2      4      3      7      8      1      5
ids sorted by T       6      2      4      3      7      8      1      5
radix_sort keeps sorted data sorted: T
4th smallest     301.2
its index       7, value  301.2

program stdlib_containers_demo
  use stdlib_kinds, only: dp
  use stdlib_hashmaps, only: chaining_hashmap_type, open_hashmap_type
  use stdlib_hashmap_wrappers, only: key_type, set, get
  use stdlib_bitsets, only: bitset_large
  implicit none
  type(chaining_hashmap_type) :: registry
  type(open_hashmap_type) :: seen
  type(key_type), allocatable :: keys(:)
  class(*), allocatable :: value
  character(len=:), allocatable :: name
  logical :: conflict, exists
  type(bitset_large) :: fixed
  integer :: i, nx

  ! a map from solver names to relaxation factors; values are class(*)
  call registry%init(slots_bits=6)
  call registry%map_entry('jacobi', 1.0_dp, conflict)
  call registry%map_entry('ftcs', 0.25_dp, conflict)
  call registry%map_entry('jacobi', 2.0_dp, conflict)        ! duplicate key
  print '(a, l1, a, i0)', 'duplicate rejected: ', conflict, ', entries = ', registry%entries()

  call registry%get_other_data('ftcs', value, exists)
  select type (value)
  type is (real(dp))
    print '(a, f5.2, a, l1)', 'ftcs -> ', value, ' found: ', exists
  end select

  call registry%set_other_data('ftcs', 0.20_dp)
  call registry%get_all_keys(keys)
  print '(a, i0)', 'keys stored: ', size(keys)
  do i = 1, size(keys)
    call get(keys(i), name)
    print '(2a)', '  key ', name
  end do
  call registry%remove('jacobi', exists)
  print '(a, l1, a, i0)', 'removed: ', exists, ', entries now ', registry%entries()

  ! the open-addressing variant has the same API
  call seen%init()
  call seen%map_entry([3, 4], conflict=conflict)
  call seen%key_test([3, 4], exists)
  print '(a, l1)', 'integer key present: ', exists

  ! a bitset marking the fixed (boundary) cells of a 12-cell row (bit positions start at 0)
  nx = 12
  call fixed%init(nx)
  call fixed%set(0)
  call fixed%set(nx - 1)
  print '(a, i0, a, i0, a)', 'fixed cells: ', fixed%bit_count(), ' of ', fixed%bits(), ' bits'
  print '(a, l1, l1)', 'bit 0 set, bit 4 set: ', fixed%test(0), fixed%test(4)
  call fixed%not()
  print '(a, i0)', 'interior cells after not: ', fixed%bit_count()
end program stdlib_containers_demo

duplicate rejected: T, entries = 2
ftcs ->  0.25 found: T
keys stored: 2
  key jacobi
  key ftcs
removed: T, entries now 1
integer key present: T
fixed cells: 2 of 12 bits
bit 0 set, bit 4 set: TF
interior cells after not: 10

program stdlib_numerics_demo
  use stdlib_kinds, only: dp, int32
  use stdlib_math, only: linspace, arange, clip, is_close, all_close, diff, deg2rad
  use stdlib_stats, only: mean, var, median, corr
  use stdlib_random, only: random_seed
  use stdlib_stats_distribution_normal, only: rvs_normal, pdf_normal, cdf_normal
  use stdlib_quadrature, only: trapz, simps, gauss_legendre
  use stdlib_specialfunctions_gamma, only: gamma, log_gamma
  implicit none
  integer, parameter :: n = 9
  real(dp) :: x(n), y(n), noisy(1000), nodes(5), weights(5)
  real(dp) :: pi_, pair(n, 2), c(2, 2)
  integer :: previous

  ! grids and sequences
  x = linspace(0.0_dp, 1.0_dp, n)
  print '(a, 9f6.3)', 'linspace        ', x
  print '(a, 5i3)', 'arange(1, 10, 2)', arange(1, 10, 2)
  print '(a, 3f6.2)', 'clip            ', clip(-0.5_dp, 0.0_dp, 1.0_dp), clip(0.4_dp, 0.0_dp, 1.0_dp), clip(1.7_dp, 0.0_dp, 1.0_dp)
  print '(a, l1, l1)', 'is_close, <=0   ', is_close(0.1_dp + 0.2_dp, 0.3_dp), abs(0.1_dp + 0.2_dp - 0.3_dp) <= 0.0_dp
  print '(a, 4f6.3)', 'diff            ', diff(x(1:5))

  ! a temperature profile and its statistics
  y = 100.0_dp * sin(deg2rad(180.0_dp) * x)
  print '(a, f8.3, a, f8.3, a, f8.3)', 'mean ', mean(y), ' var ', var(y), ' median ', median(y)
  pair(:, 1) = x
  pair(:, 2) = 100.0_dp * x**2
  c = corr(pair, 1)                     ! 2 x 2 correlation matrix of the columns
  print '(a, f8.4)', 'corr(x, 100 x^2)', c(1, 2)

  ! integrals of the profile: trapezoid, Simpson, Gauss-Legendre
  print '(a, f10.5)', 'trapz           ', trapz(y, x)
  print '(a, f10.5)', 'simps           ', simps(y, x)
  print '(a, f10.5)', 'exact           ', 200.0_dp / acos(-1.0_dp)
  call gauss_legendre(nodes, weights, [0.0_dp, 1.0_dp])
  pi_ = acos(-1.0_dp)
  print '(a, f10.5)', 'gauss_legendre  ', sum(100.0_dp * sin(pi_ * nodes) * weights)

  ! reproducible noise: fix the seed, then draw a normal sample
  call random_seed(12345_int32, previous)
  noisy = rvs_normal(10.0_dp, 2.0_dp, 1000)
  print '(a, f8.3, a, f8.3)', 'sample mean ', mean(noisy), ', sample sd ', sqrt(var(noisy))
  print '(a, f8.5, a, f8.5)', 'pdf at mean ', pdf_normal(10.0_dp, 10.0_dp, 2.0_dp), ', cdf at mean ', cdf_normal(10.0_dp, 10.0_dp, 2.0_dp)

  ! special functions
  print '(a, f10.5, f10.5)', 'gamma(5), gamma(0.5)', gamma(5.0_dp), gamma(0.5_dp)
  print '(a, f10.5)', 'log_gamma(10)   ', log_gamma(10.0_dp)
  print '(a, l1)', 'all_close       ', all_close(x, linspace(0.0_dp, 1.0_dp, n))
end program stdlib_numerics_demo

linspace         0.000 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000
arange(1, 10, 2)  1  3  5  7  9
clip              0.00  0.40  1.00
is_close, <=0   TF
diff             0.125 0.125 0.125 0.125
mean   55.859 var 1489.702 median   70.711
corr(x, 100 x^2)  0.9622
trapz             62.84174
simps             63.67055
exact             63.66198
gauss_legendre    63.66198
sample mean    9.950, sample sd    2.054
pdf at mean  0.19947, cdf at mean  0.50000
gamma(5), gamma(0.5)  24.00000   1.77245
log_gamma(10)     12.80183
all_close       T

program stdlib_sparse_demo
  use stdlib_kinds, only: dp
  use stdlib_array, only: trueloc
  use stdlib_sparse, only: CSR_dp_type, from_ijv, spmv
  use stdlib_specialmatrices, only: tridiagonal_dp_type, tridiagonal, dense, spmv
  use stdlib_linalg_iterative_solvers, only: stdlib_solve_cg
  implicit none
  integer, parameter :: n = 8
  real(dp), parameter :: h = 1.0_dp / real(n + 1, dp)
  integer :: row(3 * n - 2), col(3 * n - 2), k, i
  real(dp) :: val(3 * n - 2), b(n), x(n), y(n), res
  type(CSR_dp_type) :: a
  type(tridiagonal_dp_type) :: t
  real(dp), allocatable :: hot(:)

  ! steady 1D heat flow with uniform source: -u'' = 1, u(0) = u(1) = 0
  ! finite differences give the tridiagonal matrix tridiag(-1, 2, -1) / h**2
  k = 0
  do i = 1, n
    if (i > 1) then
      k = k + 1; row(k) = i; col(k) = i - 1; val(k) = -1.0_dp / h**2
    end if
    k = k + 1; row(k) = i; col(k) = i; val(k) = 2.0_dp / h**2
    if (i < n) then
      k = k + 1; row(k) = i; col(k) = i + 1; val(k) = -1.0_dp / h**2
    end if
  end do
  call from_ijv(a, row, col, val)          ! coordinate triplets to compressed sparse row
  print '(a, i0, a, i0, a, i0)', 'CSR: ', a%nrows, ' rows, ', a%nnz, ' stored entries instead of ', n * n

  b = 1.0_dp
  x = 0.0_dp
  call stdlib_solve_cg(a, b, x, rtol=1.0e-12_dp, restart=.false.)   ! conjugate gradients
  print '(a, 8f7.4)', 'u (CG)       =', x

  y = 0.0_dp
  call spmv(a, x, y)                       ! y = A x
  res = maxval(abs(y - b))
  print '(a, es9.2)', 'max |A u - b| =', res

  ! the exact solution u(x) = x (1 - x) / 2 at the grid points
  print '(a, 8f7.4)', 'exact        =', [(0.5_dp * real(i, dp) * h * (1.0_dp - real(i, dp) * h), i = 1, n)]

  ! the same operator as a tridiagonal special matrix: three vectors, no sparse structure needed
  t = tridiagonal(spread(-1.0_dp / h**2, 1, n - 1), spread(2.0_dp / h**2, 1, n), spread(-1.0_dp / h**2, 1, n - 1))
  y = 0.0_dp
  call spmv(t, x, y)
  print '(a, es9.2)', 'special matrix agrees: ', maxval(abs(y - b))
  print '(a, i0, a, i0)', 'dense() shape ', size(dense(t), 1), ' x ', size(dense(t), 2)

  ! trueloc turns a logical mask into the index list of its .true. elements
  hot = x(trueloc(x > 0.1_dp))
  print '(a, i0, a, 4f7.4)', 'cells above 0.1: ', size(hot), ' ', hot
end program stdlib_sparse_demo

CSR: 8 rows, 22 stored entries instead of 64
u (CG)       = 0.0494 0.0864 0.1111 0.1235 0.1235 0.1111 0.0864 0.0494
max |A u - b| = 0.00E+00
exact        = 0.0494 0.0864 0.1111 0.1235 0.1235 0.1111 0.0864 0.0494
special matrix agrees:  0.00E+00
dense() shape 8 x 8
cells above 0.1: 4  0.1111 0.1235 0.1235 0.1111

/usr/bin/ld: warning: .._stdlib-src_stdlib-fpm_src_stdlib_linalg_iterative_solvers_cg.f90.o: requires executable stack (because the .note.GNU-stack section is executable)

program stdlib_intrinsics_hash
  use stdlib_kinds, only: sp, dp
  use stdlib_intrinsics, only: stdlib_sum, stdlib_sum_kahan
  use stdlib_hash_32bit, only: fnv_1a_hash
  use stdlib_hash_64bit, only: fnv_1a_hash_64 => fnv_1a_hash
  implicit none
  integer, parameter :: n = 2000000
  real(sp), allocatable :: x(:)
  real(dp) :: reference
  integer :: i

  ! many small single-precision values: the naive running sum loses the small terms
  allocate (x(n))
  x = 1.0e-3_sp
  x(1) = 1.0e4_sp
  reference = 1.0e4_dp + real(n - 1, dp) * real(1.0e-3_sp, dp)
  print '(a, f14.4)', 'reference (dp)   ', reference
  print '(a, f14.4)', 'sum              ', sum(x)
  print '(a, f14.4)', 'stdlib_sum       ', stdlib_sum(x)
  print '(a, f14.4)', 'stdlib_sum_kahan ', stdlib_sum_kahan(x)

  ! non-cryptographic hashes of integer vectors and bytes, e.g. for cache keys
  print '(a, i0)', 'fnv_1a 32-bit    ', fnv_1a_hash([(i, i = 1, 8)])
  print '(a, i0)', 'fnv_1a 64-bit    ', fnv_1a_hash_64([(i, i = 1, 8)])
end program stdlib_intrinsics_hash

reference (dp)       11999.9991
sum                  11953.1240
stdlib_sum           11998.8975
stdlib_sum_kahan     11999.9990
fnv_1a 32-bit    -1706177971
fnv_1a 64-bit    7244529863066876845

reference (dp)       11999.9991
sum                  11999.9990
stdlib_sum           11998.8975
stdlib_sum_kahan     11999.9990
fnv_1a 32-bit    -1706177971
fnv_1a 64-bit    7244529863066876845

program stdlib_system_demo
  use stdlib_error, only: check, state_type
  use stdlib_system, only: make_directory_all, remove_directory, delete_file, exists, is_directory, &
                           join_path, base_name, dir_name, run, process_type, is_windows, operator(/)
  use stdlib_ansi, only: fg_color_green, fg_color_red, style_bold, style_reset, ansi_code, &
                         operator(//), operator(+)
  implicit none
  type(state_type) :: err
  type(process_type) :: p
  character(len=:), allocatable :: out_dir, file
  type(ansi_code) :: ok_style
  integer :: u

  out_dir = 'results' / 'run01'
  call make_directory_all(out_dir, err)
  print '(a, l1, a, l1)', 'created: ', err%ok(), ', is_directory: ', is_directory(out_dir)

  file = join_path(out_dir, 'heat.csv')
  open(newunit=u, file=file, status='replace', action='write')
  write (u, '(a)') '0.0,1.0'
  close(u)
  print '(2a)', 'file    : ', file
  print '(2a)', 'base    : ', base_name(file)
  print '(2a)', 'dir     : ', dir_name(file)
  print '(a, l1)', 'exists  : ', exists(file)

  ! run an external program and capture its standard output
  if (.not. is_windows()) then
    p = run('wc -l ' // file, want_stdout=.true.)
    print '(a, i0, 2a)', 'wc exit code ', p%exit_code, ': ', trim(adjustl(p%stdout))
  end if

  ! state-based error handling instead of error stop
  call delete_file(file, err)
  print '(a, l1)', 'first delete ok: ', err%ok()
  call delete_file(file, err)                    ! the file is gone: reported as ok, not as an error
  print '(a, l1, 2a)', 'second delete ok: ', err%ok(), ' -> ', err%print()
  call delete_file(out_dir, err)                 ! a directory cannot be deleted as a file
  print '(a, l1, 2a)', 'delete directory failed: ', err%error(), ' -> ', err%print()
  call remove_directory(out_dir, err)
  call remove_directory('results', err)
  print '(a, l1)', 'cleaned up: ', err%ok()

  ! check(): a failed condition with warn=.true. prints and continues
  call check(1 + 1 == 3, msg='arithmetic is broken', warn=.true.)

  ! ANSI escape sequences for coloured terminal output
  ok_style = fg_color_green + style_bold
  print '(a)', ok_style // 'PASS' // style_reset
  print '(a)', fg_color_red // 'FAIL' // style_reset
end program stdlib_system_demo

created: T, is_directory: T
file    : results/run01/heat.csv
base    : heat.csv
dir     : results/run01
exists  : T
wc exit code 0: 1 results/run01/heat.csv
first delete ok: T
second delete ok: T -> Success!
delete directory failed: T -> Error encountered: Filesystem Error: Cannot delete results/run01 - is a directory
 arithmetic is broken
cleaned up: T
ESC[0;1;32mPASSESC[0m
ESC[0;31mFAILESC[0m