21 real(kind=
kreal),
allocatable :: ajad(:)
22 integer(kind=kint),
allocatable :: JAJAD(:)
23 integer(kind=kint),
allocatable :: JADORD(:)
24 integer(kind=kint),
allocatable :: IAJAD(:)
25 integer(kind=kint) :: MJAD
26 real(kind=
kreal),
allocatable :: wp1(:), wp2(:), wp3(:), wp4(:)
27 integer(kind=kint) :: INITIALIZED = 0
33 allocate(wp1(hecmat%NP), wp2(hecmat%NP), wp3(hecmat%NP), wp4(hecmat%NP))
34 allocate(ajad((hecmat%NPL+hecmat%NPU)*16))
35 allocate(jajad(hecmat%NPL+hecmat%NPU))
36 allocate(jadord(hecmat%NP))
37 allocate(iajad(hecmat%NP+1))
38 call repack(hecmat%N, hecmat, mjad, ajad, jajad, iajad, jadord)
47 deallocate(wp1,wp2,wp3,wp4)
59 real(kind=
kreal),
intent(in) :: x(:)
60 real(kind=
kreal),
intent(out) :: y(:)
61 real(kind=
kreal),
intent(inout) :: commtime
62 real(kind=
kreal) :: start_time, end_time
63 real(kind=
kreal),
pointer :: d(:)
64 integer(kind=kint) :: i
65 real(kind=
kreal) :: x1, x2, x3, x4
70 commtime = commtime + end_time - start_time
81 y(4*i-3)= d(16*i-15)*x1 + d(16*i-14)*x2 + d(16*i-13)*x3 + d(16*i-12)*x4
82 y(4*i-2)= d(16*i-11)*x1 + d(16*i-10)*x2 + d(16*i- 9)*x3 + d(16*i- 8)*x4
83 y(4*i-1)= d(16*i- 7)*x1 + d(16*i- 6)*x2 + d(16*i- 5)*x3 + d(16*i- 4)*x4
84 y(4*i )= d(16*i- 3)*x1 + d(16*i- 2)*x2 + d(16*i- 1)*x3 + d(16*i- 0)*x4
88 call matjad(hecmat%N, mjad, iajad, jajad, ajad, jadord, x, y, wp1, wp2, wp3, wp4)
91 subroutine repack(N, hecMAT, MJAD, AJAD, JAJAD, IAJAD, JADORD)
96 integer(kind = kint) :: n, mjad
97 real(kind =
kreal),
dimension(*) :: ajad
98 integer(kind = kint),
dimension(*) :: jajad
99 integer(kind = kint),
dimension(*) :: iajad
100 integer(kind = kint),
dimension(*) :: jadord
102 integer(kind = kint) :: ijad, maxnz, minnz
103 integer(kind = kint) :: i, j, js, je, in, jc
104 integer(kind = kint),
allocatable :: len(:), lenz(:), jadreord(:)
107 allocate(jadreord(n))
109 len(i)= hecmat%indexL(i) - hecmat%indexL(i-1) &
110 & + hecmat%indexU(i) - hecmat%indexU(i-1)
112 maxnz=maxval(len(1:n))
113 minnz=minval(len(1:n))
115 allocate(lenz(0:mjad))
118 lenz(len(i))=lenz(len(i))+1
120 do i=maxnz-1,minnz,-1
121 lenz(i)=lenz(i)+lenz(i+1)
124 jadord(i)=lenz(len(i))
125 lenz(len(i))=lenz(len(i))-1
128 jadreord(jadord(i))=i
131 lenz(len(jadreord(i)))=i
134 lenz(i)=max(lenz(i+1),lenz(i))
138 iajad(i+1)=iajad(i)+lenz(i)
143 js= hecmat%indexL(i-1) + 1
144 je= hecmat%indexL(i )
147 len(ijad)=len(ijad)+1
148 jc=iajad(len(ijad))+ijad-1
149 ajad(jc*16-15:jc*16) = hecmat%AL(16*j-15:16*j)
155 js= hecmat%indexU(i-1) + 1
156 je= hecmat%indexU(i )
159 len(ijad)=len(ijad)+1
160 jc=iajad(len(ijad))+ijad-1
161 ajad(jc*16-15:jc*16) = hecmat%AU(16*j-15:16*j)
168 end subroutine repack
170 subroutine matjad(N, MJAD, IAJAD, JAJAD, AJAD, JADORD, X, Y, W1, W2, W3, W4)
172 integer(kind=kint) :: N, MJAD
173 integer(kind=kint) :: IAJAD(*), JAJAD(*), JADORD(*)
174 real(kind=
kreal) :: ajad(*), x(*), y(*), w1(*), w2(*), w3(*), w4(*)
176 integer(kind=kint) :: I, K, NZ, IXX
177 real(kind=
kreal) :: x1, x2, x3, x4
191 do k=iajad(nz),iajad(nz+1)-1
197 w1(ixx)=w1(ixx) + ajad(k*16-15)*x1 + ajad(k*16-14)*x2 + ajad(k*16-13)*x3 + ajad(k*16-12)*x4
198 w2(ixx)=w2(ixx) + ajad(k*16-11)*x1 + ajad(k*16-10)*x2 + ajad(k*16- 9)*x3 + ajad(k*16- 8)*x4
199 w3(ixx)=w3(ixx) + ajad(k*16- 7)*x1 + ajad(k*16- 6)*x2 + ajad(k*16- 5)*x3 + ajad(k*16- 4)*x4
200 w4(ixx)=w4(ixx) + ajad(k*16- 3)*x1 + ajad(k*16- 2)*x2 + ajad(k*16- 1)*x3 + ajad(k*16- 0)*x4
207 y(4*i-3)=y(4*i-3)+w1(jadord(i))
208 y(4*i-2)=y(4*i-2)+w2(jadord(i))
209 y(4*i-1)=y(4*i-1)+w3(jadord(i))
210 y(4*i )=y(4*i )+w4(jadord(i))