Numerical: Linear Algebra (BLAS/LAPACK)

Code

program linalg
  implicit none
  real(8), parameter :: alpha = 1.0, beta = 0.0
  integer, parameter :: m=3, n=2, k=4
  real(8) :: A(m,k), B(k,n), C(m,n)
  integer :: i, j

  ! Initialize matrices
  A = reshape([(real(i,8), i=1,m*k)], [m,k])
  B = reshape([(real(i,8), i=1,k*n)], [k,n])

  ! DGEMM: C = alpha * A * B + beta * C
  ! Matrix multiply
  call dgemm('N','N', m,n,k, alpha, A,m, B,k, beta, C,m)

  print *, 'C(1,:):', C(1,:)
  print *, 'C(2,:):', C(2,:)

  ! LAPACK example: solve Ax = b
  real(8) :: A2(3,3), b2(3)
  integer :: ipiv(3), info
  A2 = reshape([4.0,2.0,1.0, 2.0,3.0,1.0, 1.0,1.0,2.0], [3,3])
  b2 = [10.0, 7.0, 5.0]
  call dgesv(3, 1, A2, 3, ipiv, b2, 3, info)
  if (info == 0) print *, 'Solution:', b2
end program

Explanation

DGEMM (double-precision general matrix multiply) is the BLAS workhorse — uses optimized CPU kernels. LAPACK's DGESV solves Ax=b via LU decomposition. Leading dimensions (m, k) specify memory layout. Compile with -llapack -lblas. Always check info=0 (success). For production, use MKL or OpenBLAS for maximum performance.

More Fortran Snippets