Skip to main content

Basic Linear Algebra Subprograms

Install BLAS and LAPACK libraries:

apt install blaslib-dev lapacklib-dev

# compile program
gfortran filename.f90 -lblas

# in case libs are not in system path
gfortran filename.f90 -L/usr/lib -lblas

Scaler product (SDOT)​

src/22_blas_sdot.f90
! SDOT: This function performs a vector-vector operation of computing a scalar
! product of two single-precision real vectors x and y.
! https://netlib.org/blas/#_level_1
! Compile: gfortran filename.f90 -lblas
! If BLAS is not in system path: gfortran filename.f90 -L/usr/lib -lblas
PROGRAM MAIN
    IMPLICIT NONE

    INTEGER :: n = 3 ! length of vectors
    REAL, DIMENSION(:), ALLOCATABLE :: x
    REAL, DIMENSION(:), ALLOCATABLE :: y
    INTEGER :: incx = 1, incy = 1 
    
    ! incx and incy can be used to manipulate data when both x and y are stored
    ! in the same array. See the alternate array z
    REAL z(6)
    INTEGER :: incxz= 2, incyz = 1

    REAL res, SDOT

    EXTERNAL SDOT

    ALLOCATE(x(n))
    ALLOCATE(y(n))

    x(1) = 2.0
    x(2) = 1.0
    x(3) = 0.0

    y(1) = 1.0
    y(2) = 3.0
    y(3) = 5.0

    z(1) = x(1)
    z(2) = y(1)
    z(3) = x(2)
    z(4) = y(2)
    z(5) = x(3)
    z(6) = y(3)

    res = SDOT(n, x, incx, y, incy)

    PRINT *, "Result = ", res
    PRINT *, "Alt. result = ", SDOT(n, z, incxz, z, incyz)

    DEALLOCATE(x)
    DEALLOCATE(y)
END

Matrix multiplication (DGEMM)​

DGEMM calculates the product of double precision matrices:

C⇐αA∗B+βCC \Leftarrow \alpha A \ast B + \beta C
src/22_blas_dgemm.f90
PROGRAM MAIN
    IMPLICIT NONE

    DOUBLE PRECISION ALPHA, BETA
    INTEGER M, K, N, I, J
    PARAMETER (M=2, K=3, N=2)
    DOUBLE PRECISION A(M,K), B(K,N), C(M,N)

    ALPHA = 1.0 
    BETA = 0.0

    ! initialize matrices
    A(1,1) = 2.0
    A(1,2) = 3.0
    A(1,3) = 9.0
    A(2,1) = 0.0
    A(2,2) = 4.0
    A(2,3) = 1.0

    B(1,1) = 4.0
    B(1,2) = 1.0
    B(2,1) = 4.0
    B(2,2) = 5.0
    B(3,1) = 9.0
    B(3,2) = 6.0

    ! optionally initialize C matrix with zeros
    DO I = 1, M
        DO J = 1, N
          C(I,J) = 0.0
        END DO
    END DO

    CALL DGEMM('N','N', M, N, K, ALPHA, A, M, B, K, BETA, C, M)

    ! print matrix C
    DO I = 1, M
        DO J = 1, N
          PRINT *, "C(", I, ",", J, ") =", C(I,J)
        END DO
    END DO

END

Resources​