1!======================================================================!
2!                                                                      !
3!======================================================================!
4SUBROUTINE D6DOT(T,A,B,N)
5  IMPLICIT NONE
6
7  INTEGER :: jj
8  INTEGER :: l
9  INTEGER :: N
10  DOUBLE PRECISION :: T(9)
11  DOUBLE PRECISION :: A(9,*)
12  DOUBLE PRECISION :: B(9,*)
13  !----------------------------------------------------------------------
14  !
15  !      spdot1 performs inner product of sparse vectors
16  !
17  !
18  !      #coded by t.arakawa of RIST on 040510
19  !
20  !----------------------------------------------------------------------
21  DO l = 1, 9
22     T(l) = 0.0D0
23  ENDDO
24  DO jj = 1, N
25     T(1) = T(1) + A(1,jj)*B(1,jj) + A(4,jj)*B(4,jj) + A(7,jj)*B(7,jj)
26     T(2) = T(2) + A(2,jj)*B(1,jj) + A(5,jj)*B(4,jj) + A(8,jj)*B(7,jj)
27     T(3) = T(3) + A(3,jj)*B(1,jj) + A(6,jj)*B(4,jj) + A(9,jj)*B(7,jj)
28     T(4) = T(4) + A(1,jj)*B(2,jj) + A(4,jj)*B(5,jj) + A(7,jj)*B(8,jj)
29     T(5) = T(5) + A(2,jj)*B(2,jj) + A(5,jj)*B(5,jj) + A(8,jj)*B(8,jj)
30     T(6) = T(6) + A(3,jj)*B(2,jj) + A(6,jj)*B(5,jj) + A(9,jj)*B(8,jj)
31     T(7) = T(7) + A(1,jj)*B(3,jj) + A(4,jj)*B(6,jj) + A(7,jj)*B(9,jj)
32     T(8) = T(8) + A(2,jj)*B(3,jj) + A(5,jj)*B(6,jj) + A(8,jj)*B(9,jj)
33     T(9) = T(9) + A(3,jj)*B(3,jj) + A(6,jj)*B(6,jj) + A(9,jj)*B(9,jj)
34  ENDDO
35END SUBROUTINE D6DOT
36!======================================================================!
37!                                                                      !
38!======================================================================!
39SUBROUTINE D6DOTL(T,A,B,N)
40  IMPLICIT NONE
41
42  INTEGER :: jj
43  INTEGER :: l
44  INTEGER :: N
45  DOUBLE PRECISION :: T(6)
46  DOUBLE PRECISION :: A(9,*)
47  DOUBLE PRECISION :: B(9,*)
48  !----------------------------------------------------------------------
49  !
50  !      spdot1 performs inner product of sparse vectors
51  !
52  !
53  !      #coded by t.arakawa of RIST on 040510
54  !
55  !----------------------------------------------------------------------
56  !$dir max_trips(6)
57  DO l = 1, 6
58     T(l) = 0.0D0
59  ENDDO
60  DO jj = 1, N
61     T(1) = T(1) + A(1,jj)*B(1,jj) + A(4,jj)*B(4,jj) + A(7,jj)*B(7,jj)
62     T(2) = T(2) + A(2,jj)*B(1,jj) + A(5,jj)*B(4,jj) + A(8,jj)*B(7,jj)
63     T(3) = T(3) + A(2,jj)*B(2,jj) + A(5,jj)*B(5,jj) + A(8,jj)*B(8,jj)
64     T(4) = T(4) + A(3,jj)*B(1,jj) + A(6,jj)*B(4,jj) + A(9,jj)*B(7,jj)
65     T(5) = T(5) + A(3,jj)*B(2,jj) + A(6,jj)*B(5,jj) + A(9,jj)*B(8,jj)
66     T(6) = T(6) + A(3,jj)*B(3,jj) + A(6,jj)*B(6,jj) + A(9,jj)*B(9,jj)
67  ENDDO
68END SUBROUTINE D6DOTL
69!======================================================================!
70!                                                                      !
71!======================================================================!
72SUBROUTINE D6SDOT(Wi,A,B,N)
73  IMPLICIT NONE
74
75  INTEGER :: jj
76  INTEGER :: N
77  DOUBLE PRECISION :: Wi(3)
78  DOUBLE PRECISION :: A(3,*)
79  DOUBLE PRECISION :: B(9,*)
80  !----------------------------------------------------------------------
81  !
82  !      spdot1 performs inner product of sparse vectors
83  !
84  !
85  !      #coded by t.arakawa of RIST on 040510
86  !
87  !----------------------------------------------------------------------
88  DO jj = 1, N
89     Wi(1) = Wi(1) - A(1,jj)*B(1,jj) - A(2,jj)*B(4,jj) - A(3,jj)*B(7,jj)
90     Wi(2) = Wi(2) - A(1,jj)*B(2,jj) - A(2,jj)*B(5,jj) - A(3,jj)*B(8,jj)
91     Wi(3) = Wi(3) - A(1,jj)*B(3,jj) - A(2,jj)*B(6,jj) - A(3,jj)*B(9,jj)
92  ENDDO
93END SUBROUTINE D6SDOT
94!======================================================================!
95!                                                                      !
96!======================================================================!
97SUBROUTINE IDNTTY(Neqns,Invp,Iperm)
98  IMPLICIT NONE
99
100  INTEGER :: i
101  INTEGER :: IDBg1
102  INTEGER :: Neqns
103  INTEGER :: Invp(*)
104  INTEGER :: Iperm(*)
105  COMMON /DEBUG / IDBg1
106
107  i = 1
108  DO WHILE ( i<=Neqns )
109     WRITE (6,*) 'invp(', i, ')'
110     READ (5,*) Invp(i)
111     IF ( Invp(i)==0 ) THEN
112        DO i = 1, Neqns
113           Invp(i) = i
114           Iperm(i) = i
115        ENDDO
116        RETURN
117     ELSEIF ( Invp(i)<0 ) THEN
118        READ (11,*) (Invp(i),i=1,Neqns)
119        DO i = 1, Neqns
120           Iperm(Invp(i)) = i
121        ENDDO
122        GOTO 99999
123     ELSE
124        i = i + 1
125     ENDIF
126  ENDDO
127  DO i = 1, Neqns
128     Iperm(Invp(i)) = i
129  ENDDO
130  RETURN
13199999 END SUBROUTINE IDNTTY
132  !======================================================================!
133  !                                                                      !
134  !======================================================================!
135SUBROUTINE NUSOL6(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
136  IMPLICIT NONE
137
138  INTEGER :: i
139  INTEGER :: j
140  INTEGER :: joc
141  INTEGER :: k
142  INTEGER :: ke
143  INTEGER :: ks
144  INTEGER :: Neqns
145  INTEGER :: Nstop
146  INTEGER :: Xlnzr(*)
147  INTEGER :: Colno(*)
148  INTEGER :: Iperm(*)
149  !GP: DEBUG 13May04  wk(3 ---> wk(6, b(3 ---> b(6, diag(6 ---> diag(21,
150  !GP: DEBUG 13May04  zln(9 ---> zln(36, dsln(9 ---> dsln(36
151  !      double precision zln(9,*),diag(6,*),b(3,*),wk(3,*),dsln(9,*)
152  DOUBLE PRECISION :: Zln(36,*)
153  DOUBLE PRECISION :: Diag(21,*)
154  DOUBLE PRECISION :: B(6,*)
155  DOUBLE PRECISION :: Wk(6,*)
156  DOUBLE PRECISION :: Dsln(36,*)
157  ! forward
158  DO i = 1, Neqns
159     Wk(1,i) = B(1,Iperm(i))
160     Wk(2,i) = B(2,Iperm(i))
161     Wk(3,i) = B(3,Iperm(i))
162     Wk(4,i) = B(4,Iperm(i))
163     Wk(5,i) = B(5,Iperm(i))
164     Wk(6,i) = B(6,Iperm(i))
165  ENDDO
166  joc = 1
167  DO i = 1, Neqns
168     ks = Xlnzr(i)
169     ke = Xlnzr(i+1) - 1
170     IF ( ke>=ks ) CALL S6PDOT(Wk(1,i),Wk,Zln,Colno,ks,ke)
171     IF ( i>Nstop ) THEN
172        !        call d6sdot(wk(1,i),wk(1,nstop),dsln(1,joc),i-nstop)
173        CALL DXSDOT(6,Wk(1,i),Wk(1,Nstop),Dsln(1,joc),i-Nstop)
174        joc = joc + i - Nstop
175     ENDIF
176  ENDDO
177  DO i = 1, Neqns
178     Wk(2,i) = Wk(2,i) - Wk(1,i)*Diag(2,i)
179     Wk(3,i) = Wk(3,i) - Wk(1,i)*Diag(4,i) - Wk(2,i)*Diag(5,i)
180     Wk(4,i) = Wk(4,i) - Wk(1,i)*Diag(7,i) - Wk(2,i)*Diag(8,i) - Wk(3,i)*Diag(9,i)
181     Wk(5,i) = Wk(5,i) - Wk(1,i)*Diag(11,i) - Wk(2,i)*Diag(12,i) - Wk(3,i)*Diag(13,i) - Wk(4,i)*Diag(14,i)
182     Wk(6,i) = Wk(6,i) - Wk(1,i)*Diag(16,i) - Wk(2,i)*Diag(17,i)&
183          - Wk(3,i)*Diag(18,i) - Wk(4,i)*Diag(19,i) - Wk(6,i)&
184          *Diag(20,i)
185     Wk(1,i) = Wk(1,i)*Diag(1,i)
186     Wk(2,i) = Wk(2,i)*Diag(3,i)
187     Wk(3,i) = Wk(3,i)*Diag(6,i)
188     Wk(4,i) = Wk(4,i)*Diag(10,i)
189     Wk(5,i) = Wk(5,i)*Diag(15,i)
190     Wk(6,i) = Wk(6,i)*Diag(21,i)
191     Wk(5,i) = Wk(5,i) - Wk(6,i)*Diag(20,i)
192     Wk(4,i) = Wk(4,i) - Wk(6,i)*Diag(19,i) - Wk(5,i)*Diag(14,i)
193     Wk(3,i) = Wk(3,i) - Wk(6,i)*Diag(18,i) - Wk(5,i)*Diag(13,i) - Wk(4,i)*Diag(9,i)
194     Wk(2,i) = Wk(2,i) - Wk(6,i)*Diag(17,i) - Wk(5,i)*Diag(12,i) - Wk(4,i)*Diag(8,i) - Wk(3,i)*Diag(5,i)
195     Wk(1,i) = Wk(1,i) - Wk(6,i)*Diag(16,i) - Wk(5,i)*Diag(11,i)&
196          - Wk(4,i)*Diag(7,i) - Wk(3,i)*Diag(4,i) - Wk(2,i)&
197          *Diag(2,i)
198  ENDDO
199  ! back ward
200  DO i = Neqns, 1, -1
201     IF ( i>=Nstop ) THEN
202        DO j = i - 1, Nstop, -1
203           joc = joc - 1
204           Wk(1,j) = Wk(1,j) - Wk(1,i)*Dsln(1,joc) - Wk(2,i)&
205                *Dsln(2,joc) - Wk(3,i)*Dsln(3,joc) - Wk(4,i)&
206                *Dsln(4,joc) - Wk(5,i)*Dsln(5,joc) - Wk(6,i)&
207                *Dsln(6,joc)
208           Wk(2,j) = Wk(2,j) - Wk(1,i)*Dsln(7,joc) - Wk(2,i)&
209                *Dsln(8,joc) - Wk(3,i)*Dsln(9,joc) - Wk(4,i)&
210                *Dsln(10,joc) - Wk(5,i)*Dsln(11,joc) - Wk(6,i)&
211                *Dsln(12,joc)
212           Wk(3,j) = Wk(3,j) - Wk(1,i)*Dsln(13,joc) - Wk(2,i)&
213                *Dsln(14,joc) - Wk(3,i)*Dsln(15,joc) - Wk(4,i)&
214                *Dsln(16,joc) - Wk(5,i)*Dsln(17,joc) - Wk(6,i)&
215                *Dsln(18,joc)
216           Wk(4,j) = Wk(4,j) - Wk(1,i)*Dsln(19,joc) - Wk(2,i)&
217                *Dsln(20,joc) - Wk(3,i)*Dsln(21,joc) - Wk(4,i)&
218                *Dsln(22,joc) - Wk(5,i)*Dsln(23,joc) - Wk(6,i)&
219                *Dsln(24,joc)
220           Wk(5,j) = Wk(5,j) - Wk(1,i)*Dsln(25,joc) - Wk(2,i)&
221                *Dsln(26,joc) - Wk(3,i)*Dsln(27,joc) - Wk(4,i)&
222                *Dsln(28,joc) - Wk(5,i)*Dsln(29,joc) - Wk(6,i)&
223                *Dsln(30,joc)
224           Wk(6,j) = Wk(6,j) - Wk(1,i)*Dsln(31,joc) - Wk(2,i)&
225                *Dsln(32,joc) - Wk(3,i)*Dsln(33,joc) - Wk(4,i)&
226                *Dsln(34,joc) - Wk(5,i)*Dsln(35,joc) - Wk(6,i)&
227                *Dsln(36,joc)
228        ENDDO
229     ENDIF
230     ks = Xlnzr(i)
231     ke = Xlnzr(i+1) - 1
232     IF ( ke>=ks ) THEN
233        DO k = ks, ke
234           j = Colno(k)
235           Wk(1,j) = Wk(1,j) - Wk(1,i)*Zln(1,joc) - Wk(2,i)&
236                *Zln(2,joc) - Wk(3,i)*Zln(3,joc) - Wk(4,i)&
237                *Zln(4,joc) - Wk(5,i)*Zln(5,joc) - Wk(6,i)&
238                *Zln(6,joc)
239           Wk(2,j) = Wk(2,j) - Wk(1,i)*Zln(7,joc) - Wk(2,i)&
240                *Zln(8,joc) - Wk(3,i)*Zln(9,joc) - Wk(4,i)&
241                *Zln(10,joc) - Wk(5,i)*Zln(11,joc) - Wk(6,i)&
242                *Zln(12,joc)
243           Wk(3,j) = Wk(3,j) - Wk(1,i)*Zln(13,joc) - Wk(2,i)&
244                *Zln(14,joc) - Wk(3,i)*Zln(15,joc) - Wk(4,i)&
245                *Zln(16,joc) - Wk(5,i)*Zln(17,joc) - Wk(6,i)&
246                *Zln(18,joc)
247           Wk(4,j) = Wk(4,j) - Wk(1,i)*Zln(19,joc) - Wk(2,i)&
248                *Zln(20,joc) - Wk(3,i)*Zln(21,joc) - Wk(4,i)&
249                *Zln(22,joc) - Wk(5,i)*Zln(23,joc) - Wk(6,i)&
250                *Zln(24,joc)
251           Wk(5,j) = Wk(5,j) - Wk(1,i)*Zln(25,joc) - Wk(2,i)&
252                *Zln(26,joc) - Wk(3,i)*Zln(27,joc) - Wk(4,i)&
253                *Zln(28,joc) - Wk(5,i)*Zln(29,joc) - Wk(6,i)&
254                *Zln(30,joc)
255           Wk(6,j) = Wk(6,j) - Wk(1,i)*Zln(31,joc) - Wk(2,i)&
256                *Zln(32,joc) - Wk(3,i)*Zln(33,joc) - Wk(4,i)&
257                *Zln(34,joc) - Wk(5,i)*Zln(35,joc) - Wk(6,i)&
258                *Zln(36,joc)
259        ENDDO
260     ENDIF
261  ENDDO
262  ! permutaion
263  DO i = 1, Neqns
264     B(1,Iperm(i)) = Wk(1,i)
265     B(2,Iperm(i)) = Wk(2,i)
266     B(3,Iperm(i)) = Wk(3,i)
267     B(4,Iperm(i)) = Wk(4,i)
268     B(5,Iperm(i)) = Wk(5,i)
269     B(6,Iperm(i)) = Wk(6,i)
270  ENDDO
271END SUBROUTINE NUSOL6
272!======================================================================!
273!                                                                      !
274!======================================================================!
275SUBROUTINE PRT(Ip,N)
276  IMPLICIT NONE
277
278  INTEGER :: i
279  INTEGER :: Ip
280  INTEGER :: N
281  DIMENSION Ip(N)
282  WRITE (6,99001) (Ip(i),i=1,N)
28399001 FORMAT (10(2x,i4))
284END SUBROUTINE PRT
285
286SUBROUTINE VLCPY1(A,N)
287  IMPLICIT NONE
288
289  DOUBLE PRECISION :: A
290  INTEGER :: N
291  DIMENSION A(N)
292  INTEGER :: i
293  INTEGER :: j
294  A(N) = 0
295END SUBROUTINE VLCPY1
296
297!======================================================================!
298!                                                                      !
299!======================================================================!
300SUBROUTINE VERIF0(Neqns,Ndeg,Nttbr,Irow,Jcol,Val,Rhs,X)
301  IMPLICIT NONE
302
303  DOUBLE PRECISION :: err
304  DOUBLE PRECISION :: rel
305  DOUBLE PRECISION :: Rhs
306  DOUBLE PRECISION :: Val
307  DOUBLE PRECISION :: X
308  INTEGER :: i
309  INTEGER :: Irow
310  INTEGER :: j
311  INTEGER :: Jcol
312  INTEGER :: k
313  INTEGER :: l
314  INTEGER :: m
315  INTEGER :: Ndeg
316  INTEGER :: Neqns
317  INTEGER :: Nttbr
318  DIMENSION Irow(*), Jcol(*), Val(Ndeg,Ndeg,*), Rhs(Ndeg,*), X(Ndeg,*)
319  !----------------------------------------------------------------------
320  !
321  !     verify the solution(symmetric matrix)
322  !
323  !----------------------------------------------------------------------
324  rel = 0.0D0
325  DO i = 1, Neqns
326     DO l = 1, Ndeg
327        rel = rel + DABS(Rhs(l,i))
328     ENDDO
329  ENDDO
330  DO k = 1, Nttbr
331     i = Irow(k)
332     j = Jcol(k)
333     DO l = 1, Ndeg
334        DO m = 1, Ndeg
335           Rhs(l,i) = Rhs(l,i) - Val(l,m,k)*X(m,j)
336           IF ( i/=j ) Rhs(l,j) = Rhs(l,j) - Val(m,l,k)*X(m,i)
337        ENDDO
338     ENDDO
339  ENDDO
340  err = 0.0D0
341  DO i = 1, Neqns
342     DO l = 1, Ndeg
343        err = err + DABS(Rhs(l,i))
344     ENDDO
345  ENDDO
346  WRITE (6,99001) err, rel, err/rel
347  !WINDEBUG
348  !      write(16,6000) err,rel,err/rel
34999001 FORMAT (' ***verification***(symmetric)'/&
350       'norm(Ax-b)            =  ',&
351       1pd20.10/'norm(b)               =  ',&
352       1pd20.10/'norm(Ax-b)/norm(b)    =  ',1pd20.10)
353END SUBROUTINE VERIF0
354
355!======================================================================!
356!                                                                      !
357!======================================================================!
358SUBROUTINE V6PROD(Zln,Diag,Zz,N)
359  IMPLICIT NONE
360
361  DOUBLE PRECISION :: Diag
362  DOUBLE PRECISION :: Zln
363  DOUBLE PRECISION :: Zz
364  INTEGER :: i
365  INTEGER :: N
366  DIMENSION Zln(9,N), Diag(6,N), Zz(9,N)
367  DO i = 1, N
368     Zz(4,i) = Zln(4,i) - Zln(1,i)*Diag(2,i)
369     Zz(7,i) = Zln(7,i) - Zln(1,i)*Diag(4,i) - Zz(4,i)*Diag(5,i)
370     Zz(1,i) = Zln(1,i)*Diag(1,i)
371     Zz(4,i) = Zz(4,i)*Diag(3,i)
372     Zz(7,i) = Zz(7,i)*Diag(6,i)
373     Zz(4,i) = Zz(4,i) - Zz(7,i)*Diag(5,i)
374     Zz(1,i) = Zz(1,i) - Zz(4,i)*Diag(2,i) - Zz(7,i)*Diag(4,i)
375     !
376     Zz(5,i) = Zln(5,i) - Zln(2,i)*Diag(2,i)
377     Zz(8,i) = Zln(8,i) - Zln(2,i)*Diag(4,i) - Zz(5,i)*Diag(5,i)
378     Zz(2,i) = Zln(2,i)*Diag(1,i)
379     Zz(5,i) = Zz(5,i)*Diag(3,i)
380     Zz(8,i) = Zz(8,i)*Diag(6,i)
381     Zz(5,i) = Zz(5,i) - Zz(8,i)*Diag(5,i)
382     Zz(2,i) = Zz(2,i) - Zz(5,i)*Diag(2,i) - Zz(8,i)*Diag(4,i)
383     !
384     Zz(6,i) = Zln(6,i) - Zln(3,i)*Diag(2,i)
385     Zz(9,i) = Zln(9,i) - Zln(3,i)*Diag(4,i) - Zz(6,i)*Diag(5,i)
386     Zz(3,i) = Zln(3,i)*Diag(1,i)
387     Zz(6,i) = Zz(6,i)*Diag(3,i)
388     Zz(9,i) = Zz(9,i)*Diag(6,i)
389     Zz(6,i) = Zz(6,i) - Zz(9,i)*Diag(5,i)
390     Zz(3,i) = Zz(3,i) - Zz(6,i)*Diag(2,i) - Zz(9,i)*Diag(4,i)
391  ENDDO
392END SUBROUTINE V6PROD
393