symmetric-tridiagonal-eigensolver/stack_mod.f90 at master · xinye83/symmetric-tridiagonal-eigensolver · GitHub

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
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
module stack_mod
  ! a stack data structure for the main program

  implicit none

  type stack
     integer             :: len
     real(8),allocatable :: p(:,:)
     integer,allocatable :: c(:,:)
  end type stack

contains

  subroutine stack_copy(S1,S2)
    ! copy stack S1 to S2

    type(stack) :: S1,S2

    S2%len=S1%len
    allocate(S2%p(S2%len,2),S2%c(S2%len,2))

    S2%p=S1%p
    S2%c=S1%c
  end subroutine stack_copy

  subroutine stack_allo(S)
    ! allocate p and c in stack S
    type(stack) :: S
    allocate(S%p(S%len,2),S%c(S%len,2))
  end subroutine stack_allo

  subroutine stack_deal(S)
    ! deallocate stack S
    type(stack) :: S
    deallocate(S%p,S%c)
  end subroutine stack_deal

  subroutine stack_pop(S,p1,p2,c1,c2)
    ! pop one entry from stack

    type(stack) :: S
    real(8)     :: p1,p2
    integer     :: c1,c2

    type(stack) :: temp

    if (S%len==1) then
       p1=S%p(1,1)
       p2=S%p(1,2)

       c1=S%c(1,1)
       c2=S%c(1,2)

       call stack_deal(S)
       S%len=0
    else
       call stack_copy(S,temp)
       call stack_deal(S)

       S%len=temp%len-1
       call stack_allo(S)

       S%p=temp%p(2:temp%len,:)
       S%c=temp%c(2:temp%len,:)

       p1=temp%p(1,1)
       p2=temp%p(1,2)

       c1=temp%c(1,1)
       c2=temp%c(1,2)

       call stack_deal(temp)
    end if

  end subroutine stack_pop

  subroutine stack_push(S,p1,p2,c1,c2)
    ! push one entry into stack

    type(stack) :: S
    real(8)     :: p1,p2
    integer     :: c1,c2

    type(stack) :: temp

    if (S%len==0) then
       S%len=1
       call stack_allo(S)
       S%p(1,1)=p1
       S%p(1,2)=p2

       S%c(1,1)=c1
       S%c(1,2)=c2
    else
       call stack_copy(S,temp)
       call stack_deal(S)

       S%len=temp%len+1
       call stack_allo(S)

       S%p(1,1)=p1
       S%p(1,2)=p2
       S%p(2:S%len,:)=temp%p

       S%c(1,1)=c1
       S%c(1,2)=c2
       S%c(2:S%len,:)=temp%c

       call stack_deal(temp)
    end if

  end subroutine stack_push

  subroutine stack_pick(S,k,p1,p2,c1,c2)
    ! extract k-th entry of stack S
    type(stack) :: S
    integer     :: k,c1,c2
    real(8)     :: p1,p2

    p1=S%p(k,1)
    p2=S%p(k,2)

    c1=S%c(k,1)
    c2=S%c(k,2)

  end subroutine stack_pick

end module stack_mod