shoot/parder.f90 at main · sscollis/shoot · 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
129
130
131
!=============================================================================!
      subroutine parder( neq, Q, y, dQ )
!=============================================================================!
      implicit none

      integer :: neq
      real    :: y
      complex :: Q(neq), dQ(neq)

!.... local variables

      complex :: Eh(neq,neq), Fh(neq,neq), Ehinv(neq,neq)

!.... local variables for Lapack routines

      integer :: ipiv(neq), info

      real :: cpu
      real, external :: second
!=============================================================================!
      call parallel( y, Eh, Fh, 0 )

!.... Multiply Fh by Eh^{-1}

      call inverse( neq, Eh, Ehinv)
      Fh = matmul( Ehinv, Fh )

!     call CGETRF( neq, neq, Eh, neq, ipiv, info)
!     if (info.ne.0) write(*,*) 'CGETRF: ', info
!     call CGETRS('N', neq, neq, Eh, neq, ipiv, Fh, neq, info)
!     if (info.ne.0) write(*,*) 'CGETRS: ',info

!.... Form the derivative

      dQ = -matmul( Fh, Q )

      return
      end

!=============================================================================!
      subroutine inverse( neq, Eh, Ehinv )
!
!     Hardwired inverse for Eh.  Note that this takes advantage of the
!     sparce structure of Eh.  If, for some wierd reason, the definition
!     of Eh changes then this routine WILL produce incorrect results.
!
!=============================================================================!
      implicit none
      integer :: neq
      complex :: Eh(neq,neq), Ehinv(neq,neq), fact
      real, parameter :: zero = 0.0, one = 1.0
!=============================================================================!
      fact = one / ( Eh(3,1) * Eh(4,7) )
      Ehinv(1,1) = ( Eh(3,7) * Eh(4,3) - Eh(3,3) * Eh(4,7) ) / &
                   ( Eh(1,3) * Eh(3,1) * Eh(4,7) )
      Ehinv(1,2) = zero
      Ehinv(1,3) = one / Eh(3,1)
      Ehinv(1,4) = -Eh(3,7) * fact
      Ehinv(1,5) = zero
      Ehinv(1,6) = -Eh(3,2) / Eh(3,1)
      Ehinv(1,7) = ( Eh(3,7) * Eh(4,4) - Eh(3,4) * Eh(4,7) ) * fact
      Ehinv(1,8) = ( Eh(3,7) * Eh(4,5) - Eh(3,5) * Eh(4,7) ) * fact

      Ehinv(2,1) = zero
      Ehinv(2,2) = zero
      Ehinv(2,3) = zero
      Ehinv(2,4) = zero
      Ehinv(2,5) = zero
      Ehinv(2,6) = one
      Ehinv(2,7) = zero
      Ehinv(2,8) = zero

      Ehinv(3,1) = one / Eh(1,3)
      Ehinv(3,2) = zero
      Ehinv(3,3) = zero
      Ehinv(3,4) = zero
      Ehinv(3,5) = zero
      Ehinv(3,6) = zero
      Ehinv(3,7) = zero
      Ehinv(3,8) = zero

      Ehinv(4,1) = zero
      Ehinv(4,2) = zero
      Ehinv(4,3) = zero
      Ehinv(4,4) = zero
      Ehinv(4,5) = zero
      Ehinv(4,6) = zero
      Ehinv(4,7) = one
      Ehinv(4,8) = zero

      Ehinv(5,1) = zero
      Ehinv(5,2) = zero
      Ehinv(5,3) = zero
      Ehinv(5,4) = zero
      Ehinv(5,5) = zero
      Ehinv(5,6) = zero
      Ehinv(5,7) = zero
      Ehinv(5,8) = one

      fact = one / Eh(2,6)
      Ehinv(6,1) = -Eh(2,3) / ( Eh(1,3) * Eh(2,6) )
      Ehinv(6,2) = fact
      Ehinv(6,3) = zero
      Ehinv(6,4) = zero
      Ehinv(6,5) = zero
      Ehinv(6,6) = -Eh(2,2) * fact
      Ehinv(6,7) = zero
      Ehinv(6,8) = -Eh(2,5) * fact

      fact = one / Eh(4,7)
      Ehinv(7,1) = -Eh(4,3) / ( Eh(1,3) * Eh(4,7) )
      Ehinv(7,2) = zero
      Ehinv(7,3) = zero
      Ehinv(7,4) = fact
      Ehinv(7,5) = zero
      Ehinv(7,6) = zero
      Ehinv(7,7) = -Eh(4,4) * fact
      Ehinv(7,8) = -Eh(4,5) * fact

      fact = one / Eh(5,8)
      Ehinv(8,1) = -Eh(5,3) / ( Eh(1,3) * Eh(5,8) )
      Ehinv(8,2) = zero
      Ehinv(8,3) = zero
      Ehinv(8,4) = zero
      Ehinv(8,5) = fact
      Ehinv(8,6) = -Eh(5,2) * fact
      Ehinv(8,7) = -Eh(5,4) * fact
      Ehinv(8,8) = -Eh(5,5) * fact

      return
      end