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