Numerik: Lineare Algebra (BLAS/LAPACK)

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