create_nd_to_2d_mapping Subroutine

   SUBROUTINE create_nd_to_2d_mapping(map, dims, map1_2d, map2_2d, base, col_major)
      !! Create all data needed to quickly map between nd index and 2d index.

      TYPE(nd_to_2d_mapping), INTENT(OUT)                :: map
         !! index mapping data
      INTEGER, DIMENSION(:), INTENT(IN)                  :: dims, map1_2d, map2_2d
         !! nd sizes
         !! which nd-indices map to first matrix index and in which order
         !! which nd-indices map to second matrix index and in which order
      INTEGER, INTENT(IN), OPTIONAL                      :: base
         !! base index (1 for Fortran-style, 0 for C-style, default is 1)
      LOGICAL, INTENT(IN), OPTIONAL                      :: col_major
         !! whether index should be column major order (.TRUE. for Fortran-style, .FALSE. for C-style, default is .TRUE.).

      INTEGER                                            :: i

      IF (PRESENT(col_major)) THEN
         map%col_major = col_major
      ELSE
         map%col_major = .TRUE.
      END IF

      IF (PRESENT(base)) THEN
         map%base = base
      ELSE
         map%base = 1
      END IF

      map%ndim1_2d = SIZE(map1_2d)
      map%ndim2_2d = SIZE(map2_2d)
      map%ndim_nd = SIZE(dims)

      CALL allocate_any(map%map1_2d, source=map1_2d)
      CALL allocate_any(map%map2_2d, source=map2_2d)
      CALL allocate_any(map%dims_nd, source=dims)
      CALL allocate_any(map%dims1_2d, source=dims(map1_2d))
      CALL allocate_any(map%dims2_2d, source=dims(map2_2d))

      ALLOCATE (map%map_nd(map%ndim_nd))
      map%map_nd(map1_2d) = (/(i, i=1, SIZE(map1_2d))/)
      map%map_nd(map2_2d) = (/(i + SIZE(map1_2d), i=1, SIZE(map2_2d))/)

      map%dims_2d = [PRODUCT(INT(map%dims1_2d, KIND=int_8)), PRODUCT(INT(map%dims2_2d, KIND=int_8))]

   END SUBROUTINE create_nd_to_2d_mapping