pure function mat_mat_rel(A, B, transA, transB, option) result(C)
real(rk), intent(in), contiguous :: A(:,:), B(:,:)
real(rk), allocatable :: C(:,:)
character(*), intent(in), optional :: option
logical, intent(in), optional :: transA, transB
if (present(transA) .and. present(transB)) then
if (.not.transA .and. .not.transB) then
! AB
allocate(C(size(A,1), size(B,2)), source=0.0_rk)
call mat_mat_rel_AB(A, B, C, option)
else if (transA .and. transB) then
! ATBT
allocate(C(size(A,2), size(B,1)), source=0.0_rk)
call mat_mat_rel_ATBT(A, B, C, option)
else if (transA .and. .not.transB) then
! ATB
allocate(C(size(A,2), size(B,2)), source=0.0_rk)
call mat_mat_rel_ATB(A, B, C, option)
else if (.not.transA .and. transB) then
! ABT
allocate(C(size(A,1), size(B,1)), source=0.0_rk)
call mat_mat_rel_ABT(A, B, C, option)
end if
else if (present(transA) .or. present(transB)) then
if (present(transA)) then
if (transA) then
! ATB
allocate(C(size(A,2), size(B,2)), source=0.0_rk)
call mat_mat_rel_ATB(A, B, C, option)
else if (.not.transA) then
! ABT
allocate(C(size(A,1), size(B,1)), source=0.0_rk)
call mat_mat_rel_ABT(A, B, C, option)
end if
else if (present(transB)) then
if (transB) then
! ABT
allocate(C(size(A,1), size(B,1)), source=0.0_rk)
call mat_mat_rel_ABT(A, B, C, option)
else if (.not.transB) then
! ATB
allocate(C(size(A,2), size(B,2)), source=0.0_rk)
call mat_mat_rel_ATB(A, B, C, option)
end if
end if
else if (.not.present(transA) .and. .not.present(transB)) then
! AB
allocate(C(size(A,1), size(B,2)), source=0.0_rk)
call mat_mat_rel_AB(A, B, C, option)
end if
end function mat_mat_rel