-
Notifications
You must be signed in to change notification settings - Fork 0
/
insitu.f90
73 lines (53 loc) · 1.64 KB
/
insitu.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
module insitu
implicit none
type, public :: insitu_Array
real(4), pointer :: storage(:)
integer, dimension(2) :: shp
real(4), pointer :: A(:,:) => null()
contains
procedure :: init => insitu_Array_init
procedure :: t => insitu_Array_transpose
end type insitu_Array
contains
subroutine insitu_Array_init(this)
use, intrinsic :: iso_c_binding, only : c_loc, c_f_pointer
class(insitu_Array), intent(inout) :: this
call c_f_pointer(c_loc(this%storage), this%A, this%shp)
end subroutine insitu_Array_init
subroutine insitu_Array_transpose(this)
use, intrinsic :: iso_c_binding
class(insitu_Array), intent(inout) :: this
call mipt(this%shp(2), this%shp(1), this%storage)
nullify(this%A)
call c_f_pointer(c_loc(this%storage), this%A, [this%shp(2), this%shp(1)])
this%shp = shape(this%A)
end subroutine insitu_Array_transpose
subroutine mipt(n, m, A)
integer, intent(in) :: n, m
real(4), pointer, intent(in out) :: A(:)
integer :: cnt, q, k, i, j
real(4) :: tmp
logical(kind=1), allocatable :: flg(:)
allocate(flg(m*n))
flg = .false.
q = m * n - 1
do cnt = 1, m*n-2
k = m * cnt - q * (cnt / n)
do while (k > cnt)
k = m * k - q * (k / n)
end do
if (k == cnt) then
i = k
tmp = A(i+1)
do
j = m * i - q * (i / n)
if (j == k) exit
A(i+1) = A(j+1)
i = j
end do
A(i+1) = tmp
end if
end do
deallocate(flg)
end subroutine mipt
end module insitu