Create all data needed to quickly map between nd index and 2d index.
Type  Intent  Optional  Attributes  Name  

type(nd_to_2d_mapping),  intent(out)  ::  map 
index mapping data 

integer,  intent(in),  DIMENSION(:)  ::  dims 
nd sizes which ndindices map to first matrix index and in which order which ndindices map to second matrix index and in which order 

integer,  intent(in),  DIMENSION(:)  ::  map1_2d 
nd sizes which ndindices map to first matrix index and in which order which ndindices map to second matrix index and in which order 

integer,  intent(in),  DIMENSION(:)  ::  map2_2d 
nd sizes which ndindices map to first matrix index and in which order which ndindices map to second matrix index and in which order 

integer,  intent(in),  optional  ::  base 
base index (1 for Fortranstyle, 0 for Cstyle, default is 1) 

logical,  intent(in),  optional  ::  col_major 
whether index should be column major order (.TRUE. for Fortranstyle, .FALSE. for Cstyle, default is .TRUE.). 
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 ndindices map to first matrix index and in which order !! which ndindices map to second matrix index and in which order INTEGER, INTENT(IN), OPTIONAL :: base !! base index (1 for Fortranstyle, 0 for Cstyle, default is 1) LOGICAL, INTENT(IN), OPTIONAL :: col_major !! whether index should be column major order (.TRUE. for Fortranstyle, .FALSE. for Cstyle, 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