FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_solver_direct.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5 
6 !----------------------------------------------------------------------
8 ! SOLVER=DIRECT
9 !----------------------------------------------------------------------
10 module hecmw_solver_direct
11  use hecmw_util
12  implicit none
13 
14  private
15  public :: hecmw_solve_direct
16 
17  real(kind=kreal), parameter :: rmin = 4.941d-300
18 
19  integer(kind=kint), parameter :: IDBg_ini = 0
20  integer(kind=kint), parameter :: IDBg_sym = 0
21  integer(kind=kint), parameter :: IDBg_num = 0
22 
23  type cholesky_factor
24  integer(kind=kint) :: LEN_colno
25  integer(kind=kint) :: NSTop
26  integer(kind=kint) :: STAge
27  integer(kind=kint) :: NEQns
28  integer(kind=kint) :: NDEg
29  integer(kind=kint) :: NTTbr
30  integer(kind=kint) :: ISYm
31  integer(kind=kint) :: LEN_dsln
32  integer(kind=kint), pointer :: JCOl(:)
33  integer(kind=kint), pointer :: IROw(:)
34  integer(kind=kint), pointer :: IPErm(:)
35  integer(kind=kint), pointer :: INVp(:)
36  integer(kind=kint), pointer :: PARent(:)
37  integer(kind=kint), pointer :: NCH(:)
38  integer(kind=kint), pointer :: XLNzr(:)
39  integer(kind=kint), pointer :: COLno(:)
40  real(kind=kreal), pointer :: diag(:)
41  real(kind=kreal), pointer :: zln(:)
42  real(kind=kreal), pointer :: dsln(:)
43  !*Allocation variables
44  integer(kind=kint) :: ialoc
45  integer(kind=kint) :: raloc
46  end type cholesky_factor
47 
48 contains
49  !----------------------------------------------------------------------
51  ! and is also available for performance tests
52  ! to solve Ax=b for x, nozero pattern and values must be give.
53  ! arrays
54  ! jcol column entries
55  ! irow row entries
56  ! val its value
57  ! v interface array used through matini, staijx, nufctx,
58  ! nusolx
59  !----------------------------------------------------------------------
60  subroutine hecmw_solve_direct(hecMESH,hecMAT,Ifmsg)
61  use hecmw_matrix_ass
62  use hecmw_matrix_dump
63  implicit none
64  !------
65  type (hecmwst_local_mesh), intent(in)::hecmesh
66  type (hecmwst_matrix), intent(inout)::hecmat
67  integer(kind=kint), intent(in):: ifmsg
68  !------
69  type (cholesky_factor), save :: fct
70  !------
71  integer(kind=kint):: i98
72  integer(kind=kint):: i97
73  integer(kind=kint):: timelog
74  integer(kind=kint):: iterlog
75  integer(kind=kint):: ordering
76  integer(kind=kint):: loglevel
77  integer(kind=kint):: ir
78  integer(kind=kint):: i
79  real(kind=kreal):: t1
80  real(kind=kreal):: t2
81  real(kind=kreal):: t3
82  real(kind=kreal):: t4
83  real(kind=kreal):: t5
84 
85  ir = 0
86 
87  timelog = hecmat%IARRAY(22)
88  iterlog = hecmat%IARRAY(21)
89  ordering = hecmat%IARRAY(41)
90  loglevel = max(timelog,iterlog)
91 
92  call hecmw_mat_dump(hecmat,hecmesh)
93 
94  call ptime(t1)
95  t2 = t1
96 
97  !*EHM HECMW June 7 2004
98  i98 = hecmat%IARRAY(98)
99  if ( hecmat%IARRAY(98)==1 ) then
100  !* Interface to symbolic factorization
101  call setij(hecmesh,hecmat,fct)
102 
103  !* Symbolic factorization
104  call matini(fct,ordering,loglevel,ir)
105  hecmat%IARRAY(98) = 0
106 
107  if ( loglevel > 0 ) write (*,*) "[DIRECT]: symbolic fct done"
108  endif
109  call ptime(t2)
110  t3 = t2
111 
112  i97 = hecmat%IARRAY(97)
113  if ( hecmat%IARRAY(97)==1 ) then
114  !* Interface to numeric factorization
115  call nuform(hecmesh,hecmat,fct,ir)
116  call ptime(t3)
117 
118  !* Numeric factorization
119  call nufct0(fct,ir)
120  hecmat%IARRAY(97) = 0
121 
122  if ( loglevel > 0 ) write (*,*) "[DIRECT]: numeric fct done"
123 
124  !*Memory Details
125  if ( loglevel > 1 ) then
126  write (*,*) '*-----------------------------------*'
127  write (*,*) '| Direct Solver Memory Usage |'
128  write (*,*) '*-----------------------------------*'
129  write (*,*) 'INTEGER memory: ', real(fct%IALoc*4)/real(1048576), 'MB'
130  write (*,*) 'REAL*8 memory: ', real(fct%RALoc*8)/real(1048576), 'MB'
131  write (*,*) 'TOTAL memory: ', real((fct%RALoc*2+fct%IALoc)*4)/real(1048576), 'MB'
132  write (*,*) '*-----------------------------------*'
133  endif
134  endif
135  call ptime(t4)
136 
137  !* Finalize
138  !* Errors 1
139  if ( i98/=0 .and. i98/=1 ) then
140  write (ifmsg,*) 'ERROR in symb. fact. flag: Should be 1 or 0'
141  stop 'ERROR in symb. fact. flag: Should be 1 or 0'
142  endif
143  if ( i97/=0 .and. i97/=1 ) then
144  write (ifmsg,*) 'ERROR in numer. fact. flag: Should be 1 or 0'
145  stop 'ERROR in numer. fact. flag: Should be 1 or 0'
146  endif
147  if ( i98==1 .and. i97==0 ) then
148  write (ifmsg,*) 'WARNING: Numeric factorization not performed!'
149  stop 'WARNING: Numeric factorization not performed! Solve will not be performed'
150  endif
151  !* Errors 2
152  if ( ir/=0 ) then
153  write (ifmsg,*) 'ERROR in nufct0. ir = ', ir
154  stop
155  endif
156 
157  !* Solve
158  do i=1,hecmat%NP*hecmesh%n_dof
159  hecmat%X(i) = hecmat%B(i)
160  end do
161  call nusol0(hecmat%X,fct,ir)
162  call ptime(t5)
163  !* Errors 4
164  if ( ir/=0 ) then
165  write (ifmsg,*) 'error in nusol0. irr = ', ir
166  stop
167  endif
168 
169  if ( timelog > 1 ) then
170  write(*,*) 'sym fct time : ',t2-t1
171  write(*,*) 'nuform time : ',t3-t2
172  write(*,*) 'num fct time : ',t4-t3
173  write(*,*) 'solve time : ',t5-t4
174  elseif ( timelog > 0 ) then
175  write(*,*) 'setup time : ',t4-t1
176  write(*,*) 'solve time : ',t5-t4
177  endif
178 
179  call hecmw_mat_dump_solution(hecmat)
180  end subroutine hecmw_solve_direct
181 
182  !======================================================================!
184  !======================================================================!
185  subroutine ptime(Cputim)
186  implicit none
187  !------
188  real(kind=kreal), intent(out):: cputim
189  !------
190  ! cpu time by hour
191  cputim = hecmw_wtime()
192  end subroutine ptime
193 
194  !======================================================================!
196  ! sets NEQns, NDEg, NTTbr, ISYm, JCOl, IROw
197  !======================================================================!
198  subroutine setij(hecMESH,hecMAT,FCT)
199  implicit none
200  !------
201  type (HECMWST_LOCAL_MESH), intent(in)::hecMESH
202  type (HECMWST_MATRIX), intent(in)::hecMAT
203  type (cholesky_factor), intent(inout) :: FCT
204  !------
205  integer(kind=kint):: i
206  integer(kind=kint):: ierr
207  integer(kind=kint):: j
208  integer(kind=kint):: k
209  integer(kind=kint):: kk
210  integer(kind=kint):: ntotal
211  integer(kind=kint):: numnp
212  integer(kind=kint):: ndof
213  integer(kind=kint):: ndof2
214 
215  numnp = hecmat%NP
216  ndof = hecmesh%N_DOF
217  ntotal = numnp*ndof
218 
219  !*NUFACT variables
220  fct%NEQns = numnp
221  fct%NDEg = ndof
222  fct%NTTbr = hecmat%NP + hecmat%NPL
223  !+hecMAT%NPU if unsymmetric
224  fct%ISYm = 0
225 
226  !*Allocations
227  allocate (fct%IROw(fct%NTTbr),stat=ierr)
228  if ( ierr/=0 ) stop "Allocation error: irow"
229  allocate (fct%JCOl(fct%NTTbr),stat=ierr)
230  if ( ierr/=0 ) stop "Allocation error: jcol"
231 
232  kk = 0
233  ndof2 = ndof*ndof
234  do j = 1, numnp
235  !*Diagonal
236  kk = kk + 1
237  fct%IROw(kk) = j
238  fct%JCOl(kk) = j
239  !*Lower
240  do k = hecmat%INDEXL(j-1) + 1, hecmat%INDEXL(j)
241  i = hecmat%ITEML(k)
242  kk = kk + 1
243  fct%IROw(kk) = j
244  fct%JCOl(kk) = i
245  enddo
246  enddo
247  end subroutine setij
248 
249  !======================================================================!
251  ! this routine is used for both symmetric and asymmetric matrices
252  ! and must be called once at the beginning
253  ! (i)
254  ! neqns number of unknowns
255  ! nttbr number of non0s, pattern of non-zero elements are
256  ! given like following.
257  ! nonz(A)={(i,j);i=irow(l),j=jcol(l); 1<= l <= nttbr}
258  ! irow
259  ! jcol to define non-zero pattern
260  ! lenv length of the array v (iv)
261  ! (o)
262  ! iv communication array. v is the original name
263  ! ir return code
264  ! =0 normal
265  ! =-1 non positive index
266  ! =1 too big index
267  ! =10 insufficient storage
268  ! contents of iv
269  ! pointers 1 &zpiv(1) 2 &iperm(1) 3 &invp(1)
270  ! 4 &parent(1)5 &nch(1) 6 &xlnzr(1)
271  ! 7 &colno(1) 8 &diag(1) 9 &zln(1)
272  ! 10 &dsln(1)
273  !
274  ! scalars 21 len(colno) 22 nstop 23 stage
275  ! 24 neqns 25 len(iv) 26 len(dsln)
276  ! 27 total
277  ! stage 10 after initialization
278  ! 20 building up matrix
279  ! 30 after LU decomposition
280  ! 40 after solving
281  !
282  ! sets LEN_colno, NSTop, LEN_dsln, IPErm, INVp, PARent, NCH, XLNzr, COLno
283  ! IALoc, RALoc, STAge
284  !======================================================================!
285  subroutine matini(FCT,ordering,loglevel,Ir)
286  use hecmw_ordering
287  implicit none
288  !------
289  type (cholesky_factor), intent(inout) :: FCT
290  integer(kind=kint), intent(in):: ordering
291  integer(kind=kint), intent(in):: loglevel
292  integer(kind=kint), intent(out):: Ir
293  !------
294  !*Work arrays
295  integer(kind=kint), allocatable :: zpiv(:)
296  integer(kind=kint), allocatable :: jcpt(:)
297  integer(kind=kint), allocatable :: jcolno(:)
298  integer(kind=kint), allocatable :: ia(:)
299  integer(kind=kint), allocatable :: ja(:)
300  integer(kind=kint), allocatable :: quarent(:)
301  integer(kind=kint), allocatable :: btree(:)
302  integer(kind=kint), allocatable :: xleaf(:)
303  integer(kind=kint), allocatable :: leaf(:)
304  integer(kind=kint):: ir1
305  integer(kind=kint):: irr
306  integer(kind=kint):: izz
307  integer(kind=kint):: izz0
308  integer(kind=kint):: lncol
309  integer(kind=kint):: neqnsz
310  integer(kind=kint):: ierror
311 
312  izz0 = 0
313 
314  ir = 0
315 
316  !*Initialize allocation measure variables
317  fct%IALoc = 0
318  fct%RALoc = 0
319  !
320  ! set z pivot
321  !
322  allocate (zpiv(fct%NEQns),stat=ierror)
323  if ( ierror/=0 ) stop "ALLOCATION ERROR, zpiv: SUB. matini"
324  call zpivot(fct%NEQns,neqnsz,fct%NTTbr,fct%JCOl,fct%IROw,zpiv,ir1)
325  if ( ir1/=0 ) then
326  ir = ir1
327  return
328  endif
329 
330  !
331  ! build jcpt,jcolno
332  !
333  allocate (jcpt(2*fct%NTTbr),stat=ierror)
334  if ( ierror/=0 ) stop "ALLOCATION ERROR, jcpt: SUB. matini"
335  allocate (jcolno(2*fct%NTTbr),stat=ierror)
336  if ( ierror/=0 ) stop "ALLOCATION ERROR, jcolno: SUB. matini"
337  call stsmat(fct%NEQns,fct%NTTbr,fct%IROw,fct%JCOl,jcpt,jcolno)
338  !
339  ! build ia,ja
340  !
341  allocate (ia(fct%NEQns+1),stat=ierror)
342  if ( ierror/=0 ) stop "ALLOCATION ERROR, ia: SUB. matini"
343  allocate (ja(2*fct%NTTbr),stat=ierror)
344  if ( ierror/=0 ) stop "ALLOCATION ERROR, ja: SUB. matini"
345  call stiaja(fct%NEQns,ia,ja,jcpt,jcolno)
346 
347  !*Deallocation of work array
348  deallocate (jcpt)
349  deallocate (jcolno)
350  !
351  ! get permutation vector iperm,invp
352  !
353  allocate (fct%IPErm(fct%NEQns),stat=ierror)
354  if ( ierror/=0 ) stop "ALLOCATION ERROR, iperm: SUB. matini"
355  allocate (fct%INVp(fct%NEQns),stat=ierror)
356  if ( ierror/=0 ) stop "ALLOCATION ERROR, invp: SUB. matini"
357  allocate (quarent(fct%NEQns+1),stat=ierror)
358  if ( ierror/=0 ) stop "ALLOCATION ERROR, quarent: SUB. matini"
359  call hecmw_ordering_gen(neqnsz,fct%NTTbr,ia,ja,fct%IPErm,fct%INVp,ordering,loglevel)
360  do
361  ! build up the parent vector
362  ! parent vector will be saved in Quarent for a while
363  call genpaq(ia,ja,fct%INVp,fct%IPErm,quarent,fct%NEQns)
364  !
365  ! build up the binary tree
366  !
367  allocate (btree(2*(fct%NEQns+1)),stat=ierror)
368  if ( ierror/=0 ) stop "ALLOCATION ERROR, btree: SUB. matini"
369  call genbtq(fct%INVp,quarent,btree,zpiv,izz,fct%NEQns)
370  !
371  ! rotate the binary tree to avoid a zero pivot
372  !
373  if ( izz==0 ) then
374  !
375  ! post ordering
376  !
377  allocate (fct%PARent(fct%NEQns),stat=ierror)
378  if ( ierror/=0 ) stop "ALLOCATION ERROR, parent: SUB. matini.f"
379  allocate (fct%NCH(fct%NEQns+1),stat=ierror)
380  if ( ierror/=0 ) stop "ALLOCATION ERROR, nch: SUB. matini.f"
381  call posord(fct%PARent,btree,fct%INVp,fct%IPErm,fct%NCH,fct%NEQns,quarent)
382  !
383  ! generate skeleton graph
384  !
385  allocate (xleaf(fct%NEQns+1),stat=ierror)
386  if ( ierror/=0 ) stop "ALLOCATION ERROR, xleaf: SUB. matini.f"
387  allocate (leaf(fct%NTTbr),stat=ierror)
388  if ( ierror/=0 ) stop "ALLOCATION ERROR, leaf: SUB. matini.f"
389  call gnleaf(ia,ja,fct%INVp,fct%IPErm,fct%NCH,xleaf,leaf,fct%NEQns)
390  call forpar(fct%NEQns,fct%PARent,fct%NCH,fct%NSTop)
391  !*Deallocation of work arrays
392  deallocate (ia)
393  deallocate (ja)
394  deallocate (quarent)
395  deallocate (zpiv)
396  !
397  ! build up xlnzr,colno (this is the symbolic fct.)
398  !
399  allocate (fct%XLNzr(fct%NEQns+1),stat=ierror)
400  if ( ierror/=0 ) stop "ALLOCATION ERROR, xlnzr: SUB. matini.f"
401  call pre_gnclno(fct%PARent,xleaf,leaf,fct%XLNzr,fct%NEQns,fct%NSTop,lncol,ir1)
402  allocate (fct%COLno(lncol),stat=ierror)
403  if ( ierror/=0 ) stop "ALLOCATION ERROR, colno: SUB. matini.f"
404  call gnclno(fct%PARent,xleaf,leaf,fct%XLNzr,fct%COLno,fct%NEQns,fct%NSTop,lncol,ir1)
405  !*Deallocate work arrays
406  deallocate (xleaf)
407  deallocate (leaf)
408  deallocate (btree)
409  fct%LEN_dsln = (fct%NEQns-fct%NSTop+1)*(fct%NEQns-fct%NSTop)/2
410 
411  !Scalar assignments
412  fct%LEN_colno = lncol
413 
414  fct%STAge = 10
415  fct%IALoc = 5*fct%NEQns + lncol + 1
416  exit
417  else
418  if ( izz0==0 ) izz0 = izz
419  if ( izz0/=izz ) then
420  call bringu(zpiv,fct%IPErm,fct%INVp,quarent,izz,fct%NEQns,irr)
421  else
422  call rotate(ia,ja,fct%INVp,fct%IPErm,quarent,btree,izz,fct%NEQns,irr)
423  endif
424  endif
425  enddo
426  end subroutine matini
427 
428  !======================================================================!
430  !======================================================================!
431  subroutine zpivot(Neqns,Neqnsz,Nttbr,Jcol,Irow,Zpiv,Ir)
432  implicit none
433  !------
434  integer(kind=kint), intent(in):: Neqns
435  integer(kind=kint), intent(in):: Nttbr
436  integer(kind=kint), intent(in):: Jcol(:)
437  integer(kind=kint), intent(in):: Irow(:)
438  integer(kind=kint), intent(out):: Ir
439  integer(kind=kint), intent(out):: Neqnsz
440  integer(kind=kint), intent(out):: Zpiv(:)
441  !------
442  integer(kind=kint):: i
443  integer(kind=kint):: j
444  integer(kind=kint):: l
445 
446  ir = 0
447  zpiv(1:neqns) = 1
448 
449  loop1: do
450  do l = 1, nttbr
451  i = irow(l)
452  j = jcol(l)
453  if ( i<=0 .or. j<=0 ) then
454  ir = -1
455  exit loop1
456  elseif ( i>neqns .or. j>neqns ) then
457  ir = 1
458  exit loop1
459  endif
460  if ( i==j ) zpiv(i) = 0
461  enddo
462  do i = neqns, 1, -1
463  if ( zpiv(i)==0 ) then
464  neqnsz = i
465  exit
466  endif
467  enddo
468  exit
469  enddo loop1
470  if ( idbg_ini/=0 ) write (6,"(20I3)") (zpiv(i),i=1,neqns)
471  end subroutine zpivot
472 
473  !======================================================================!
475  !======================================================================!
476  subroutine stsmat(Neqns,Nttbr,Irow,Jcol,Jcpt,Jcolno)
477  implicit none
478  !------
479  integer(kind=kint), intent(in):: Neqns
480  integer(kind=kint), intent(in):: Nttbr
481  integer(kind=kint), intent(in):: Irow(:)
482  integer(kind=kint), intent(in):: Jcol(:)
483  integer(kind=kint), intent(out):: Jcpt(:)
484  integer(kind=kint), intent(out):: Jcolno(:)
485  !------
486  integer(kind=kint):: i
487  integer(kind=kint):: j
488  integer(kind=kint):: joc
489  integer(kind=kint):: k
490  integer(kind=kint):: l
491  integer(kind=kint):: locr
492  logical:: found
493 
494  jcpt(1:2*nttbr) = 0
495  jcolno(1:2*nttbr) = 0
496  do i = 1, neqns
497  jcpt(i) = i + neqns
498  jcolno(i+neqns) = i
499  enddo
500 
501  k = 2*neqns
502 
503  loop1: do l = 1, nttbr
504  i = irow(l)
505  j = jcol(l)
506  if ( i/=j ) then
507  joc = jcpt(i)
508  locr = i
509  found = .false.
510  do while ( joc/=0 )
511  if ( jcolno(joc)==j ) cycle loop1
512  if ( jcolno(joc)>j ) then
513  k = k + 1
514  jcpt(locr) = k
515  jcpt(k) = joc
516  jcolno(k) = j
517  found = .true.
518  exit
519  endif
520  locr = joc
521  joc = jcpt(joc)
522  enddo
523  if (.not. found) then
524  k = k + 1
525  jcpt(locr) = k
526  jcolno(k) = j
527  endif
528 
529  joc = jcpt(j)
530  locr = j
531  do while ( joc/=0 )
532  if ( jcolno(joc)==i ) cycle loop1
533  if ( jcolno(joc)>i ) then
534  k = k + 1
535  jcpt(locr) = k
536  jcpt(k) = joc
537  jcolno(k) = i
538  cycle loop1
539  endif
540  locr = joc
541  joc = jcpt(joc)
542  enddo
543  k = k + 1
544  jcpt(locr) = k
545  jcolno(k) = i
546  endif
547  enddo loop1
548  if ( idbg_ini/=0 ) then
549  write (6,*) 'jcolno'
550  write (6,"(10I7)") (jcolno(i),i=1,k)
551  write (6,*) 'jcpt'
552  write (6,"(10I7)") (jcpt(i),i=1,k)
553  endif
554  end subroutine stsmat
555 
556  !======================================================================!
558  ! (asymmetric version)
559  !======================================================================!
560  subroutine stiaja(Neqns,Ia,Ja,Jcpt,Jcolno)
561  implicit none
562  !------
563  integer(kind=kint), intent(in):: Neqns
564  integer(kind=kint), intent(in):: Jcpt(:)
565  integer(kind=kint), intent(in):: Jcolno(:)
566  integer(kind=kint), intent(out):: Ia(:)
567  integer(kind=kint), intent(out):: Ja(:)
568  !------
569  integer(kind=kint):: i
570  integer(kind=kint):: ii
571  integer(kind=kint):: joc
572  integer(kind=kint):: k
573  integer(kind=kint):: l
574 
575  ia(1) = 1
576  l = 0
577  do k = 1, neqns
578  joc = jcpt(k)
579  do while ( joc/=0 )
580  ii = jcolno(joc)
581  if ( ii/=k ) then
582  l = l + 1
583  ja(l) = ii
584  endif
585  joc = jcpt(joc)
586  enddo
587  ia(k+1) = l + 1
588  enddo
589  if ( idbg_ini/=0 ) then
590  write (6,*) 'ia '
591  write (6,"(10I7)") (ia(i),i=1,neqns)
592  write (6,*) 'ja '
593  write (6,"(10I7)") (ja(i),i=1,ia(neqns+1))
594  endif
595  end subroutine stiaja
596 
597  !======================================================================!
599  !======================================================================!
600  subroutine genpaq(Xadj,Adjncy,Invp,Iperm,Parent,Neqns)
601  implicit none
602  !------
603  integer(kind=kint), intent(in):: Neqns
604  integer(kind=kint), intent(in):: Xadj(:)
605  integer(kind=kint), intent(in):: Adjncy(:)
606  integer(kind=kint), intent(in):: Invp(:)
607  integer(kind=kint), intent(in):: Iperm(:)
608  integer(kind=kint), intent(out):: Parent(:)
609  !------
610  integer(kind=kint), allocatable:: Ancstr(:)
611  integer(kind=kint):: i
612  integer(kind=kint):: ip
613  integer(kind=kint):: it
614  integer(kind=kint):: k
615  integer(kind=kint):: l
616  integer(kind=kint):: ierror
617 
618  allocate (ancstr(neqns+1),stat=ierror)
619  if ( ierror/=0 ) stop "ALLOCATION ERROR, ancstr: SUB. genpaq"
620 
621  do i = 1, neqns
622  parent(i) = 0
623  ancstr(i) = 0
624  ip = iperm(i)
625  loop1: do k = xadj(ip), xadj(ip+1) - 1
626  l = invp(adjncy(k))
627  if ( l<i ) then
628  do while ( ancstr(l)/=0 )
629  if ( ancstr(l)==i ) cycle loop1
630  it = ancstr(l)
631  ancstr(l) = i
632  l = it
633  enddo
634  ancstr(l) = i
635  parent(l) = i
636  endif
637  enddo loop1
638  enddo
639  do i = 1, neqns
640  if ( parent(i)==0 ) parent(i) = neqns + 1
641  enddo
642  parent(neqns+1) = 0
643 
644  if ( idbg_sym/=0 ) then
645  write (6,"(' parent')")
646  write (6,"(2I6)") (i,parent(i),i=1,neqns)
647  endif
648 
649  deallocate(ancstr)
650  end subroutine genpaq
651 
652  !======================================================================!
654  !======================================================================!
655  subroutine genbtq(Invp,Parent,Btree,Zpiv,Izz,Neqns)
656  implicit none
657  !------
658  integer(kind=kint), intent(in):: Neqns
659  integer(kind=kint), intent(in):: Parent(:)
660  integer(kind=kint), intent(in):: Invp(:)
661  integer(kind=kint), intent(in):: Zpiv(:)
662  integer(kind=kint), intent(out):: Btree(2,*)
663  !------
664  integer(kind=kint):: i
665  integer(kind=kint):: ib
666  integer(kind=kint):: inext
667  integer(kind=kint):: ip
668  integer(kind=kint):: Izz
669 
670  btree(1,1:neqns + 1) = 0
671  btree(2,1:neqns + 1) = 0
672  do i = 1, neqns + 1
673  ip = parent(i)
674  if ( ip>0 ) then
675  ib = btree(1,ip)
676  if ( ib==0 ) then
677  btree(1,ip) = i
678  else
679  do
680  inext = btree(2,ib)
681  if ( inext==0 ) then
682  btree(2,ib) = i
683  else
684  ib = inext
685  cycle
686  endif
687  exit
688  enddo
689  endif
690  endif
691  enddo
692  !
693  ! find zeropivot
694  !
695  izz = 0
696  do i = 1, neqns
697  if ( zpiv(i)/=0 ) then
698  if ( btree(1,invp(i))==0 ) then
699  izz = i
700  exit
701  endif
702  endif
703  enddo
704 
705  if ( idbg_sym/=0 ) then
706  write (6,"(' binary tree')")
707  write (6,"(i6,'(',2I6,')')") (i,btree(1,i),btree(2,i),i=1,neqns)
708  write (6,"(' the first zero pivot is ',i4)") izz
709  endif
710  end subroutine genbtq
711 
712  !======================================================================!
714  !======================================================================!
715  subroutine posord(Parent,Btree,Invp,Iperm,Nch,Neqns,Qarent)
716  implicit none
717  !------
718  integer(kind=kint), intent(in):: Neqns
719  integer(kind=kint), intent(in):: Btree(2,*)
720  integer(kind=kint), intent(in):: Qarent(:)
721  integer(kind=kint), intent(out):: Parent(:)
722  integer(kind=kint), intent(inout):: Invp(:)
723  integer(kind=kint), intent(out):: Iperm(:)
724  integer(kind=kint), intent(out):: Nch(:)
725  !------
726  integer(kind=kint), allocatable:: Pordr(:)
727  integer(kind=kint), allocatable:: Iw(:)
728  integer(kind=kint), allocatable:: Mch(:)
729  integer(kind=kint):: i
730  integer(kind=kint):: ii
731  integer(kind=kint):: invpos
732  integer(kind=kint):: ipinv
733  integer(kind=kint):: joc
734  integer(kind=kint):: l
735  integer(kind=kint):: locc
736  integer(kind=kint):: locp
737  integer(kind=kint):: ierror
738 
739  allocate (pordr(neqns+1),stat=ierror)
740  if ( ierror/=0 ) stop "ALLOCATION ERROR, pordr: SUB. posord"
741  allocate (iw(neqns+1),stat=ierror)
742  if ( ierror/=0 ) stop "ALLOCATION ERROR, iw: SUB. posord"
743  allocate (mch(0:neqns+1),stat=ierror)
744  if ( ierror/=0 ) stop "ALLOCATION ERROR, mch: SUB. posord"
745 
746  mch(1:neqns) = 0
747  pordr(1:neqns) = 0
748  l = 1
749  locc = neqns + 1
750  do
751  joc = locc
752  locc = btree(1,joc)
753  if ( locc==0 ) then
754  locp = qarent(joc)
755  mch(locp) = mch(locp) + 1
756  do
757  pordr(joc) = l
758  if ( l>=neqns ) then
759  do i = 1, neqns
760  ipinv = pordr(invp(i))
761  invp(i) = ipinv
762  iperm(ipinv) = i
763  iw(pordr(i)) = i
764  enddo
765  do i = 1, neqns
766  invpos = iw(i)
767  nch(i) = mch(invpos)
768  ii = qarent(invpos)
769  if ( ii>0 .and. ii<=neqns ) then
770  parent(i) = pordr(ii)
771  else
772  parent(i) = qarent(invpos)
773  endif
774  enddo
775  if ( idbg_sym/=0 ) then
776  write (6,"(' post order')")
777  write (6,"(10I6)") (pordr(i),i=1,neqns)
778  write (6,"(/' invp ')")
779  write (6,"(/' parent')")
780  write (6,"(10I6)") (parent(i),i=1,neqns)
781  write (6,"(10I6)") (invp(i),i=1,neqns)
782  write (6,"(/' iperm ')")
783  write (6,"(10I6)") (iperm(i),i=1,neqns)
784  write (6,"(' nch')")
785  write (6,"(10I6)") (nch(i),i=1,neqns)
786  endif
787  return
788  else
789  l = l + 1
790  locc = btree(2,joc)
791  if ( locc/=0 ) exit
792  joc = qarent(joc)
793  locp = qarent(joc)
794  mch(locp) = mch(locp) + mch(joc) + 1
795  endif
796  enddo
797  endif
798  enddo
799 
800  deallocate (pordr)
801  deallocate (iw)
802  deallocate (mch)
803  end subroutine posord
804 
805  !======================================================================!
807  !======================================================================!
808  subroutine gnleaf(Xadj,Adjncy,Invp,Iperm,Nch,Xleaf,Leaf,Neqns)
809  implicit none
810  !------
811  integer(kind=kint), intent(in):: Neqns
812  integer(kind=kint), intent(in):: Xadj(:)
813  integer(kind=kint), intent(in):: Adjncy(:)
814  integer(kind=kint), intent(in):: Nch(:)
815  integer(kind=kint), intent(in):: Invp(:)
816  integer(kind=kint), intent(in):: Iperm(:)
817  integer(kind=kint), intent(out):: Xleaf(:)
818  integer(kind=kint), intent(out):: Leaf(:)
819  !------
820  integer(kind=kint), allocatable:: Adjncp(:)
821  integer(kind=kint):: Lnleaf
822  integer(kind=kint):: i
823  integer(kind=kint):: ik
824  integer(kind=kint):: ip
825  integer(kind=kint):: iq
826  integer(kind=kint):: istart
827  integer(kind=kint):: k
828  integer(kind=kint):: l
829  integer(kind=kint):: lc
830  integer(kind=kint):: lc1
831  integer(kind=kint):: m
832  integer(kind=kint):: ierror
833 
834  allocate (adjncp(neqns+1),stat=ierror)
835  if ( ierror/=0 ) stop "ALLOCATION ERROR, adjncp: SUB. gnleaf"
836 
837  l = 1
838  ik = 0
839  istart = 0
840  do i = 1, neqns
841  xleaf(i) = l
842  ip = iperm(i)
843  do k = xadj(ip), xadj(ip+1) - 1
844  iq = invp(adjncy(k))
845  if ( iq<i ) then
846  ik = ik + 1
847  adjncp(ik) = iq
848  endif
849  enddo
850  m = ik - istart
851  if ( m/=0 ) then
852  call qqsort(adjncp(istart+1:),m)
853  lc1 = adjncp(istart+1)
854  if ( lc1<i ) then
855  leaf(l) = lc1
856  l = l + 1
857  do k = istart + 2, ik
858  lc = adjncp(k)
859  if ( lc1<lc-nch(lc) ) then
860  leaf(l) = lc
861  l = l + 1
862  endif
863 
864  lc1 = lc
865  enddo
866  ik = 1
867  istart = ik
868  endif
869  endif
870  enddo
871  xleaf(neqns+1) = l
872  lnleaf = l - 1
873 
874  if ( idbg_sym/=0 ) then
875  write (6,"(' xleaf')")
876  write (6,"(10I6)") (xleaf(i),i=1,neqns+1)
877  write (6,"(' leaf (len = ',i6,')')") lnleaf
878  write (6,"(10I6)") (leaf(i),i=1,lnleaf)
879  endif
880 
881  deallocate (adjncp)
882  return
883  end subroutine gnleaf
884 
885  !======================================================================!
887  ! sort in increasing order up to i
888  ! iw array
889  ! ik number of input/output
890  ! i deal with numbers less than this numberi
891  !======================================================================!
892  subroutine qqsort(Iw,Ik)
893  implicit none
894  !------
895  integer(kind=kint), intent(in):: Ik
896  integer(kind=kint), intent(inout):: Iw(:)
897  !------
898  integer(kind=kint):: itemp
899  integer(kind=kint):: l
900  integer(kind=kint):: m
901 
902  if ( ik<=1 ) return
903  do l = 1, ik - 1
904  do m = l + 1, ik
905  if ( iw(l)>=iw(m) ) then
906  itemp = iw(l)
907  iw(l) = iw(m)
908  iw(m) = itemp
909  endif
910  enddo
911  enddo
912  end subroutine qqsort
913 
914  !======================================================================!
916  !======================================================================!
917  subroutine forpar(Neqns,Parent,Nch,Nstop)
918  implicit none
919  !------
920  integer(kind=kint), intent(in):: Neqns
921  integer(kind=kint), intent(in):: Parent(:)
922  integer(kind=kint), intent(out):: Nstop
923  integer(kind=kint), intent(out):: Nch(:)
924  !------
925  integer(kind=kint):: i
926  integer(kind=kint):: idens
927  integer(kind=kint):: ii
928 
929  nch(1:neqns) = 0
930  nch(neqns+1) = 0
931  do i = 1, neqns
932  ii = parent(i)
933  nch(ii) = nch(ii) + 1
934  enddo
935  do i = neqns, 1, -1
936  if ( nch(i)/=1 ) exit
937  enddo
938 
939  idens = 0
940  if ( idens==1 ) then
941  nstop = i
942  else
943  nstop = neqns + 1
944  endif
945  end subroutine forpar
946 
947  !======================================================================!
949  !======================================================================!
950  subroutine pre_gnclno(Parent,Xleaf,Leaf,Xlnzr,Neqns,Nstop,Lncol,Ir)
951  implicit none
952  !------
953  integer(kind=kint), intent(in):: Neqns
954  integer(kind=kint), intent(in):: Nstop
955  integer(kind=kint), intent(in):: Parent(:)
956  integer(kind=kint), intent(in):: Xleaf(:)
957  integer(kind=kint), intent(in):: Leaf(:)
958  integer(kind=kint), intent(out):: Lncol
959  integer(kind=kint), intent(out):: Xlnzr(:)
960  !------
961  integer(kind=kint):: i
962  integer(kind=kint):: Ir
963  integer(kind=kint):: j
964  integer(kind=kint):: k
965  integer(kind=kint):: ke
966  integer(kind=kint):: ks
967  integer(kind=kint):: l
968  integer(kind=kint):: nc
969  integer(kind=kint):: nxleaf
970 
971  nc = 0
972  ir = 0
973  l = 1
974  loop1: do i = 1, neqns
975  xlnzr(i) = l
976  ks = xleaf(i)
977  ke = xleaf(i+1) - 1
978  if ( ke>=ks ) then
979  nxleaf = leaf(ks)
980  do k = ks, ke - 1
981  j = nxleaf
982  nxleaf = leaf(k+1)
983  do while ( j<nxleaf )
984  if ( j>=nstop ) cycle loop1
985  l = l + 1
986  j = parent(j)
987  enddo
988  enddo
989  j = leaf(ke)
990  do while ( j<nstop )
991  if ( j>=i .or. j==0 ) exit
992  l = l + 1
993  j = parent(j)
994  enddo
995  endif
996  enddo loop1
997  xlnzr(neqns+1) = l
998  lncol = l - 1
999  end subroutine pre_gnclno
1000 
1001  !======================================================================!
1003  !======================================================================!
1004  subroutine gnclno(Parent,Xleaf,Leaf,Xlnzr,Colno,Neqns,Nstop,Lncol,Ir)
1005  implicit none
1006  !------
1007  integer(kind=kint), intent(in):: Neqns
1008  integer(kind=kint), intent(in):: Nstop
1009  integer(kind=kint), intent(in):: Parent(:)
1010  integer(kind=kint), intent(in):: Xleaf(:)
1011  integer(kind=kint), intent(in):: Leaf(:)
1012  integer(kind=kint), intent(out):: Ir
1013  integer(kind=kint), intent(out):: Lncol
1014  integer(kind=kint), intent(out):: Xlnzr(:)
1015  integer(kind=kint), intent(out):: Colno(:)
1016  !------
1017  integer(kind=kint):: i
1018  integer(kind=kint):: j
1019  integer(kind=kint):: k
1020  integer(kind=kint):: ke
1021  integer(kind=kint):: ks
1022  integer(kind=kint):: l
1023  integer(kind=kint):: nc
1024  integer(kind=kint):: nxleaf
1025 
1026  nc = 0
1027  ir = 0
1028  l = 1
1029  loop1: do i = 1, neqns
1030  xlnzr(i) = l
1031  ks = xleaf(i)
1032  ke = xleaf(i+1) - 1
1033  if ( ke>=ks ) then
1034  nxleaf = leaf(ks)
1035  do k = ks, ke - 1
1036  j = nxleaf
1037  nxleaf = leaf(k+1)
1038  do while ( j<nxleaf )
1039  if ( j>=nstop ) cycle loop1
1040  colno(l) = j
1041  l = l + 1
1042  j = parent(j)
1043  enddo
1044  enddo
1045  j = leaf(ke)
1046  do while ( j<nstop )
1047  if ( j>=i .or. j==0 ) exit
1048  colno(l) = j
1049  l = l + 1
1050  j = parent(j)
1051  enddo
1052  endif
1053  enddo loop1
1054  xlnzr(neqns+1) = l
1055  lncol = l - 1
1056 
1057  if ( idbg_sym/=0 ) then
1058  write (6,"(' xlnzr')")
1059  write (6,"(' colno (lncol =',i10,')')") lncol
1060  do k = 1, neqns
1061  write (6,"(/' row = ',i6)") k
1062  write (6,"(10I4)") (colno(i),i=xlnzr(k),xlnzr(k+1)-1)
1063  enddo
1064  endif
1065  end subroutine gnclno
1066 
1067  !======================================================================!
1069  ! irr = 0 complete
1070  ! = 1 impossible
1071  !======================================================================!
1072  subroutine bringu(Zpiv,Iperm,Invp,Parent,Izz,Neqns,Irr)
1073  implicit none
1074  !------
1075  integer(kind=kint), intent(in):: Izz
1076  integer(kind=kint), intent(in):: Neqns
1077  integer(kind=kint), intent(in):: Zpiv(:)
1078  integer(kind=kint), intent(in):: Parent(:)
1079  integer(kind=kint), intent(out):: Irr
1080  integer(kind=kint), intent(inout):: Iperm(:)
1081  integer(kind=kint), intent(inout):: Invp(:)
1082  !------
1083  integer(kind=kint):: i
1084  integer(kind=kint):: ib
1085  integer(kind=kint):: ib0
1086  integer(kind=kint):: ibp
1087  integer(kind=kint):: idbg
1088  integer(kind=kint):: izzp
1089 
1090  idbg = 0
1091  irr = 0
1092  ib0 = invp(izz)
1093  ib = ib0
1094  do while ( ib>0 )
1095  ibp = parent(ib)
1096  izzp = iperm(ibp)
1097  if ( zpiv(izzp)==0 ) then
1098  invp(izz) = ibp
1099  invp(izzp) = ib0
1100  iperm(ibp) = izz
1101  iperm(ib0) = izzp
1102  if ( idbg/=0 ) then
1103  do i = 1, neqns
1104  if ( invp(iperm(i))/=i .or. iperm(invp(i))/=i) then
1105  write (6,*) 'permutation error'
1106  stop
1107  endif
1108  enddo
1109  return
1110  endif
1111  return
1112  else
1113  ib = ibp
1114  endif
1115  enddo
1116  irr = 1
1117  end subroutine bringu
1118 
1119  !======================================================================!
1121  ! irr return code irr=0 node izz is not a bottom node
1122  ! irr=1 is a bottom node then rotation is
1123  ! performed
1124  !======================================================================!
1125  subroutine rotate(Xadj,Adjncy,Invp,Iperm,Parent,Btree,Izz,Neqns,Irr)
1126  implicit none
1127  !------
1128  integer(kind=kint), intent(in):: Izz
1129  integer(kind=kint), intent(in):: Neqns
1130  integer(kind=kint), intent(in):: Xadj(:)
1131  integer(kind=kint), intent(in):: Adjncy(:)
1132  integer(kind=kint), intent(in):: Parent(:)
1133  integer(kind=kint), intent(in):: Btree(2,*)
1134  integer(kind=kint), intent(out):: Irr
1135  integer(kind=kint), intent(inout):: Invp(:)
1136  integer(kind=kint), intent(inout):: Iperm(:)
1137  !------
1138  integer(kind=kint), allocatable:: Anc(:)
1139  integer(kind=kint), allocatable:: Adjt(:)
1140  integer(kind=kint):: i
1141  integer(kind=kint):: iy
1142  integer(kind=kint):: izzz
1143  integer(kind=kint):: joc
1144  integer(kind=kint):: k
1145  integer(kind=kint):: l
1146  integer(kind=kint):: ll
1147  integer(kind=kint):: locc
1148  integer(kind=kint):: nanc
1149  integer(kind=kint):: ierror
1150 
1151  allocate (anc(neqns+1),stat=ierror)
1152  if ( ierror/=0 ) stop "ALLOCATION ERROR, anc: SUB. rotate"
1153  allocate (adjt(neqns+1),stat=ierror)
1154  if ( ierror/=0 ) stop "ALLOCATION ERROR, adjt: SUB. rotate"
1155 
1156  if ( izz==0 ) then
1157  irr = 0
1158  return
1159  endif
1160  izzz = invp(izz)
1161 
1162  if ( btree(1,izzz)/=0 ) irr = 0
1163  irr = 1
1164  !
1165  ! ancestors of izzz
1166  !
1167  nanc = 0
1168  joc = izzz
1169  do
1170  nanc = nanc + 1
1171  anc(nanc) = joc
1172  joc = parent(joc)
1173  if ( joc==0 ) then
1174  !
1175  ! to find the eligible node from ancestors of izz
1176  !
1177  l = 1
1178  exit
1179  endif
1180  enddo
1181 
1182  loop1: do
1183  adjt(1:neqns) = 0
1184  locc = anc(l)
1185  do
1186  joc = locc
1187  locc = btree(1,joc)
1188  if ( locc==0 ) then
1189  do
1190  adjt(invp(adjncy(xadj(iperm(joc)):xadj(iperm(joc)+1) - 1))) = 1
1191  if ( joc>=anc(l) ) then
1192  do ll = l + 1, nanc
1193  if ( adjt(anc(ll))==0 ) then
1194  l = l + 1
1195  cycle loop1
1196  endif
1197  enddo
1198  if ( l==1 ) then
1199  !
1200  ! izz can be numbered last
1201  !
1202  k = 0
1203  do i = 1, neqns
1204  if ( i/=izzz ) then
1205  k = k + 1
1206  invp(iperm(i)) = k
1207  endif
1208  enddo
1209  invp(iperm(izzz)) = neqns
1210  else
1211  !
1212  ! anc(l-1) is the eligible node
1213  !
1214  ! (1) number the node not in Ancestor(iy)
1215  iy = anc(l-1)
1216  adjt(1:neqns) = 0
1217  adjt(anc(l:nanc)) = 1
1218  k = 0
1219  do ll = 1, neqns
1220  if ( adjt(ll)==0 ) then
1221  k = k + 1
1222  invp(iperm(ll)) = k
1223  endif
1224  enddo
1225  ! (2) followed by nodes in Ancestor(iy)-Adj(T(iy))
1226  adjt(1:neqns) = 0
1227  locc = iy
1228  loop2: do
1229  joc = locc
1230  locc = btree(1,joc)
1231  if ( locc==0 ) then
1232  do
1233  adjt(invp(adjncy(xadj(iperm(joc)):xadj(iperm(joc)+1)-1))) = 1
1234  if ( joc>=iy ) then
1235  do ll = l, nanc
1236  if ( adjt(anc(ll))==0 ) then
1237  k = k + 1
1238  invp(iperm(anc(ll))) = k
1239  endif
1240  enddo
1241  ! (3) and finally number the node in Adj(t(iy))
1242  do ll = l, nanc
1243  if ( adjt(anc(ll))/=0 ) then
1244  k = k + 1
1245  invp(iperm(anc(ll))) = k
1246  endif
1247  enddo
1248  exit loop2
1249  else
1250  locc = btree(2,joc)
1251  if ( locc/=0 ) exit
1252  joc = parent(joc)
1253  endif
1254  enddo
1255  endif
1256  enddo loop2
1257  endif
1258  !
1259  ! set iperm
1260  !
1261  do i = 1, neqns
1262  iperm(invp(i)) = i
1263  enddo
1264 
1265  if ( idbg_sym/=0 ) write (6,"(10I6)") (invp(i),i=1,neqns)
1266  return
1267  else
1268  locc = btree(2,joc)
1269  if ( locc/=0 ) exit
1270  joc = parent(joc)
1271  endif
1272  enddo
1273  endif
1274  enddo
1275  exit
1276  enddo loop1
1277 
1278  deallocate (anc)
1279  deallocate (adjt)
1280  end subroutine rotate
1281 
1282  !======================================================================!
1284  ! sets DIAg, ZLN, DSLn
1285  ! updates RALoc, STAge
1286  !======================================================================!
1287  subroutine nuform(hecMESH,hecMAT,FCT,Ir)
1288  implicit none
1289  !------
1290  type (HECMWST_LOCAL_MESH), intent(in)::hecMESH
1291  type (HECMWST_MATRIX), intent(in)::hecMAT
1292  type (cholesky_factor), intent(inout):: FCT
1293  integer(kind=kint), intent(out):: Ir
1294  !------
1295  integer(kind=kint):: i
1296  integer(kind=kint):: idbg
1297  integer(kind=kint):: ierr
1298  integer(kind=kint):: j
1299  integer(kind=kint):: k
1300  integer(kind=kint):: kk
1301  integer(kind=kint):: ntotal
1302  integer(kind=kint):: numnp
1303  integer(kind=kint):: ndof
1304  integer(kind=kint):: ndof2
1305  real(kind=kreal), allocatable :: val(:)
1306 
1307  idbg = 0
1308  numnp = hecmat%NP
1309  ndof = hecmesh%N_DOF
1310  ntotal = numnp*ndof
1311 
1312  !*Allocations
1313  allocate (val(fct%NDEg*fct%NDEg),stat=ierr)
1314  if ( ierr/=0 ) stop "Allocation error:val"
1315  if ( idbg_num/= 0 ) write (6,*) "nuform:stage = ", fct%STAge
1316  kk = 0
1317  ndof2 = ndof*ndof
1318 
1319  do j = 1, numnp
1320  !*Diagonal
1321  kk = kk + 1
1322  call vlcpy(val,hecmat%D(ndof2*(j-1)+1:ndof2*j),ndof)
1323  call staij1(0,j,j,val,fct,ir)
1324 
1325  do i = 1, ndof
1326  if ( val((i-1)*ndof+i)<=0 ) write (idbg,*) 'j,j,val:', j, i, val((i-1)*ndof+i)
1327  enddo
1328 
1329  !*Lower
1330  do k = hecmat%INDEXL(j-1) + 1, hecmat%INDEXL(j)
1331  i = hecmat%ITEML(k)
1332  kk = kk + 1
1333  call vlcpy(val,hecmat%AL(ndof2*(k-1)+1:ndof2*k),ndof)
1334  call staij1(0,j,i,val,fct,ir)
1335  enddo
1336  enddo
1337 
1338  deallocate (val)
1339  end subroutine nuform
1340 
1341  !======================================================================!
1343  !======================================================================!
1344  subroutine vlcpy(A,B,N)
1345  implicit none
1346  !------
1347  integer(kind=kint), intent(in):: N
1348  real(kind=kreal), intent(in):: b(n,n)
1349  real(kind=kreal), intent(out):: a(n,n)
1350  !------
1351  integer(kind=kint):: i
1352 
1353  do i = 1, n
1354  a(1:n,i) = b(i,1:n)
1355  enddo
1356  end subroutine vlcpy
1357 
1358  !======================================================================!
1360  ! (symmetric version)
1361  ! (i)
1362  ! isw =0 set the value
1363  ! =1 add the value
1364  ! i row entry
1365  ! j column entry
1366  ! aij value
1367  ! (o)
1368  ! iv communication array
1369  ! sets DIAg, ZLN, DSLn
1370  ! updates RALoc, STAge
1371  !======================================================================!
1372  subroutine staij1(Isw,I,J,Aij,FCT,Ir)
1373  implicit none
1374  !------
1375  integer(kind=kint), intent(in):: I
1376  integer(kind=kint), intent(in):: Isw
1377  integer(kind=kint), intent(in):: J
1378  real(kind=kreal), intent(in):: aij(:)
1379  type (cholesky_factor), intent(inout):: FCT
1380  integer(kind=kint), intent(out):: Ir
1381  !------
1382  integer(kind=kint):: ndeg2
1383  integer(kind=kint):: ndeg2l
1384  integer(kind=kint):: ierror
1385 
1386  ir = 0
1387  ndeg2 = fct%NDEg*fct%NDEg
1388  ndeg2l = fct%NDEg*(fct%NDEg+1)/2
1389  if ( fct%STAge==30 ) write (6,*) 'warning a matrix was build up '//'but never solved.'
1390  if ( fct%STAge==10 ) then
1391  allocate (fct%DIAg(fct%NEQns*ndeg2l),stat=ierror)
1392  if ( ierror/=0 ) stop "Allocation error diag"
1393  fct%RALoc = fct%RALoc + fct%NEQns*ndeg2l
1394 
1395  allocate (fct%ZLN(fct%LEN_colno*ndeg2),stat=ierror)
1396  if ( ierror/=0 ) stop "Allocation error zln"
1397  fct%RALoc = fct%RALoc + fct%LEN_colno*ndeg2
1398 
1399  allocate (fct%DSLn(fct%LEN_dsln*ndeg2),stat=ierror)
1400  if ( ierror/=0 ) stop "Allocation error dsln"
1401  fct%RALoc = fct%RALoc + fct%LEN_dsln*ndeg2
1402  endif
1403  if ( fct%STAge/=20 ) then
1404  !
1405  ! for diagonal
1406  !
1407  fct%DIAg = 0.
1408  !
1409  ! for lower triangle
1410  !
1411  fct%ZLN = 0.
1412  !
1413  ! for dense window
1414  !
1415  fct%DSLn = 0.
1416 
1417  fct%STAge = 20
1418  endif
1419  ! Print *,'********Set Stage 20 *********'
1420  !
1421  if ( fct%NDEg<=2 ) then
1422  call addr0(isw,i,j,aij,fct%INVp,fct%XLNzr,fct%COLno,fct%DIAg,fct%ZLN,fct%DSLn,fct%NSTop,ndeg2,ndeg2l,ir)
1423  elseif ( fct%NDEg==3 ) then
1424  call addr3(i,j,aij,fct%INVp,fct%XLNzr,fct%COLno,fct%DIAg,fct%ZLN,fct%DSLn,fct%NSTop,ir)
1425  else
1426  call addrx(i,j,aij,fct%INVp,fct%XLNzr,fct%COLno,fct%DIAg,fct%ZLN,fct%DSLn,fct%NSTop,fct%NDEg,ndeg2l,ir)
1427  endif
1428  end subroutine staij1
1429 
1430  !======================================================================!
1432  !======================================================================!
1433  subroutine addr0(Isw,I,J,Aij,Invp,Xlnzr,Colno,Diag,Zln,Dsln,Nstop,Ndeg2,Ndeg2l,Ir)
1434  implicit none
1435  !------
1436  integer(kind=kint), intent(in):: I
1437  integer(kind=kint), intent(in):: Isw
1438  integer(kind=kint), intent(in):: J
1439  integer(kind=kint), intent(in):: Ndeg2
1440  integer(kind=kint), intent(in):: Ndeg2l
1441  integer(kind=kint), intent(in):: Nstop
1442  integer(kind=kint), intent(in):: Invp(:)
1443  integer(kind=kint), intent(in):: Xlnzr(:)
1444  integer(kind=kint), intent(in):: Colno(:)
1445  real(kind=kreal), intent(in):: aij(ndeg2)
1446  integer(kind=kint), intent(out):: Ir
1447  real(kind=kreal), intent(inout):: zln(ndeg2,*)
1448  real(kind=kreal), intent(inout):: diag(ndeg2l,*)
1449  real(kind=kreal), intent(inout):: dsln(ndeg2,*)
1450  !------
1451  integer(kind=kint):: i0
1452  integer(kind=kint):: ii
1453  integer(kind=kint):: itrans
1454  integer(kind=kint):: j0
1455  integer(kind=kint):: jj
1456  integer(kind=kint):: k
1457  integer(kind=kint):: ke
1458  integer(kind=kint):: ks
1459  integer(kind=kint), parameter:: idbg = 0
1460 
1461  ir = 0
1462  ii = invp(i)
1463  jj = invp(j)
1464  if ( idbg/=0 ) write (6,*) ii, jj, aij
1465  if ( ii==jj ) then
1466  if ( ndeg2==1 ) then
1467  if ( isw==0 ) then
1468  diag(1,ii) = aij(1)
1469  else
1470  diag(1,ii) = diag(1,ii) + aij(1)
1471  endif
1472  elseif ( ndeg2==4 ) then
1473  if ( isw==0 ) then
1474  diag(1,ii) = aij(1)
1475  diag(2,ii) = aij(2)
1476  diag(3,ii) = aij(4)
1477  else
1478  diag(1,ii) = diag(1,ii) + aij(1)
1479  diag(2,ii) = diag(2,ii) + aij(2)
1480  diag(3,ii) = diag(3,ii) + aij(4)
1481  endif
1482  endif
1483  return
1484  endif
1485  itrans = 0
1486  if ( jj>ii ) then
1487  k = jj
1488  jj = ii
1489  ii = k
1490  itrans = 1
1491  endif
1492  if ( jj>=nstop ) then
1493  i0 = ii-nstop
1494  j0 = jj-nstop + 1
1495  k = i0*(i0-1)/2 + j0
1496  if ( ndeg2==1 ) then
1497  dsln(1,k) = aij(1)
1498  return
1499  elseif ( ndeg2==4 ) then
1500  if ( itrans==0 ) then
1501  dsln(1:ndeg2,k) = aij(1:ndeg2)
1502  else
1503  dsln(1,k) = aij(1)
1504  dsln(2,k) = aij(3)
1505  dsln(3,k) = aij(2)
1506  dsln(4,k) = aij(4)
1507  endif
1508  return
1509  endif
1510  endif
1511  ks = xlnzr(ii)
1512  ke = xlnzr(ii+1) - 1
1513  do k = ks, ke
1514  if ( colno(k)==jj ) then
1515  if ( isw==0 ) then
1516  if ( ndeg2==1 ) then
1517  zln(1,k) = aij(1)
1518  elseif ( ndeg2==4 ) then
1519  if ( itrans==0 ) then
1520  zln(1:ndeg2,k) = aij(1:ndeg2)
1521  else
1522  zln(1,k) = aij(1)
1523  zln(2,k) = aij(3)
1524  zln(3,k) = aij(2)
1525  zln(4,k) = aij(4)
1526  endif
1527  endif
1528  elseif ( ndeg2==1 ) then
1529  zln(1,k) = zln(1,k) + aij(1)
1530  elseif ( ndeg2==4 ) then
1531  if ( itrans==0 ) then
1532  zln(1:ndeg2,k) = zln(1:ndeg2,k) + aij(1:ndeg2)
1533  else
1534  zln(1,k) = zln(1,k) + aij(1)
1535  zln(2,k) = zln(2,k) + aij(3)
1536  zln(3,k) = zln(3,k) + aij(2)
1537  zln(4,k) = zln(4,k) + aij(4)
1538  endif
1539  endif
1540  return
1541  endif
1542  enddo
1543  ir = 20
1544  end subroutine addr0
1545 
1546  !======================================================================!
1548  !======================================================================!
1549  subroutine addr3(I,J,Aij,Invp,Xlnzr,Colno,Diag,Zln,Dsln,Nstop,Ir)
1550  implicit none
1551  !------
1552  integer(kind=kint), intent(in):: I
1553  integer(kind=kint), intent(in):: J
1554  integer(kind=kint), intent(in):: Nstop
1555  integer(kind=kint), intent(in):: Invp(:)
1556  integer(kind=kint), intent(in):: Xlnzr(:)
1557  integer(kind=kint), intent(in):: Colno(:)
1558  real(kind=kreal), intent(in):: aij(9)
1559  integer(kind=kint), intent(out):: Ir
1560  real(kind=kreal), intent(inout):: zln(9,*)
1561  real(kind=kreal), intent(inout):: diag(6,*)
1562  real(kind=kreal), intent(inout):: dsln(9,*)
1563  !------
1564  integer(kind=kint):: i0
1565  integer(kind=kint):: ii
1566  integer(kind=kint):: itrans
1567  integer(kind=kint):: j0
1568  integer(kind=kint):: jj
1569  integer(kind=kint):: k
1570  integer(kind=kint):: ke
1571  integer(kind=kint):: ks
1572  integer(kind=kint):: l
1573  integer(kind=kint), parameter:: idbg = 0
1574  integer(kind=kint), parameter:: ndeg2 = 9
1575  integer(kind=kint), parameter:: ndeg2l = 6
1576 
1577  ir = 0
1578  ii = invp(i)
1579  jj = invp(j)
1580  if ( idbg/=0 ) write (6,*) ii, jj, aij
1581  if ( ii==jj ) then
1582  diag(1,ii) = aij(1)
1583  diag(2,ii) = aij(2)
1584  diag(3,ii) = aij(5)
1585  diag(4,ii) = aij(3)
1586  diag(5,ii) = aij(6)
1587  diag(6,ii) = aij(9)
1588  return
1589  endif
1590  itrans = 0
1591  if ( jj>ii ) then
1592  k = jj
1593  jj = ii
1594  ii = k
1595  itrans = 1
1596  endif
1597  if ( jj>=nstop ) then
1598  i0 = ii - nstop
1599  j0 = jj - nstop + 1
1600  k = i0*(i0-1)/2 + j0
1601  if ( itrans==0 ) then
1602  dsln(1:ndeg2,k) = aij(1:ndeg2)
1603  else
1604  dsln(1,k) = aij(1)
1605  dsln(2,k) = aij(4)
1606  dsln(3,k) = aij(7)
1607  dsln(4,k) = aij(2)
1608  dsln(5,k) = aij(5)
1609  dsln(6,k) = aij(8)
1610  dsln(7,k) = aij(3)
1611  dsln(8,k) = aij(6)
1612  dsln(9,k) = aij(9)
1613  endif
1614  return
1615  endif
1616  ks = xlnzr(ii)
1617  ke = xlnzr(ii+1) - 1
1618  do k = ks, ke
1619  if ( colno(k)==jj ) then
1620  if ( itrans==0 ) then
1621  do l = 1, ndeg2
1622  zln(l,k) = aij(l)
1623  enddo
1624  else
1625  zln(1,k) = aij(1)
1626  zln(2,k) = aij(4)
1627  zln(3,k) = aij(7)
1628  zln(4,k) = aij(2)
1629  zln(5,k) = aij(5)
1630  zln(6,k) = aij(8)
1631  zln(7,k) = aij(3)
1632  zln(8,k) = aij(6)
1633  zln(9,k) = aij(9)
1634  endif
1635  return
1636  endif
1637  enddo
1638  ir = 20
1639  end subroutine addr3
1640 
1641  !======================================================================!
1643  !======================================================================!
1644  subroutine addrx(I,J,Aij,Invp,Xlnzr,Colno,Diag,Zln,Dsln,Nstop,Ndeg,Ndeg2l,Ir)
1645  implicit none
1646  !------
1647  integer(kind=kint), intent(in):: I
1648  integer(kind=kint), intent(in):: J
1649  integer(kind=kint), intent(in):: Ndeg
1650  integer(kind=kint), intent(in):: Ndeg2l
1651  integer(kind=kint), intent(in):: Nstop
1652  integer(kind=kint), intent(in):: Invp(:)
1653  integer(kind=kint), intent(in):: Xlnzr(:)
1654  integer(kind=kint), intent(in):: Colno(:)
1655  real(kind=kreal), intent(in):: aij(ndeg,ndeg)
1656  integer(kind=kint), intent(out):: Ir
1657  real(kind=kreal), intent(inout):: zln(ndeg,ndeg,*)
1658  real(kind=kreal), intent(inout):: diag(ndeg2l,*)
1659  real(kind=kreal), intent(inout):: dsln(ndeg,ndeg,*)
1660  !------
1661  integer(kind=kint):: i0
1662  integer(kind=kint):: ii
1663  integer(kind=kint):: itrans
1664  integer(kind=kint):: j0
1665  integer(kind=kint):: jj
1666  integer(kind=kint):: k
1667  integer(kind=kint):: ke
1668  integer(kind=kint):: ks
1669  integer(kind=kint):: l
1670  integer(kind=kint):: m
1671  integer(kind=kint):: n
1672  integer(kind=kint), parameter:: idbg = 0
1673 
1674  ir = 0
1675  ii = invp(i)
1676  jj = invp(j)
1677  if ( idbg/=0 ) write (6,*) ii, jj, aij
1678  if ( ii==jj ) then
1679  l = 0
1680  do n = 1, ndeg
1681  do m = 1, n
1682  l = l + 1
1683  diag(l,ii) = aij(n,m)
1684  enddo
1685  enddo
1686  return
1687  endif
1688  itrans = 0
1689  if ( jj>ii ) then
1690  k = jj
1691  jj = ii
1692  ii = k
1693  itrans = 1
1694  endif
1695  if ( jj>=nstop ) then
1696  i0 = ii - nstop
1697  j0 = jj - nstop + 1
1698  k = i0*(i0-1)/2 + j0
1699  if ( itrans==0 ) then
1700  do m = 1, ndeg
1701  dsln(1:ndeg,m,k) = aij(1:ndeg,m)
1702  enddo
1703  else
1704  do m = 1, ndeg
1705  dsln(1:ndeg,m,k) = aij(m,1:ndeg)
1706  enddo
1707  endif
1708  return
1709  endif
1710  ks = xlnzr(ii)
1711  ke = xlnzr(ii+1) - 1
1712  do k = ks, ke
1713  if ( colno(k)==jj ) then
1714  if ( itrans==0 ) then
1715  do m = 1, ndeg
1716  zln(1:ndeg,m,k) = aij(1:ndeg,m)
1717  enddo
1718  else
1719  do m = 1, ndeg
1720  zln(1:ndeg,m,k) = aij(m,1:ndeg)
1721  enddo
1722  endif
1723  return
1724  endif
1725  enddo
1726  ir = 20
1727  end subroutine addrx
1728 
1729  !======================================================================!
1731  ! if(iv(22).eq.0) normal type
1732  ! if(iv(22).gt.0) code generation type
1733  ! updates NCH, DIAg, ZLN, DSLn, STAge
1734  !======================================================================!
1735  subroutine nufct0(FCT,Ir)
1736  implicit none
1737  !------
1738  type (cholesky_factor), intent(inout):: FCT
1739  integer(kind=kint), intent(out):: Ir
1740  !------
1741  integer(kind=kint):: irr
1742  integer(kind=kint):: ndeg2
1743  integer(kind=kint):: ndegl
1744  integer(kind=kint), allocatable :: INDx(:)
1745  real(kind=kreal), allocatable :: temp(:)
1746 
1747  if ( fct%STAge/=20 ) then
1748  print *, '*********Setting Stage 40!*********'
1749  ir = 40
1750  return
1751  else
1752  ir = 0
1753  endif
1754 
1755  allocate (temp(fct%NDEg*fct%NDEg*fct%NEQns),stat=irr)
1756  if ( irr/=0 ) then
1757  write (*,*) '##Error : Not enough memory'
1758  call hecmw_abort(hecmw_comm_get_comm())
1759  !stop
1760  endif
1761  allocate (indx(fct%NEQns),stat=irr)
1762  if ( irr/=0 ) then
1763  write (*,*) '##Error : Not enough memory'
1764  call hecmw_abort(hecmw_comm_get_comm())
1765  !stop
1766  endif
1767  !rmiv
1768  ndegl = fct%NDEg*(fct%NDEg+1)
1769  ndegl = ndegl/2
1770  ndeg2 = fct%NDEg*fct%NDEg
1771  !rmiv
1772  select case (fct%NDEg)
1773  case (1)
1774  call nufct(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,indx,temp,fct%NEQns,fct%PARent,fct%NCH,fct%NSTop,ir)
1775  case (2)
1776  call nufct2(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,indx,temp,fct%NEQns,fct%PARent,fct%NCH,fct%NSTop,ir)
1777  case (3)
1778  call nufct3(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,indx,temp,fct%NEQns,fct%PARent,fct%NCH,fct%NSTop,ir)
1779  case (6)
1780  call nufct6(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,indx,temp,fct%NEQns,fct%PARent,fct%NCH,fct%NSTop,ir)
1781  case default
1782  call nufctx(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,indx,temp,fct%NEQns,fct%PARent,fct%NCH,fct%NSTop,&
1783  fct%NDEg,ndegl,ir)
1784  endselect
1785 
1786  fct%STAge = 30
1787  deallocate (temp)
1788  deallocate (indx)
1789  end subroutine nufct0
1790 
1791  !======================================================================!
1793  ! (i) xlnzr,colno,zln,diag
1794  ! symbolicaly factorized
1795  ! (o) zln,diag,dsln
1796  !======================================================================!
1797  subroutine nufct(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
1798  implicit none
1799  !------
1800  integer(kind=kint), intent(in):: Neqns
1801  integer(kind=kint), intent(in):: Nstop
1802  integer(kind=kint), intent(in):: Xlnzr(:)
1803  integer(kind=kint), intent(in):: Colno(:)
1804  integer(kind=kint), intent(in):: Parent(:)
1805  integer(kind=kint), intent(out):: Ir
1806  integer(kind=kint), intent(out):: Indx(:)
1807  integer(kind=kint), intent(inout):: Nch(:)
1808  real(kind=kreal), intent(inout):: zln(:)
1809  real(kind=kreal), intent(inout):: diag(:)
1810  real(kind=kreal), intent(out):: temp(:)
1811  real(kind=kreal), intent(inout):: dsln(:)
1812  !------
1813  integer(kind=kint):: ic
1814  integer(kind=kint):: l
1815  integer(kind=kint):: ISEm
1816  real(kind=kreal):: t1
1817  real(kind=kreal):: t2
1818  real(kind=kreal):: t3
1819  real(kind=kreal):: t4
1820  real(kind=kreal):: t5
1821  real(kind=kreal):: tt
1822 
1823  isem = 1
1824  !
1825  ! phase I
1826  !
1827  call ptime(t1)
1828  diag(1) = 1.0d0/diag(1)
1829  l = parent(1)
1830  nch(l) = nch(l) - 1
1831  nch(1) = -1
1832  do ic = 2, nstop - 1
1833  call sum(ic,xlnzr,colno,zln,diag,nch,parent,temp,indx,isem)
1834  enddo
1835  !
1836  ! phase II
1837  !
1838  call ptime(t2)
1839  do ic = nstop, neqns
1840  call sum1(ic,xlnzr,colno,zln,temp,indx)
1841  enddo
1842  !
1843  ! phase III
1844  !
1845  call ptime(t3)
1846  call sum2(neqns,nstop,xlnzr,colno,zln,diag,dsln,temp,indx)
1847  !
1848  ! phase IV
1849  !
1850  call ptime(t4)
1851  call sum3(neqns-nstop+1,dsln,diag(nstop:),indx,temp)
1852  call ptime(t5)
1853  tt = t5 - t1
1854  t1 = t2 - t1
1855  t2 = t3 - t2
1856  t3 = t4 - t3
1857  t4 = t5 - t4
1858  return
1859  ir = 30
1860  end subroutine nufct
1861 
1862  !======================================================================!
1864  ! (i) xlnzr,colno,zln,diag
1865  ! symbolicaly factorized
1866  ! (o) zln,diag,dsln
1867  !======================================================================!
1868  subroutine nufct2(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
1869  implicit none
1870  !------
1871  integer(kind=kint), intent(in):: Neqns
1872  integer(kind=kint), intent(in):: Nstop
1873  integer(kind=kint), intent(in):: Xlnzr(:)
1874  integer(kind=kint), intent(in):: Colno(:)
1875  integer(kind=kint), intent(in):: Parent(:)
1876  integer(kind=kint), intent(out):: Ir
1877  integer(kind=kint), intent(out):: Indx(:)
1878  integer(kind=kint), intent(inout):: Nch(:)
1879  real(kind=kreal), intent(inout):: zln(4,*)
1880  real(kind=kreal), intent(inout):: diag(3,*)
1881  real(kind=kreal), intent(out):: temp(4,*)
1882  real(kind=kreal), intent(inout):: dsln(4,*)
1883  !------
1884  integer(kind=kint):: ic
1885  integer(kind=kint):: l
1886  real(kind=kreal):: t1
1887  real(kind=kreal):: t2
1888  real(kind=kreal):: t3
1889  real(kind=kreal):: t4
1890  real(kind=kreal):: t5
1891  real(kind=kreal):: tt
1892 
1893  !
1894  ! phase I
1895  !
1896  ir = 0
1897  call ptime(t1)
1898  if ( nstop>1 ) call inv2(diag(1,1),ir)
1899  l = parent(1)
1900  nch(l) = nch(l) - 1
1901  nch(1) = -1
1902  do ic = 2, nstop - 1
1903  call s2um(ic,xlnzr,colno,zln,diag,nch,parent,temp,indx)
1904  enddo
1905  !
1906  ! phase II
1907  !
1908  call ptime(t2)
1909  do ic = nstop, neqns
1910  call s2um1(ic,xlnzr,colno,zln,temp,indx)
1911  enddo
1912  !
1913  ! phase III
1914  !
1915  call ptime(t3)
1916  call s2um2(neqns,nstop,xlnzr,colno,zln,diag,dsln,temp,indx)
1917  !
1918  ! phase IV
1919  !
1920  call ptime(t4)
1921  call s2um3(neqns-nstop+1,dsln,diag(1,nstop),indx,temp)
1922  call ptime(t5)
1923  tt = t5 - t1
1924  t1 = t2 - t1
1925  t2 = t3 - t2
1926  t3 = t4 - t3
1927  t4 = t5 - t4
1928  end subroutine nufct2
1929 
1930  !======================================================================!
1932  ! (i) xlnzr,colno,zln,diag
1933  ! symbolicaly factorized
1934  ! (o) zln,diag,dsln
1935  !======================================================================!
1936  subroutine nufct3(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
1937  implicit none
1938  !------
1939  integer(kind=kint), intent(in):: Neqns
1940  integer(kind=kint), intent(in):: Nstop
1941  integer(kind=kint), intent(in):: Xlnzr(:)
1942  integer(kind=kint), intent(in):: Colno(:)
1943  integer(kind=kint), intent(in):: Parent(:)
1944  integer(kind=kint), intent(out):: Ir
1945  integer(kind=kint), intent(out):: Indx(:)
1946  integer(kind=kint), intent(inout):: Nch(:)
1947  real(kind=kreal), intent(inout):: zln(9,*)
1948  real(kind=kreal), intent(inout):: diag(6,*)
1949  real(kind=kreal), intent(out):: temp(:)
1950  real(kind=kreal), intent(inout):: dsln(9,*)
1951  !------
1952  integer(kind=kint):: ic
1953  integer(kind=kint):: l
1954  real(kind=kreal):: t1
1955  real(kind=kreal):: t2
1956  real(kind=kreal):: t3
1957  real(kind=kreal):: t4
1958  real(kind=kreal):: t5
1959  real(kind=kreal):: tt
1960 
1961  !
1962  ! phase I
1963  !
1964  call ptime(t1)
1965  if ( nstop>1 ) call inv3(diag(1,1),ir)
1966  l = parent(1)
1967  nch(l) = nch(l) - 1
1968  nch(1) = -1
1969  do ic = 2, nstop - 1
1970  call s3um(ic,xlnzr,colno,zln,diag,nch,parent,temp,indx)
1971  enddo
1972  !
1973  ! phase II
1974  !
1975  call ptime(t2)
1976  do ic = nstop, neqns
1977  call s3um1(ic,xlnzr,colno,zln,temp,indx)
1978  enddo
1979  !
1980  ! phase III
1981  !
1982  call ptime(t3)
1983  call s3um2(neqns,nstop,xlnzr,colno,zln,diag,dsln,temp,indx)
1984  !
1985  ! phase IV
1986  !
1987  call ptime(t4)
1988  call s3um3(neqns-nstop+1,dsln,diag(1,nstop),indx,temp)
1989  call ptime(t5)
1990  tt = t5 - t1
1991  t1 = t2 - t1
1992  t2 = t3 - t2
1993  t3 = t4 - t3
1994  t4 = t5 - t4
1995  end subroutine nufct3
1996 
1997  !======================================================================!
1999  ! (i) xlnzr,colno,zln,diag
2000  ! symbolicaly factorized
2001  ! (o) zln,diag,dsln
2002  !======================================================================!
2003  subroutine nufct6(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ir)
2004  implicit none
2005  !------
2006  integer(kind=kint), intent(in):: Neqns
2007  integer(kind=kint), intent(in):: Nstop
2008  integer(kind=kint), intent(in):: Xlnzr(:)
2009  integer(kind=kint), intent(in):: Colno(:)
2010  integer(kind=kint), intent(in):: Parent(:)
2011  integer(kind=kint), intent(out):: Indx(:)
2012  integer(kind=kint), intent(inout):: Nch(:)
2013  real(kind=kreal), intent(inout):: zln(36,*)
2014  real(kind=kreal), intent(inout):: diag(21,*)
2015  real(kind=kreal), intent(out):: temp(:)
2016  real(kind=kreal), intent(inout):: dsln(36,*)
2017  !------
2018  integer(kind=kint):: ic
2019  integer(kind=kint):: Ir
2020  integer(kind=kint):: l
2021  real(kind=kreal):: t1
2022  real(kind=kreal):: t2
2023  real(kind=kreal):: t3
2024  real(kind=kreal):: t4
2025  real(kind=kreal):: t5
2026  real(kind=kreal):: tt
2027 
2028  !
2029  ! phase I
2030  !
2031  call ptime(t1)
2032  if ( nstop>1 ) call inv6(diag(1,1),ir)
2033  l = parent(1)
2034  nch(l) = nch(l) - 1
2035  nch(1) = -1
2036  do ic = 2, nstop - 1
2037  call s6um(ic,xlnzr,colno,zln,diag,nch,parent,temp,indx)
2038  enddo
2039  !
2040  ! phase II
2041  !
2042  call ptime(t2)
2043  do ic = nstop, neqns
2044  call s6um1(ic,xlnzr,colno,zln,temp,indx)
2045  enddo
2046  !
2047  ! phase III
2048  !
2049  call ptime(t3)
2050  call s6um2(neqns,nstop,xlnzr,colno,zln,diag,dsln,temp,indx)
2051  !
2052  ! phase IV
2053  !
2054  call ptime(t4)
2055  call s6um3(neqns-nstop+1,dsln,diag(1,nstop),indx,temp)
2056  call ptime(t5)
2057  tt = t5 - t1
2058  t1 = t2 - t1
2059  t2 = t3 - t2
2060  t3 = t4 - t3
2061  t4 = t5 - t4
2062  end subroutine nufct6
2063 
2064  !======================================================================!
2066  ! (i) xlnzr,colno,zln,diag
2067  ! symbolicaly factorized
2068  ! (o) zln,diag,dsln
2069  !======================================================================!
2070  subroutine nufctx(Xlnzr,Colno,Dsln,Zln,Diag,Indx,Temp,Neqns,Parent,Nch,Nstop,Ndeg,Ndegl,Ir)
2071  implicit none
2072  !------
2073  integer(kind=kint), intent(in):: Ndeg
2074  integer(kind=kint), intent(in):: Ndegl
2075  integer(kind=kint), intent(in):: Neqns
2076  integer(kind=kint), intent(in):: Nstop
2077  integer(kind=kint), intent(in):: Xlnzr(:)
2078  integer(kind=kint), intent(in):: Colno(:)
2079  integer(kind=kint), intent(in):: Parent(:)
2080  integer(kind=kint), intent(out):: Indx(:)
2081  integer(kind=kint), intent(inout):: Nch(:)
2082  real(kind=kreal), intent(inout):: zln(ndeg*ndeg,*)
2083  real(kind=kreal), intent(inout):: diag(ndegl,*)
2084  real(kind=kreal), intent(out):: temp(ndeg*ndeg,*)
2085  real(kind=kreal), intent(inout):: dsln(ndeg*ndeg,*)
2086  !------
2087  integer(kind=kint):: ic
2088  integer(kind=kint):: Ir
2089  integer(kind=kint):: l
2090  real(kind=kreal):: t1
2091  real(kind=kreal):: t2
2092  real(kind=kreal):: t3
2093  real(kind=kreal):: t4
2094  real(kind=kreal):: t5
2095  real(kind=kreal):: tt
2096  real(kind=kreal):: zz(100)
2097  real(kind=kreal):: t(100)
2098 
2099  !
2100  ! phase I
2101  !
2102  call ptime(t1)
2103  if ( nstop>1 ) call invx(diag(1,1),ndeg,ir)
2104  l = parent(1)
2105  nch(l) = nch(l) - 1
2106  nch(1) = -1
2107  do ic = 2, nstop - 1
2108  call sxum(ic,xlnzr,colno,zln,diag,nch,parent,temp,indx,ndeg,ndegl,zz,t)
2109  enddo
2110  !
2111  ! phase II
2112  !
2113  call ptime(t2)
2114  do ic = nstop, neqns
2115  call sxum1(ic,xlnzr,colno,zln,temp,indx,ndeg,t)
2116  enddo
2117  !
2118  ! phase III
2119  !
2120  call ptime(t3)
2121  call sxum2(neqns,nstop,xlnzr,colno,zln,diag,dsln,temp,indx,ndeg,ndegl)
2122  !
2123  ! phase IV
2124  !
2125  call ptime(t4)
2126  call sxum3(neqns-nstop+1,dsln,diag(1,nstop),indx,temp,ndeg,ndegl,t)
2127  call ptime(t5)
2128  tt = t5 - t1
2129  t1 = t2 - t1
2130  t2 = t3 - t2
2131  t3 = t4 - t3
2132  t4 = t5 - t4
2133  end subroutine nufctx
2134 
2135  !======================================================================!
2137  !======================================================================!
2138  subroutine sum(Ic,Xlnzr,Colno,Zln,Diag,Nch,Par,Temp,Indx,ISEm)
2139  implicit none
2140  !------
2141  integer(kind=kint), intent(in):: Ic
2142  integer(kind=kint), intent(in):: Xlnzr(:)
2143  integer(kind=kint), intent(in):: Colno(:)
2144  integer(kind=kint), intent(in):: Par(:)
2145  integer(kind=kint), intent(inout):: Indx(:)
2146  integer(kind=kint), intent(inout):: Nch(:)
2147  integer(kind=kint), intent(inout):: ISEm
2148  real(kind=kreal), intent(inout):: diag(:)
2149  real(kind=kreal), intent(inout):: temp(:)
2150  real(kind=kreal), intent(inout):: zln(:)
2151  !------
2152  integer(kind=kint):: j
2153  integer(kind=kint):: jc
2154  integer(kind=kint):: jj
2155  integer(kind=kint):: k
2156  integer(kind=kint):: ke
2157  integer(kind=kint):: kk
2158  integer(kind=kint):: ks
2159  real(kind=kreal):: piv
2160  real(kind=kreal):: s
2161  real(kind=kreal):: t
2162  real(kind=kreal):: zz
2163 
2164  ks = xlnzr(ic)
2165  ke = xlnzr(ic+1)
2166  t = 0.0d0
2167 
2168  do k = ks, ke - 1
2169  jc = colno(k)
2170  indx(jc) = ic
2171  s = 0.0d0
2172  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2173  j = colno(jj)
2174  if ( indx(j)==ic ) s = s + temp(j)*zln(jj)
2175  enddo
2176  zz = zln(k) - s
2177  zln(k) = zz*diag(jc)
2178  temp(jc) = zz
2179  t = t + zz*zln(k)
2180  enddo
2181  piv = diag(ic) - t
2182  if ( dabs(piv)>rmin ) diag(ic) = 1.0d0/piv
2183  do while ( isem/=1 )
2184  enddo
2185  isem = 0
2186  nch(ic) = -1
2187  kk = par(ic)
2188  nch(kk) = nch(kk) - 1
2189  isem = 1
2190  end subroutine sum
2191 
2192  !======================================================================!
2194  !======================================================================!
2195  subroutine sum1(Ic,Xlnzr,Colno,Zln,Temp,Indx)
2196  implicit none
2197  !------
2198  integer(kind=kint), intent(in):: Ic
2199  integer(kind=kint), intent(in):: Xlnzr(:)
2200  integer(kind=kint), intent(in):: Colno(:)
2201  integer(kind=kint), intent(inout):: Indx(:)
2202  real(kind=kreal), intent(inout):: temp(:)
2203  real(kind=kreal), intent(inout):: zln(:)
2204  !------
2205  integer(kind=kint):: j
2206  integer(kind=kint):: jc
2207  integer(kind=kint):: jj
2208  integer(kind=kint):: k
2209  integer(kind=kint):: ke
2210  integer(kind=kint):: ks
2211  real(kind=kreal):: s
2212  real(kind=kreal):: t
2213  real(kind=kreal):: zz
2214 
2215  ks = xlnzr(ic)
2216  ke = xlnzr(ic+1)
2217  t = 0.0d0
2218 
2219  do k = ks, ke - 1
2220  jc = colno(k)
2221  indx(jc) = ic
2222  s = 0.0d0
2223  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2224  j = colno(jj)
2225  if ( indx(j)==ic ) s = s + temp(j)*zln(jj)
2226  enddo
2227  zz = zln(k) - s
2228 
2229  zln(k) = zz
2230  temp(jc) = zz
2231  enddo
2232  end subroutine sum1
2233 
2234  !======================================================================!
2236  !======================================================================!
2237  subroutine sum2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
2238  implicit none
2239  !------
2240  integer(kind=kint), intent(in):: Neqns
2241  integer(kind=kint), intent(in):: Nstop
2242  integer(kind=kint), intent(in):: Xlnzr(:)
2243  integer(kind=kint), intent(in):: Colno(:)
2244  integer(kind=kint), intent(inout):: Indx(:)
2245  real(kind=kreal), intent(inout):: diag(:)
2246  real(kind=kreal), intent(inout):: dsln(:)
2247  real(kind=kreal), intent(inout):: temp(:)
2248  real(kind=kreal), intent(inout):: zln(:)
2249  !------
2250  integer(kind=kint):: ic
2251  integer(kind=kint):: j
2252  integer(kind=kint):: jc
2253  integer(kind=kint):: jj
2254  integer(kind=kint):: joc
2255  integer(kind=kint):: k
2256  integer(kind=kint):: ke
2257  integer(kind=kint):: ks
2258  real(kind=kreal):: s
2259 
2260  joc = 0
2261  do ic = nstop, neqns
2262  temp(1:nstop) = 0.0d0
2263  ks = xlnzr(ic)
2264  ke = xlnzr(ic+1) - 1
2265  do k = ks, ke
2266  jj = colno(k)
2267  temp(jj) = zln(k)
2268  zln(k) = temp(jj)*diag(jj)
2269  indx(jj) = ic
2270  diag(ic) = diag(ic) - temp(jj)*zln(k)
2271  enddo
2272  do jc = nstop, ic - 1
2273  s = 0.0d0
2274  joc = joc + 1
2275  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2276  j = colno(jj)
2277  if ( indx(j)==ic ) s = s + temp(j)*zln(jj)
2278  enddo
2279  if ( s==0.0d0 ) write (16,*) ic, jc
2280  dsln(joc) = dsln(joc) - s
2281  enddo
2282  enddo
2283  end subroutine sum2
2284 
2285  !======================================================================!
2287  !======================================================================!
2288  subroutine sum3(N,Dsln,Diag,Indx,Temp)
2289  implicit none
2290  !------
2291  integer(kind=kint), intent(in):: N
2292  integer(kind=kint), intent(out):: Indx(:)
2293  real(kind=kreal), intent(inout):: diag(:)
2294  real(kind=kreal), intent(inout):: dsln(:)
2295  real(kind=kreal), intent(out):: temp(:)
2296  !------
2297  integer(kind=kint):: i
2298  integer(kind=kint):: j
2299  integer(kind=kint):: joc
2300 
2301  if ( n>0 ) then
2302  indx(1) = 0
2303  joc = 1
2304  diag(1) = 1.0d0/diag(1)
2305  do i = 2, n
2306  indx(i) = joc
2307  do j = 1, i - 1
2308  dsln(joc) = dsln(joc) - ddot(dsln(indx(i):),dsln(indx(j):),j-1)
2309  joc = joc + 1
2310  enddo
2311  call vprod(dsln(indx(i):),diag,temp,i-1)
2312  diag(i) = diag(i) - ddot(temp,dsln(indx(i):),i-1)
2313  call vcopy(temp,dsln(indx(i):),i-1)
2314  diag(i) = 1.0d0/diag(i)
2315  enddo
2316  endif
2317  end subroutine sum3
2318 
2319  !======================================================================!
2321  !======================================================================!
2322  subroutine s2um(Ic,Xlnzr,Colno,Zln,Diag,Nch,Par,Temp,Indx)
2323  implicit none
2324  !------
2325  integer(kind=kint), intent(in):: Ic
2326  integer(kind=kint), intent(in):: Xlnzr(:)
2327  integer(kind=kint), intent(in):: Colno(:)
2328  integer(kind=kint), intent(in):: Par(:)
2329  integer(kind=kint), intent(inout):: Indx(:)
2330  integer(kind=kint), intent(inout):: Nch(:)
2331  real(kind=kreal), intent(inout):: diag(3,*)
2332  real(kind=kreal), intent(inout):: temp(4,*)
2333  real(kind=kreal), intent(inout):: zln(4,*)
2334  !------
2335  integer(kind=kint):: ir
2336  integer(kind=kint):: j
2337  integer(kind=kint):: jc
2338  integer(kind=kint):: jj
2339  integer(kind=kint):: k
2340  integer(kind=kint):: ke
2341  integer(kind=kint):: kk
2342  integer(kind=kint):: ks
2343  real(kind=kreal):: s(4)
2344  real(kind=kreal):: t(3)
2345  real(kind=kreal):: zz(4)
2346 
2347  ks = xlnzr(ic)
2348  ke = xlnzr(ic+1)
2349  t(1:3) = 0.0d0
2350  do k = ks, ke - 1
2351  jc = colno(k)
2352  indx(jc) = ic
2353  s(1:4) = 0.0d0
2354  zz(1:4) = zln(1:4,k)
2355  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2356  j = colno(jj)
2357  if ( indx(j)==ic ) then
2358  zz(1) = zz(1) - temp(1,j)*zln(1,jj) - temp(3,j)*zln(3,jj)
2359  zz(2) = zz(2) - temp(2,j)*zln(1,jj) - temp(4,j)*zln(3,jj)
2360  zz(3) = zz(3) - temp(1,j)*zln(2,jj) - temp(3,j)*zln(4,jj)
2361  zz(4) = zz(4) - temp(2,j)*zln(2,jj) - temp(4,j)*zln(4,jj)
2362  endif
2363  enddo
2364  call inv22(zln(1,k),zz,diag(1,jc))
2365  temp(1:4,jc) = zz(1:4)
2366  t(1) = t(1) + zz(1)*zln(1,k) + zz(3)*zln(3,k)
2367  t(2) = t(2) + zz(1)*zln(2,k) + zz(3)*zln(4,k)
2368  t(3) = t(3) + zz(2)*zln(2,k) + zz(4)*zln(4,k)
2369  enddo
2370  diag(1,ic) = diag(1,ic) - t(1)
2371  diag(2,ic) = diag(2,ic) - t(2)
2372  diag(3,ic) = diag(3,ic) - t(3)
2373  call inv2(diag(1,ic),ir)
2374  nch(ic) = -1
2375  kk = par(ic)
2376  nch(kk) = nch(kk) - 1
2377  end subroutine s2um
2378 
2379  !======================================================================!
2381  !======================================================================!
2382  subroutine s2um1(Ic,Xlnzr,Colno,Zln,Temp,Indx)
2383  implicit none
2384  !------
2385  integer(kind=kint), intent(in):: Ic
2386  integer(kind=kint), intent(in):: Xlnzr(:)
2387  integer(kind=kint), intent(in):: Colno(:)
2388  integer(kind=kint), intent(inout):: Indx(:)
2389  real(kind=kreal), intent(inout):: temp(4,*)
2390  real(kind=kreal), intent(inout):: zln(4,*)
2391  !------
2392  integer(kind=kint):: j
2393  integer(kind=kint):: jc
2394  integer(kind=kint):: jj
2395  integer(kind=kint):: k
2396  integer(kind=kint):: ke
2397  integer(kind=kint):: ks
2398  integer(kind=kint):: l
2399  real(kind=kreal):: s(4)
2400 
2401  ks = xlnzr(ic)
2402  ke = xlnzr(ic+1)
2403  s(1:4) = 0.0d0
2404  do k = ks, ke - 1
2405  jc = colno(k)
2406  indx(jc) = ic
2407  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2408  j = colno(jj)
2409  if ( indx(j)==ic ) then
2410  s(1) = s(1) + temp(1,j)*zln(1,jj) + temp(3,j)*zln(3,jj)
2411  s(2) = s(2) + temp(2,j)*zln(1,jj) + temp(4,j)*zln(3,jj)
2412  s(3) = s(3) + temp(1,j)*zln(2,jj) + temp(3,j)*zln(4,jj)
2413  s(4) = s(4) + temp(2,j)*zln(2,jj) + temp(4,j)*zln(4,jj)
2414  endif
2415  enddo
2416  do l = 1, 4
2417  temp(l,jc) = zln(l,k) - s(l)
2418  zln(l,k) = temp(l,jc)
2419  s(l) = 0.0d0
2420  enddo
2421  enddo
2422  end subroutine s2um1
2423 
2424  !======================================================================!
2426  !======================================================================!
2427  subroutine s2um2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
2428  implicit none
2429  !------
2430  integer(kind=kint), intent(in):: Neqns
2431  integer(kind=kint), intent(in):: Nstop
2432  integer(kind=kint), intent(in):: Xlnzr(:)
2433  integer(kind=kint), intent(in):: Colno(:)
2434  integer(kind=kint), intent(inout):: Indx(:)
2435  real(kind=kreal), intent(inout):: diag(3,*)
2436  real(kind=kreal), intent(inout):: dsln(4,*)
2437  real(kind=kreal), intent(inout):: temp(4,*)
2438  real(kind=kreal), intent(inout):: zln(4,*)
2439  !------
2440  integer(kind=kint):: ic
2441  integer(kind=kint):: j
2442  integer(kind=kint):: jc
2443  integer(kind=kint):: jj
2444  integer(kind=kint):: joc
2445  integer(kind=kint):: k
2446  integer(kind=kint):: ke
2447  integer(kind=kint):: ks
2448 
2449  joc = 0
2450  do ic = nstop, neqns
2451  ks = xlnzr(ic)
2452  ke = xlnzr(ic+1) - 1
2453  do k = ks, ke
2454  jj = colno(k)
2455  temp(1,jj) = zln(1,k)
2456  temp(2,jj) = zln(2,k)
2457  temp(3,jj) = zln(3,k)
2458  temp(4,jj) = zln(4,k)
2459 
2460  zln(3,k) = temp(3,jj) - temp(1,jj)*diag(2,jj)
2461  zln(1,k) = temp(1,jj)*diag(1,jj)
2462  zln(3,k) = zln(3,k)*diag(3,jj)
2463  zln(1,k) = zln(1,k) - zln(3,k)*diag(2,jj)
2464 
2465  zln(4,k) = temp(4,jj) - temp(2,jj)*diag(2,jj)
2466  zln(2,k) = temp(2,jj)*diag(1,jj)
2467  zln(4,k) = zln(4,k)*diag(3,jj)
2468  zln(2,k) = zln(2,k) - zln(4,k)*diag(2,jj)
2469 
2470  diag(1,ic) = diag(1,ic) - (temp(1,jj)*zln(1,k)+temp(3,jj)*zln(3,k))
2471  diag(2,ic) = diag(2,ic) - (temp(1,jj)*zln(2,k)+temp(3,jj)*zln(4,k))
2472  diag(3,ic) = diag(3,ic) - (temp(2,jj)*zln(2,k)+temp(4,jj)*zln(4,k))
2473  indx(jj) = ic
2474  enddo
2475  do jc = nstop, ic - 1
2476  joc = joc + 1
2477  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2478  j = colno(jj)
2479  if ( indx(j)==ic ) then
2480  dsln(1,joc) = dsln(1,joc) - (temp(1,j)*zln(1,jj)+temp(3,j)*zln(3,jj))
2481  dsln(2,joc) = dsln(2,joc) - (temp(2,j)*zln(1,jj)+temp(4,j)*zln(3,jj))
2482  dsln(3,joc) = dsln(3,joc) - (temp(1,j)*zln(2,jj)+temp(3,j)*zln(4,jj))
2483  dsln(4,joc) = dsln(4,joc) - (temp(2,j)*zln(2,jj)+temp(4,j)*zln(4,jj))
2484  endif
2485  enddo
2486  enddo
2487  enddo
2488  end subroutine s2um2
2489 
2490  !======================================================================!
2492  !======================================================================!
2493  subroutine s2um3(N,Dsln,Diag,Indx,Temp)
2494  implicit none
2495  !------
2496  integer(kind=kint), intent(in):: N
2497  integer(kind=kint), intent(out):: Indx(:)
2498  real(kind=kreal), intent(inout):: diag(3,*)
2499  real(kind=kreal), intent(inout):: dsln(4,*)
2500  real(kind=kreal), intent(out):: temp(4,*)
2501  !------
2502  integer(kind=kint):: i
2503  integer(kind=kint):: ir
2504  integer(kind=kint):: j
2505  integer(kind=kint):: joc
2506  real(kind=kreal):: t(4)
2507 
2508  if ( n>0 ) then
2509  indx(1) = 0
2510  joc = 1
2511  call inv2(diag(1,1),ir)
2512  do i = 2, n
2513  indx(i) = joc
2514  do j = 1, i - 1
2515  call d2dot(t,dsln(1,indx(i)),dsln(1,indx(j)),j-1)
2516  dsln(1,joc) = dsln(1,joc) - t(1)
2517  dsln(2,joc) = dsln(2,joc) - t(2)
2518  dsln(3,joc) = dsln(3,joc) - t(3)
2519  dsln(4,joc) = dsln(4,joc) - t(4)
2520  joc = joc + 1
2521  enddo
2522  call v2prod(dsln(1,indx(i)),diag,temp,i-1)
2523  call d2dot(t,temp,dsln(1,indx(i)),i-1)
2524  diag(1,i) = diag(1,i) - t(1)
2525  diag(2,i) = diag(2,i) - t(2)
2526  diag(3,i) = diag(3,i) - t(4)
2527  call vcopy(temp,dsln(1,indx(i)),4*(i-1))
2528  call inv2(diag(1,i),ir)
2529  enddo
2530  endif
2531  end subroutine s2um3
2532 
2533  !======================================================================!
2535  !======================================================================!
2536  subroutine s3um(Ic,Xlnzr,Colno,Zln,Diag,Nch,Par,Temp,Indx)
2537  implicit none
2538  !------
2539  integer(kind=kint), intent(in):: Ic
2540  integer(kind=kint), intent(in):: Xlnzr(:)
2541  integer(kind=kint), intent(in):: Colno(:)
2542  integer(kind=kint), intent(in):: Par(:)
2543  integer(kind=kint), intent(inout):: Indx(:)
2544  integer(kind=kint), intent(inout):: Nch(:)
2545  real(kind=kreal), intent(inout):: diag(6,*)
2546  real(kind=kreal), intent(inout):: temp(9,*)
2547  real(kind=kreal), intent(inout):: zln(9,*)
2548  !------
2549  integer(kind=kint):: ir
2550  integer(kind=kint):: j
2551  integer(kind=kint):: jc
2552  integer(kind=kint):: jj
2553  integer(kind=kint):: k
2554  integer(kind=kint):: ke
2555  integer(kind=kint):: kk
2556  integer(kind=kint):: ks
2557  integer(kind=kint):: l
2558  real(kind=kreal):: t(6)
2559  real(kind=kreal):: zz(9)
2560 
2561  ks = xlnzr(ic)
2562  ke = xlnzr(ic+1)
2563  t(1:6) = 0.0d0
2564  do k = ks, ke - 1
2565  jc = colno(k)
2566  indx(jc) = ic
2567  zz(1:9) = zln(1:9,k)
2568  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2569  j = colno(jj)
2570  if ( indx(j)==ic ) then
2571  zz(1) = zz(1) - temp(1,j)*zln(1,jj) - temp(4,j)*zln(4,jj) - temp(7,j)*zln(7,jj)
2572  zz(2) = zz(2) - temp(2,j)*zln(1,jj) - temp(5,j)*zln(4,jj) - temp(8,j)*zln(7,jj)
2573  zz(3) = zz(3) - temp(3,j)*zln(1,jj) - temp(6,j)*zln(4,jj) - temp(9,j)*zln(7,jj)
2574 
2575  zz(4) = zz(4) - temp(1,j)*zln(2,jj) - temp(4,j)*zln(5,jj) - temp(7,j)*zln(8,jj)
2576  zz(5) = zz(5) - temp(2,j)*zln(2,jj) - temp(5,j)*zln(5,jj) - temp(8,j)*zln(8,jj)
2577  zz(6) = zz(6) - temp(3,j)*zln(2,jj) - temp(6,j)*zln(5,jj) - temp(9,j)*zln(8,jj)
2578 
2579  zz(7) = zz(7) - temp(1,j)*zln(3,jj) - temp(4,j)*zln(6,jj) - temp(7,j)*zln(9,jj)
2580  zz(8) = zz(8) - temp(2,j)*zln(3,jj) - temp(5,j)*zln(6,jj) - temp(8,j)*zln(9,jj)
2581  zz(9) = zz(9) - temp(3,j)*zln(3,jj) - temp(6,j)*zln(6,jj) - temp(9,j)*zln(9,jj)
2582  endif
2583  enddo
2584  call inv33(zln(1,k),zz,diag(1,jc))
2585  temp(1:9,jc) = zz(1:9)
2586  t(1) = t(1) + zz(1)*zln(1,k) + zz(4)*zln(4,k) + zz(7)*zln(7,k)
2587  t(2) = t(2) + zz(1)*zln(2,k) + zz(4)*zln(5,k) + zz(7)*zln(8,k)
2588  t(3) = t(3) + zz(2)*zln(2,k) + zz(5)*zln(5,k) + zz(8)*zln(8,k)
2589  t(4) = t(4) + zz(1)*zln(3,k) + zz(4)*zln(6,k) + zz(7)*zln(9,k)
2590  t(5) = t(5) + zz(2)*zln(3,k) + zz(5)*zln(6,k) + zz(8)*zln(9,k)
2591  t(6) = t(6) + zz(3)*zln(3,k) + zz(6)*zln(6,k) + zz(9)*zln(9,k)
2592  enddo
2593  do l = 1, 6
2594  diag(l,ic) = diag(l,ic) - t(l)
2595  enddo
2596  call inv3(diag(1,ic),ir)
2597  nch(ic) = -1
2598  kk = par(ic)
2599  nch(kk) = nch(kk) - 1
2600  end subroutine s3um
2601 
2602  !======================================================================!
2604  !======================================================================!
2605  subroutine s3um1(Ic,Xlnzr,Colno,Zln,Temp,Indx)
2606  implicit none
2607  !------
2608  integer(kind=kint), intent(in):: Ic
2609  integer(kind=kint), intent(in):: Xlnzr(:)
2610  integer(kind=kint), intent(in):: Colno(:)
2611  integer(kind=kint), intent(inout):: Indx(:)
2612  real(kind=kreal), intent(inout):: temp(9,*)
2613  real(kind=kreal), intent(inout):: zln(9,*)
2614  !------
2615  integer(kind=kint):: j
2616  integer(kind=kint):: jc
2617  integer(kind=kint):: jj
2618  integer(kind=kint):: k
2619  integer(kind=kint):: ke
2620  integer(kind=kint):: ks
2621  integer(kind=kint):: l
2622  real(kind=kreal):: s(9)
2623 
2624  ks = xlnzr(ic)
2625  ke = xlnzr(ic+1)
2626  s(1:9) = 0.0d0
2627  do k = ks, ke - 1
2628  jc = colno(k)
2629  indx(jc) = ic
2630  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2631  j = colno(jj)
2632  if ( indx(j)==ic ) then
2633  s(1) = s(1) + temp(1,j)*zln(1,jj) + temp(4,j)*zln(4,jj) + temp(7,j)*zln(7,jj)
2634  s(2) = s(2) + temp(2,j)*zln(1,jj) + temp(5,j)*zln(4,jj) + temp(8,j)*zln(7,jj)
2635  s(3) = s(3) + temp(3,j)*zln(1,jj) + temp(6,j)*zln(4,jj) + temp(9,j)*zln(7,jj)
2636 
2637  s(4) = s(4) + temp(1,j)*zln(2,jj) + temp(4,j)*zln(5,jj) + temp(7,j)*zln(8,jj)
2638  s(5) = s(5) + temp(2,j)*zln(2,jj) + temp(5,j)*zln(5,jj) + temp(8,j)*zln(8,jj)
2639  s(6) = s(6) + temp(3,j)*zln(2,jj) + temp(6,j)*zln(5,jj) + temp(9,j)*zln(8,jj)
2640 
2641  s(7) = s(7) + temp(1,j)*zln(3,jj) + temp(4,j)*zln(6,jj) + temp(7,j)*zln(9,jj)
2642  s(8) = s(8) + temp(2,j)*zln(3,jj) + temp(5,j)*zln(6,jj) + temp(8,j)*zln(9,jj)
2643  s(9) = s(9) + temp(3,j)*zln(3,jj) + temp(6,j)*zln(6,jj) + temp(9,j)*zln(9,jj)
2644  endif
2645  enddo
2646  do l = 1, 9
2647  temp(l,jc) = zln(l,k) - s(l)
2648  zln(l,k) = temp(l,jc)
2649  s(l) = 0.0d0
2650  enddo
2651  enddo
2652  end subroutine s3um1
2653 
2654  !======================================================================!
2656  !======================================================================!
2657  subroutine s3um2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
2658  implicit none
2659  !------
2660  integer(kind=kint), intent(in):: Neqns
2661  integer(kind=kint), intent(in):: Nstop
2662  integer(kind=kint), intent(in):: Xlnzr(:)
2663  integer(kind=kint), intent(in):: Colno(:)
2664  integer(kind=kint), intent(inout):: Indx(:)
2665  real(kind=kreal), intent(inout):: diag(6,*)
2666  real(kind=kreal), intent(inout):: dsln(9,*)
2667  real(kind=kreal), intent(inout):: temp(neqns,9)
2668  real(kind=kreal), intent(inout):: zln(9,*)
2669  !------
2670  integer(kind=kint):: ic
2671  integer(kind=kint):: j
2672  integer(kind=kint):: j1
2673  integer(kind=kint):: j2
2674  integer(kind=kint):: jc
2675  integer(kind=kint):: jj
2676  integer(kind=kint):: joc
2677  integer(kind=kint):: k
2678  integer(kind=kint):: ke
2679  integer(kind=kint):: ks
2680 
2681  joc = 0
2682  do ic = nstop, neqns
2683  ks = xlnzr(ic)
2684  ke = xlnzr(ic+1) - 1
2685  do k = ks, ke
2686  jj = colno(k)
2687  temp(jj,1) = zln(1,k)
2688  temp(jj,2) = zln(2,k)
2689  temp(jj,3) = zln(3,k)
2690  temp(jj,4) = zln(4,k)
2691  temp(jj,5) = zln(5,k)
2692  temp(jj,6) = zln(6,k)
2693  temp(jj,7) = zln(7,k)
2694  temp(jj,8) = zln(8,k)
2695  temp(jj,9) = zln(9,k)
2696  indx(jj) = ic
2697  enddo
2698 
2699  do k = ks, ke
2700  jj = colno(k)
2701  zln(4,k) = temp(jj,4) - temp(jj,1)*diag(2,jj)
2702  zln(7,k) = temp(jj,7) - temp(jj,1)*diag(4,jj) - zln(4,k)*diag(5,jj)
2703  zln(1,k) = temp(jj,1)*diag(1,jj)
2704  zln(4,k) = zln(4,k)*diag(3,jj)
2705  zln(7,k) = zln(7,k)*diag(6,jj)
2706  zln(4,k) = zln(4,k) - zln(7,k)*diag(5,jj)
2707  zln(1,k) = zln(1,k) - zln(4,k)*diag(2,jj) - zln(7,k)*diag(4,jj)
2708 
2709  zln(5,k) = temp(jj,5) - temp(jj,2)*diag(2,jj)
2710  zln(8,k) = temp(jj,8) - temp(jj,2)*diag(4,jj) - zln(5,k)*diag(5,jj)
2711  zln(2,k) = temp(jj,2)*diag(1,jj)
2712  zln(5,k) = zln(5,k)*diag(3,jj)
2713  zln(8,k) = zln(8,k)*diag(6,jj)
2714  zln(5,k) = zln(5,k) - zln(8,k)*diag(5,jj)
2715  zln(2,k) = zln(2,k) - zln(5,k)*diag(2,jj) - zln(8,k)*diag(4,jj)
2716 
2717  zln(6,k) = temp(jj,6) - temp(jj,3)*diag(2,jj)
2718  zln(9,k) = temp(jj,9) - temp(jj,3)*diag(4,jj) - zln(6,k)*diag(5,jj)
2719  zln(3,k) = temp(jj,3)*diag(1,jj)
2720  zln(6,k) = zln(6,k)*diag(3,jj)
2721  zln(9,k) = zln(9,k)*diag(6,jj)
2722  zln(6,k) = zln(6,k) - zln(9,k)*diag(5,jj)
2723  zln(3,k) = zln(3,k) - zln(6,k)*diag(2,jj) - zln(9,k)*diag(4,jj)
2724  enddo
2725 
2726  do k = ks, ke
2727  jj = colno(k)
2728  diag(1,ic) = diag(1,ic) - temp(jj,1)*zln(1,k) - temp(jj,4)*zln(4,k) - temp(jj,7)*zln(7,k)
2729  diag(2,ic) = diag(2,ic) - temp(jj,1)*zln(2,k) - temp(jj,4)*zln(5,k) - temp(jj,7)*zln(8,k)
2730  diag(3,ic) = diag(3,ic) - temp(jj,2)*zln(2,k) - temp(jj,5)*zln(5,k) - temp(jj,8)*zln(8,k)
2731  diag(4,ic) = diag(4,ic) - temp(jj,1)*zln(3,k) - temp(jj,4)*zln(6,k) - temp(jj,7)*zln(9,k)
2732  diag(5,ic) = diag(5,ic) - temp(jj,2)*zln(3,k) - temp(jj,5)*zln(6,k) - temp(jj,8)*zln(9,k)
2733  diag(6,ic) = diag(6,ic) - temp(jj,3)*zln(3,k) - temp(jj,6)*zln(6,k) - temp(jj,9)*zln(9,k)
2734  enddo
2735 
2736  do jc = nstop, ic - 1
2737  joc = joc + 1
2738  j1 = xlnzr(jc)
2739  j2 = xlnzr(jc+1)
2740  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2741  j = colno(jj)
2742  if ( indx(j)==ic ) then
2743  dsln(1,joc) = dsln(1,joc) - temp(j,1)*zln(1,jj) - temp(j,4)*zln(4,jj) - temp(j,7)*zln(7,jj)
2744  dsln(2,joc) = dsln(2,joc) - temp(j,2)*zln(1,jj) - temp(j,5)*zln(4,jj) - temp(j,8)*zln(7,jj)
2745  dsln(3,joc) = dsln(3,joc) - temp(j,3)*zln(1,jj) - temp(j,6)*zln(4,jj) - temp(j,9)*zln(7,jj)
2746 
2747  dsln(4,joc) = dsln(4,joc) - temp(j,1)*zln(2,jj) - temp(j,4)*zln(5,jj) - temp(j,7)*zln(8,jj)
2748  dsln(5,joc) = dsln(5,joc) - temp(j,2)*zln(2,jj) - temp(j,5)*zln(5,jj) - temp(j,8)*zln(8,jj)
2749  dsln(6,joc) = dsln(6,joc) - temp(j,3)*zln(2,jj) - temp(j,6)*zln(5,jj) - temp(j,9)*zln(8,jj)
2750 
2751  dsln(7,joc) = dsln(7,joc) - temp(j,1)*zln(3,jj) - temp(j,4)*zln(6,jj) - temp(j,7)*zln(9,jj)
2752  dsln(8,joc) = dsln(8,joc) - temp(j,2)*zln(3,jj) - temp(j,5)*zln(6,jj) - temp(j,8)*zln(9,jj)
2753  dsln(9,joc) = dsln(9,joc) - temp(j,3)*zln(3,jj) - temp(j,6)*zln(6,jj) - temp(j,9)*zln(9,jj)
2754  endif
2755  enddo
2756  enddo
2757  enddo
2758  end subroutine s3um2
2759 
2760  !======================================================================!
2762  !======================================================================!
2763  subroutine s3um3(N,Dsln,Diag,Indx,Temp)
2764  implicit none
2765  !------
2766  integer(kind=kint), intent(in):: N
2767  integer(kind=kint), intent(out):: Indx(:)
2768  real(kind=kreal), intent(inout):: diag(6,*)
2769  real(kind=kreal), intent(inout):: dsln(9,*)
2770  real(kind=kreal), intent(out):: temp(9,*)
2771  !------
2772  integer(kind=kint):: i
2773  integer(kind=kint):: ir
2774  integer(kind=kint):: j
2775  integer(kind=kint):: joc
2776  real(kind=kreal):: t(9)
2777 
2778  if ( n>0 ) then
2779  indx(1) = 0
2780  joc = 1
2781  call inv3(diag(1,1),ir)
2782  do i = 2, n
2783  indx(i) = joc
2784  do j = 1, i - 1
2785  call d3dot(t,dsln(1,indx(i)),dsln(1,indx(j)),j-1)
2786  dsln(:,joc) = dsln(:,joc) - t(:)
2787  joc = joc + 1
2788  enddo
2789  call v3prod(dsln(1,indx(i)),diag,temp,i-1)
2790  call d3dotl(t,temp,dsln(1,indx(i)),i-1)
2791  diag(:,i) = diag(:,i) - t(1:6)
2792  call vcopy(temp,dsln(1,indx(i)),9*(i-1))
2793  call inv3(diag(1,i),ir)
2794  enddo
2795  endif
2796  end subroutine s3um3
2797 
2798  !======================================================================!
2800  !======================================================================!
2801  subroutine s6um(Ic,Xlnzr,Colno,Zln,Diag,Nch,Par,Temp,Indx)
2802  implicit none
2803  !------
2804  integer(kind=kint), intent(in):: Ic
2805  integer(kind=kint), intent(in):: Xlnzr(:)
2806  integer(kind=kint), intent(in):: Colno(:)
2807  integer(kind=kint), intent(in):: Par(:)
2808  integer(kind=kint), intent(inout):: Indx(:)
2809  integer(kind=kint), intent(inout):: Nch(:)
2810  real(kind=kreal), intent(inout):: diag(21,*)
2811  real(kind=kreal), intent(inout):: temp(36,*)
2812  real(kind=kreal), intent(inout):: zln(36,*)
2813  !------
2814  integer(kind=kint):: ir
2815  integer(kind=kint):: j
2816  integer(kind=kint):: jc
2817  integer(kind=kint):: jj
2818  integer(kind=kint):: k
2819  integer(kind=kint):: ke
2820  integer(kind=kint):: kk
2821  integer(kind=kint):: ks
2822  integer(kind=kint):: l
2823  real(kind=kreal):: t(21)
2824  real(kind=kreal):: zz(36)
2825 
2826  ks = xlnzr(ic)
2827  ke = xlnzr(ic+1)
2828  t(1:21) = 0.0d0
2829  do k = ks, ke - 1
2830  jc = colno(k)
2831  indx(jc) = ic
2832  zz(1:36) = zln(1:36,k)
2833  do jj = xlnzr(jc), xlnzr(jc+1) - 1
2834  j = colno(jj)
2835  if ( indx(j)==ic ) then
2836  zz(1) = zz(1) - temp(1,j)*zln(1,jj) - temp(7,j)*zln(7,jj)&
2837  - temp(13,j)*zln(13,jj) - temp(19,j)*zln(19,jj)&
2838  - temp(25,j)*zln(25,jj) - temp(31,j)*zln(31,jj)
2839  zz(2) = zz(2) - temp(2,j)*zln(1,jj) - temp(8,j)*zln(7,jj)&
2840  - temp(14,j)*zln(13,jj) - temp(20,j)*zln(19,jj)&
2841  - temp(26,j)*zln(25,jj) - temp(32,j)*zln(31,jj)
2842  zz(3) = zz(3) - temp(3,j)*zln(1,jj) - temp(9,j)*zln(7,jj)&
2843  - temp(15,j)*zln(13,jj) - temp(21,j)*zln(19,jj)&
2844  - temp(27,j)*zln(25,jj) - temp(33,j)*zln(31,jj)
2845  zz(4) = zz(4) - temp(4,j)*zln(1,jj) - temp(10,j)*zln(7,jj)&
2846  - temp(16,j)*zln(13,jj) - temp(22,j)*zln(19,jj)&
2847  - temp(28,j)*zln(25,jj) - temp(34,j)*zln(31,jj)
2848  zz(5) = zz(5) - temp(5,j)*zln(1,jj) - temp(11,j)*zln(7,jj)&
2849  - temp(17,j)*zln(13,jj) - temp(23,j)*zln(19,jj)&
2850  - temp(29,j)*zln(25,jj) - temp(35,j)*zln(31,jj)
2851  zz(6) = zz(6) - temp(6,j)*zln(1,jj) - temp(12,j)*zln(7,jj)&
2852  - temp(18,j)*zln(13,jj) - temp(24,j)*zln(19,jj)&
2853  - temp(30,j)*zln(25,jj) - temp(36,j)*zln(31,jj)
2854  zz(7) = zz(7) - temp(1,j)*zln(2,jj) - temp(7,j)*zln(8,jj)&
2855  - temp(13,j)*zln(14,jj) - temp(19,j)*zln(20,jj)&
2856  - temp(25,j)*zln(26,jj) - temp(31,j)*zln(32,jj)
2857  zz(8) = zz(8) - temp(2,j)*zln(2,jj) - temp(8,j)*zln(8,jj)&
2858  - temp(14,j)*zln(14,jj) - temp(20,j)*zln(20,jj)&
2859  - temp(26,j)*zln(26,jj) - temp(32,j)*zln(32,jj)
2860  zz(9) = zz(9) - temp(3,j)*zln(2,jj) - temp(9,j)*zln(8,jj)&
2861  - temp(15,j)*zln(14,jj) - temp(21,j)*zln(20,jj)&
2862  - temp(27,j)*zln(26,jj) - temp(33,j)*zln(32,jj)
2863  zz(10) = zz(10) - temp(4,j)*zln(2,jj) - temp(10,j)*zln(8,jj)&
2864  - temp(16,j)*zln(14,jj) - temp(22,j)*zln(20,jj)&
2865  - temp(28,j)*zln(26,jj) - temp(34,j)*zln(32,jj)
2866  zz(11) = zz(11) - temp(5,j)*zln(2,jj) - temp(11,j)*zln(8,jj)&
2867  - temp(17,j)*zln(14,jj) - temp(23,j)*zln(20,jj)&
2868  - temp(29,j)*zln(26,jj) - temp(35,j)*zln(32,jj)
2869  zz(12) = zz(12) - temp(6,j)*zln(2,jj) - temp(12,j)*zln(8,jj)&
2870  - temp(18,j)*zln(14,jj) - temp(24,j)*zln(20,jj)&
2871  - temp(30,j)*zln(26,jj) - temp(36,j)*zln(32,jj)
2872  zz(13) = zz(13) - temp(1,j)*zln(3,jj) - temp(7,j)*zln(9,jj)&
2873  - temp(13,j)*zln(15,jj) - temp(19,j)*zln(21,jj)&
2874  - temp(25,j)*zln(27,jj) - temp(31,j)*zln(33,jj)
2875  zz(14) = zz(14) - temp(2,j)*zln(3,jj) - temp(8,j)*zln(9,jj)&
2876  - temp(14,j)*zln(15,jj) - temp(20,j)*zln(21,jj)&
2877  - temp(26,j)*zln(27,jj) - temp(32,j)*zln(33,jj)
2878  zz(15) = zz(15) - temp(3,j)*zln(3,jj) - temp(9,j)*zln(9,jj)&
2879  - temp(15,j)*zln(15,jj) - temp(21,j)*zln(21,jj)&
2880  - temp(27,j)*zln(27,jj) - temp(33,j)*zln(33,jj)
2881  zz(16) = zz(16) - temp(4,j)*zln(3,jj) - temp(10,j)*zln(9,jj)&
2882  - temp(16,j)*zln(15,jj) - temp(22,j)*zln(21,jj)&
2883  - temp(28,j)*zln(27,jj) - temp(34,j)*zln(33,jj)
2884  zz(17) = zz(17) - temp(5,j)*zln(3,jj) - temp(11,j)*zln(9,jj)&
2885  - temp(17,j)*zln(15,jj) - temp(23,j)*zln(21,jj)&
2886  - temp(29,j)*zln(27,jj) - temp(35,j)*zln(33,jj)
2887  zz(18) = zz(18) - temp(6,j)*zln(3,jj) - temp(12,j)*zln(9,jj)&
2888  - temp(18,j)*zln(15,jj) - temp(24,j)*zln(21,jj)&
2889  - temp(30,j)*zln(27,jj) - temp(36,j)*zln(33,jj)
2890  zz(19) = zz(19) - temp(1,j)*zln(4,jj) - temp(7,j)*zln(10,jj)&
2891  - temp(13,j)*zln(16,jj) - temp(19,j)*zln(22,jj)&
2892  - temp(25,j)*zln(28,jj) - temp(31,j)*zln(34,jj)
2893  zz(20) = zz(20) - temp(2,j)*zln(4,jj) - temp(8,j)*zln(10,jj)&
2894  - temp(14,j)*zln(16,jj) - temp(20,j)*zln(22,jj)&
2895  - temp(26,j)*zln(28,jj) - temp(32,j)*zln(34,jj)
2896  zz(21) = zz(21) - temp(3,j)*zln(4,jj) - temp(9,j)*zln(10,jj)&
2897  - temp(15,j)*zln(16,jj) - temp(21,j)*zln(22,jj)&
2898  - temp(27,j)*zln(28,jj) - temp(33,j)*zln(34,jj)
2899  zz(22) = zz(22) - temp(4,j)*zln(4,jj) - temp(10,j)*zln(10,jj)&
2900  - temp(16,j)*zln(16,jj) - temp(22,j)*zln(22,jj)&
2901  - temp(28,j)*zln(28,jj) - temp(34,j)*zln(34,jj)
2902  zz(23) = zz(23) - temp(5,j)*zln(4,jj) - temp(11,j)*zln(10,jj)&
2903  - temp(17,j)*zln(16,jj) - temp(23,j)*zln(22,jj)&
2904  - temp(29,j)*zln(28,jj) - temp(35,j)*zln(34,jj)
2905  zz(24) = zz(24) - temp(6,j)*zln(4,jj) - temp(12,j)*zln(10,jj)&
2906  - temp(18,j)*zln(16,jj) - temp(24,j)*zln(22,jj)&
2907  - temp(30,j)*zln(28,jj) - temp(36,j)*zln(34,jj)
2908  zz(25) = zz(25) - temp(1,j)*zln(5,jj) - temp(7,j)*zln(11,jj)&
2909  - temp(13,j)*zln(17,jj) - temp(19,j)*zln(23,jj)&
2910  - temp(25,j)*zln(29,jj) - temp(31,j)*zln(35,jj)
2911  zz(26) = zz(26) - temp(2,j)*zln(5,jj) - temp(8,j)*zln(11,jj)&
2912  - temp(14,j)*zln(17,jj) - temp(20,j)*zln(23,jj)&
2913  - temp(26,j)*zln(29,jj) - temp(32,j)*zln(35,jj)
2914  zz(27) = zz(27) - temp(3,j)*zln(5,jj) - temp(9,j)*zln(11,jj)&
2915  - temp(15,j)*zln(17,jj) - temp(21,j)*zln(23,jj)&
2916  - temp(27,j)*zln(29,jj) - temp(33,j)*zln(35,jj)
2917  zz(28) = zz(28) - temp(4,j)*zln(5,jj) - temp(10,j)*zln(11,jj)&
2918  - temp(16,j)*zln(17,jj) - temp(22,j)*zln(23,jj)&
2919  - temp(28,j)*zln(29,jj) - temp(34,j)*zln(35,jj)
2920  zz(29) = zz(29) - temp(5,j)*zln(5,jj) - temp(11,j)*zln(11,jj)&
2921  - temp(17,j)*zln(17,jj) - temp(23,j)*zln(23,jj)&
2922  - temp(29,j)*zln(29,jj) - temp(35,j)*zln(35,jj)
2923  zz(30) = zz(30) - temp(6,j)*zln(5,jj) - temp(12,j)*zln(11,jj)&
2924  - temp(18,j)*zln(17,jj) - temp(24,j)*zln(23,jj)&
2925  - temp(30,j)*zln(29,jj) - temp(36,j)*zln(35,jj)
2926  zz(31) = zz(31) - temp(1,j)*zln(6,jj) - temp(7,j)*zln(12,jj)&
2927  - temp(13,j)*zln(18,jj) - temp(19,j)*zln(24,jj)&
2928  - temp(25,j)*zln(30,jj) - temp(31,j)*zln(36,jj)
2929  zz(32) = zz(32) - temp(2,j)*zln(6,jj) - temp(8,j)*zln(12,jj)&
2930  - temp(14,j)*zln(18,jj) - temp(20,j)*zln(24,jj)&
2931  - temp(26,j)*zln(30,jj) - temp(32,j)*zln(36,jj)
2932  zz(33) = zz(33) - temp(3,j)*zln(6,jj) - temp(9,j)*zln(12,jj)&
2933  - temp(15,j)*zln(18,jj) - temp(21,j)*zln(24,jj)&
2934  - temp(27,j)*zln(30,jj) - temp(33,j)*zln(36,jj)
2935  zz(34) = zz(34) - temp(4,j)*zln(6,jj) - temp(10,j)*zln(12,jj)&
2936  - temp(16,j)*zln(18,jj) - temp(22,j)*zln(24,jj)&
2937  - temp(28,j)*zln(30,jj) - temp(34,j)*zln(36,jj)
2938  zz(35) = zz(35) - temp(5,j)*zln(6,jj) - temp(11,j)*zln(12,jj)&
2939  - temp(17,j)*zln(18,jj) - temp(23,j)*zln(24,jj)&
2940  - temp(29,j)*zln(30,jj) - temp(35,j)*zln(36,jj)
2941  zz(36) = zz(36) - temp(6,j)*zln(6,jj) - temp(12,j)*zln(12,jj)&
2942  - temp(18,j)*zln(18,jj) - temp(24,j)*zln(24,jj)&
2943  - temp(30,j)*zln(30,jj) - temp(36,j)*zln(36,jj)
2944  endif
2945  enddo
2946  call inv66(zln(1,k),zz,diag(1,jc))
2947  temp(1:36,jc) = zz(1:36)
2948 
2949  t(1) = t(1) + zz(1)*zln(1,k) + zz(7)*zln(7,k) + zz(13)*zln(13,k)&
2950  + zz(19)*zln(19,k) + zz(25)*zln(25,k) + zz(31)*zln(31,k)
2951  t(2) = t(2) + zz(1)*zln(2,k) + zz(7)*zln(8,k) + zz(13)*zln(14,k)&
2952  + zz(19)*zln(20,k) + zz(25)*zln(26,k) + zz(31)*zln(32,k)
2953  t(3) = t(3) + zz(2)*zln(2,k) + zz(8)*zln(8,k) + zz(14)*zln(14,k)&
2954  + zz(20)*zln(20,k) + zz(26)*zln(26,k) + zz(32)*zln(32,k)
2955  t(4) = t(4) + zz(1)*zln(3,k) + zz(7)*zln(9,k) + zz(13)*zln(15,k)&
2956  + zz(19)*zln(21,k) + zz(25)*zln(27,k) + zz(31)*zln(33,k)
2957  t(5) = t(5) + zz(2)*zln(3,k) + zz(8)*zln(9,k) + zz(14)*zln(15,k)&
2958  + zz(20)*zln(21,k) + zz(26)*zln(27,k) + zz(32)*zln(33,k)
2959  t(6) = t(6) + zz(3)*zln(3,k) + zz(9)*zln(9,k) + zz(15)*zln(15,k)&
2960  + zz(21)*zln(21,k) + zz(27)*zln(27,k) + zz(33)*zln(33,k)
2961  t(7) = t(7) + zz(1)*zln(4,k) + zz(7)*zln(10,k) + zz(13)*zln(16,k)&
2962  + zz(19)*zln(22,k) + zz(25)*zln(28,k) + zz(31)*zln(34,k)
2963  t(8) = t(8) + zz(2)*zln(4,k) + zz(8)*zln(10,k) + zz(14)*zln(16,k)&
2964  + zz(20)*zln(22,k) + zz(26)*zln(28,k) + zz(32)*zln(34,k)
2965  t(9) = t(9) + zz(3)*zln(4,k) + zz(9)*zln(10,k) + zz(15)*zln(16,k)&
2966  + zz(21)*zln(22,k) + zz(27)*zln(28,k) + zz(33)*zln(34,k)
2967  t(10) = t(10) + zz(4)*zln(4,k) + zz(10)*zln(10,k) + zz(16)*zln(16,k)&
2968  + zz(22)*zln(22,k) + zz(28)*zln(28,k) + zz(34)*zln(34,k)
2969  t(11) = t(11) + zz(1)*zln(5,k) + zz(7)*zln(11,k) + zz(13)*zln(17,k)&
2970  + zz(19)*zln(23,k) + zz(25)*zln(29,k) + zz(31)*zln(35,k)
2971  t(12) = t(12) + zz(2)*zln(5,k) + zz(8)*zln(11,k) + zz(14)*zln(17,k)&
2972  + zz(20)*zln(23,k) + zz(26)*zln(29,k) + zz(32)*zln(35,k)
2973  t(13) = t(13) + zz(3)*zln(5,k) + zz(9)*zln(11,k) + zz(15)*zln(17,k)&
2974  + zz(21)*zln(23,k) + zz(27)*zln(29,k) + zz(33)*zln(35,k)
2975  t(14) = t(14) + zz(4)*zln(5,k) + zz(10)*zln(11,k) + zz(16)*zln(17,k)&
2976  + zz(22)*zln(23,k) + zz(28)*zln(29,k) + zz(34)*zln(35,k)
2977  t(15) = t(15) + zz(5)*zln(5,k) + zz(11)*zln(11,k) + zz(17)*zln(17,k)&
2978  + zz(23)*zln(23,k) + zz(29)*zln(29,k) + zz(35)*zln(35,k)
2979  t(16) = t(16) + zz(1)*zln(6,k) + zz(7)*zln(12,k) + zz(13)*zln(18,k)&
2980  + zz(19)*zln(24,k) + zz(25)*zln(30,k) + zz(31)*zln(36,k)
2981  t(17) = t(17) + zz(2)*zln(6,k) + zz(8)*zln(12,k) + zz(14)*zln(18,k)&
2982  + zz(20)*zln(24,k) + zz(26)*zln(30,k) + zz(32)*zln(36,k)
2983  t(18) = t(18) + zz(3)*zln(6,k) + zz(9)*zln(12,k) + zz(15)*zln(18,k)&
2984  + zz(21)*zln(24,k) + zz(27)*zln(30,k) + zz(33)*zln(36,k)
2985  t(19) = t(19) + zz(4)*zln(6,k) + zz(10)*zln(12,k) + zz(16)*zln(18,k)&
2986  + zz(22)*zln(24,k) + zz(28)*zln(30,k) + zz(34)*zln(36,k)
2987  t(20) = t(20) + zz(5)*zln(6,k) + zz(11)*zln(12,k) + zz(17)*zln(18,k)&
2988  + zz(23)*zln(24,k) + zz(29)*zln(30,k) + zz(35)*zln(36,k)
2989  t(21) = t(21) + zz(6)*zln(6,k) + zz(12)*zln(12,k) + zz(18)*zln(18,k)&
2990  + zz(24)*zln(24,k) + zz(30)*zln(30,k) + zz(36)*zln(36,k)
2991  enddo
2992  do l = 1, 21
2993  diag(l,ic) = diag(l,ic) - t(l)
2994  enddo
2995  call inv6(diag(1,ic),ir)
2996  nch(ic) = -1
2997  kk = par(ic)
2998  nch(kk) = nch(kk) - 1
2999  end subroutine s6um
3000 
3001  !======================================================================!
3003  !======================================================================!
3004  subroutine s6um1(Ic,Xlnzr,Colno,Zln,Temp,Indx)
3005  implicit none
3006  !------
3007  integer(kind=kint), intent(in):: Ic
3008  integer(kind=kint), intent(in):: Xlnzr(:)
3009  integer(kind=kint), intent(in):: Colno(:)
3010  integer(kind=kint), intent(inout):: Indx(:)
3011  real(kind=kreal), intent(inout):: temp(9,*)
3012  real(kind=kreal), intent(inout):: zln(9,*)
3013  !------
3014  integer(kind=kint):: j
3015  integer(kind=kint):: jc
3016  integer(kind=kint):: jj
3017  integer(kind=kint):: k
3018  integer(kind=kint):: ke
3019  integer(kind=kint):: ks
3020  integer(kind=kint):: l
3021  real(kind=kreal):: s(9)
3022 
3023  ks = xlnzr(ic)
3024  ke = xlnzr(ic+1)
3025  s(1:9) = 0.0d0
3026  do k = ks, ke - 1
3027  jc = colno(k)
3028  indx(jc) = ic
3029  do jj = xlnzr(jc), xlnzr(jc+1) - 1
3030  j = colno(jj)
3031  if ( indx(j)==ic ) then
3032  s(1) = s(1) + temp(1,j)*zln(1,jj) + temp(4,j)*zln(4,jj) + temp(7,j)*zln(7,jj)
3033  s(2) = s(2) + temp(2,j)*zln(1,jj) + temp(5,j)*zln(4,jj) + temp(8,j)*zln(7,jj)
3034  s(3) = s(3) + temp(3,j)*zln(1,jj) + temp(6,j)*zln(4,jj) + temp(9,j)*zln(7,jj)
3035 
3036  s(4) = s(4) + temp(1,j)*zln(2,jj) + temp(4,j)*zln(5,jj) + temp(7,j)*zln(8,jj)
3037  s(5) = s(5) + temp(2,j)*zln(2,jj) + temp(5,j)*zln(5,jj) + temp(8,j)*zln(8,jj)
3038  s(6) = s(6) + temp(3,j)*zln(2,jj) + temp(6,j)*zln(5,jj) + temp(9,j)*zln(8,jj)
3039 
3040  s(7) = s(7) + temp(1,j)*zln(3,jj) + temp(4,j)*zln(6,jj) + temp(7,j)*zln(9,jj)
3041  s(8) = s(8) + temp(2,j)*zln(3,jj) + temp(5,j)*zln(6,jj) + temp(8,j)*zln(9,jj)
3042  s(9) = s(9) + temp(3,j)*zln(3,jj) + temp(6,j)*zln(6,jj) + temp(9,j)*zln(9,jj)
3043  endif
3044  enddo
3045  do l = 1, 9
3046  temp(l,jc) = zln(l,k) - s(l)
3047  zln(l,k) = temp(l,jc)
3048  s(l) = 0.0d0
3049  enddo
3050  enddo
3051  end subroutine s6um1
3052 
3053  !======================================================================!
3055  !======================================================================!
3056  subroutine s6um2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx)
3057  implicit none
3058  !------
3059  integer(kind=kint), intent(in):: Neqns
3060  integer(kind=kint), intent(in):: Nstop
3061  integer(kind=kint), intent(in):: Xlnzr(:)
3062  integer(kind=kint), intent(in):: Colno(:)
3063  integer(kind=kint), intent(inout):: Indx(:)
3064  real(kind=kreal), intent(inout):: diag(21,*)
3065  real(kind=kreal), intent(inout):: dsln(36,*)
3066  real(kind=kreal), intent(inout):: temp(36,neqns)
3067  real(kind=kreal), intent(inout):: zln(36,*)
3068  !------
3069  integer(kind=kint):: ic
3070  integer(kind=kint):: j
3071  integer(kind=kint):: j1
3072  integer(kind=kint):: j2
3073  integer(kind=kint):: jc
3074  integer(kind=kint):: jj
3075  integer(kind=kint):: joc
3076  integer(kind=kint):: k
3077  integer(kind=kint):: ke
3078  integer(kind=kint):: ks
3079 
3080  joc = 0
3081  do ic = nstop, neqns
3082  temp(1:nstop,1:36) = 0.0d0
3083  ks = xlnzr(ic)
3084  ke = xlnzr(ic+1) - 1
3085  do k = ks, ke
3086  jj = colno(k)
3087  temp(:,jj) = zln(:,k)
3088  indx(jj) = ic
3089  enddo
3090  do k = ks, ke
3091  jj = colno(k)
3092  call inv66(zln(1,k),temp,diag(1,jj))
3093  enddo
3094 
3095  do k = ks, ke
3096  jj = colno(k)
3097  diag(1,ic) = diag(1,ic) - temp(jj,1)*zln(1,k) - temp(jj,4)*zln(4,k) - temp(jj,7)*zln(7,k)
3098  diag(2,ic) = diag(2,ic) - temp(jj,1)*zln(2,k) - temp(jj,4)*zln(5,k) - temp(jj,7)*zln(8,k)
3099  diag(3,ic) = diag(3,ic) - temp(jj,2)*zln(2,k) - temp(jj,5)*zln(5,k) - temp(jj,8)*zln(8,k)
3100  diag(4,ic) = diag(4,ic) - temp(jj,1)*zln(3,k) - temp(jj,4)*zln(6,k) - temp(jj,7)*zln(9,k)
3101  diag(5,ic) = diag(5,ic) - temp(jj,2)*zln(3,k) - temp(jj,5)*zln(6,k) - temp(jj,8)*zln(9,k)
3102  diag(6,ic) = diag(6,ic) - temp(jj,3)*zln(3,k) - temp(jj,6)*zln(6,k) - temp(jj,9)*zln(9,k)
3103  enddo
3104  do jc = nstop, ic - 1
3105  joc = joc + 1
3106  j1 = xlnzr(jc)
3107  j2 = xlnzr(jc+1)
3108  do jj = xlnzr(jc), xlnzr(jc+1) - 1
3109  j = colno(jj)
3110  if ( indx(j)==ic ) then
3111  dsln(1,joc) = dsln(1,joc) - temp(j,1)*zln(1,jj) - temp(j,4)*zln(4,jj) - temp(j,7)*zln(7,jj)
3112  dsln(2,joc) = dsln(2,joc) - temp(j,2)*zln(1,jj) - temp(j,5)*zln(4,jj) - temp(j,8)*zln(7,jj)
3113  dsln(3,joc) = dsln(3,joc) - temp(j,3)*zln(1,jj) - temp(j,6)*zln(4,jj) - temp(j,9)*zln(7,jj)
3114 
3115  dsln(4,joc) = dsln(4,joc) - temp(j,1)*zln(2,jj) - temp(j,4)*zln(5,jj) - temp(j,7)*zln(8,jj)
3116  dsln(5,joc) = dsln(5,joc) - temp(j,2)*zln(2,jj) - temp(j,5)*zln(5,jj) - temp(j,8)*zln(8,jj)
3117  dsln(6,joc) = dsln(6,joc) - temp(j,3)*zln(2,jj) - temp(j,6)*zln(5,jj) - temp(j,9)*zln(8,jj)
3118 
3119  dsln(7,joc) = dsln(7,joc) - temp(j,1)*zln(3,jj) - temp(j,4)*zln(6,jj) - temp(j,7)*zln(9,jj)
3120  dsln(8,joc) = dsln(8,joc) - temp(j,2)*zln(3,jj) - temp(j,5)*zln(6,jj) - temp(j,8)*zln(9,jj)
3121  dsln(9,joc) = dsln(9,joc) - temp(j,3)*zln(3,jj) - temp(j,6)*zln(6,jj) - temp(j,9)*zln(9,jj)
3122  endif
3123  enddo
3124  enddo
3125  enddo
3126  end subroutine s6um2
3127 
3128  !======================================================================!
3130  !======================================================================!
3131  subroutine s6um3(N,Dsln,Diag,Indx,Temp)
3132  implicit none
3133  !------
3134  integer(kind=kint), intent(in):: N
3135  integer(kind=kint), intent(out):: Indx(:)
3136  real(kind=kreal), intent(inout):: diag(6,*)
3137  real(kind=kreal), intent(inout):: dsln(9,*)
3138  real(kind=kreal), intent(out):: temp(9,*)
3139  !------
3140  integer(kind=kint):: i
3141  integer(kind=kint):: ir
3142  integer(kind=kint):: j
3143  integer(kind=kint):: joc
3144  integer(kind=kint):: l
3145  real(kind=kreal):: t(9)
3146 
3147  if ( n>0 ) then
3148  indx(1) = 0
3149  joc = 1
3150  call inv3(diag(1,1),ir)
3151  do i = 2, n
3152  indx(i) = joc
3153  do j = 1, i - 1
3154  call d3dot(t,dsln(1,indx(i)),dsln(1,indx(j)),j-1)
3155  do l = 1, 9
3156  dsln(l,joc) = dsln(l,joc) - t(l)
3157  enddo
3158  joc = joc + 1
3159  enddo
3160  call v3prod(dsln(1,indx(i)),diag,temp,i-1)
3161  call d3dotl(t,temp,dsln(1,indx(i)),i-1)
3162  do l = 1, 6
3163  diag(l,i) = diag(l,i) - t(l)
3164  enddo
3165  call vcopy(temp,dsln(1,indx(i)),9*(i-1))
3166  call inv3(diag(1,i),ir)
3167  enddo
3168  endif
3169  end subroutine s6um3
3170 
3171  !======================================================================!
3173  !======================================================================!
3174  subroutine sxum(Ic,Xlnzr,Colno,Zln,Diag,Nch,Par,Temp,Indx,Ndeg,Ndegl,Zz,T)
3175  implicit none
3176  !------
3177  integer(kind=kint), intent(in):: Ic
3178  integer(kind=kint), intent(in):: Ndeg
3179  integer(kind=kint), intent(in):: Ndegl
3180  integer(kind=kint), intent(in):: Xlnzr(:)
3181  integer(kind=kint), intent(in):: Colno(:)
3182  integer(kind=kint), intent(in):: Par(:)
3183  integer(kind=kint), intent(inout):: Indx(:)
3184  integer(kind=kint), intent(inout):: Nch(:)
3185  real(kind=kreal), intent(inout):: diag(ndegl,*)
3186  real(kind=kreal), intent(out):: t(ndegl)
3187  real(kind=kreal), intent(inout):: temp(ndeg,ndeg,*)
3188  real(kind=kreal), intent(inout):: zln(ndeg,ndeg,*)
3189  real(kind=kreal), intent(out):: zz(ndeg,ndeg)
3190  !------
3191  integer(kind=kint):: ir
3192  integer(kind=kint):: j
3193  integer(kind=kint):: jc
3194  integer(kind=kint):: jj
3195  integer(kind=kint):: joc
3196  integer(kind=kint):: k
3197  integer(kind=kint):: ke
3198  integer(kind=kint):: kk
3199  integer(kind=kint):: ks
3200  integer(kind=kint):: m
3201  integer(kind=kint):: n
3202  integer(kind=kint):: ndeg22
3203 
3204  ndeg22 = ndeg*ndeg
3205  ks = xlnzr(ic)
3206  ke = xlnzr(ic+1)
3207  t = 0.0
3208  do k = ks, ke - 1
3209  jc = colno(k)
3210  indx(jc) = ic
3211  zz = zln(:,:,k)
3212  do jj = xlnzr(jc), xlnzr(jc+1) - 1
3213  j = colno(jj)
3214  if ( indx(j)==ic ) then
3215  do m = 1, ndeg, 2
3216  do n = 1, ndeg, 2
3217  do kk = 1, ndeg, 2
3218  zz(n,m) = zz(n,m) - temp(n,kk,j)*zln(m,kk,jj) - temp(n,kk+1,j)*zln(m,kk+1,jj)
3219  zz(n,m+1) = zz(n,m+1) - temp(n,kk,j)*zln(m+1,kk,jj) - temp(n,kk+1,j)*zln(m+1,kk+1,jj)
3220  zz(n+1,m) = zz(n+1,m) - temp(n+1,kk,j)*zln(m,kk,jj) - temp(n+1,kk+1,j)*zln(m,kk+1,jj)
3221  zz(n+1,m+1) = zz(n+1,m+1) - temp(n+1,kk,j)*zln(m+1,kk,jj) - temp(n+1,kk+1,j)*zln(m+1,kk+1,jj)
3222  enddo
3223  enddo
3224  enddo
3225  endif
3226  enddo
3227  call invxx(zln(1,1,k),zz,diag(1,jc),ndeg)
3228 
3229  temp(:,:,jc) = zz
3230  joc = 0
3231  do n = 1, ndeg
3232  do m = 1, n
3233  joc = joc + 1
3234  do kk = 1, ndeg, 2
3235  t(joc) = t(joc) + zz(n,kk)*zln(m,kk,k) + zz(n,kk+1)*zln(m,kk+1,k)
3236  enddo
3237  enddo
3238  enddo
3239  enddo
3240 
3241  diag(:,ic) = diag(:,ic) - t
3242  call invx(diag(1,ic),ndeg,ir)
3243  nch(ic) = -1
3244  kk = par(ic)
3245  nch(kk) = nch(kk) - 1
3246  end subroutine sxum
3247 
3248  !======================================================================!
3250  !======================================================================!
3251  subroutine sxum1(Ic,Xlnzr,Colno,Zln,Temp,Indx,Ndeg,S)
3252  implicit none
3253  !------
3254  integer(kind=kint), intent(in):: Ndeg
3255  integer(kind=kint), intent(in):: Xlnzr(:)
3256  integer(kind=kint), intent(in):: Colno(:)
3257  integer(kind=kint), intent(inout):: Indx(:)
3258  real(kind=kreal), intent(inout):: s(ndeg,ndeg)
3259  real(kind=kreal), intent(inout):: temp(ndeg,ndeg,*)
3260  real(kind=kreal), intent(inout):: zln(ndeg,ndeg,*)
3261  !------
3262  integer(kind=kint):: Ic
3263  integer(kind=kint):: j
3264  integer(kind=kint):: jc
3265  integer(kind=kint):: jj
3266  integer(kind=kint):: k
3267  integer(kind=kint):: ke
3268  integer(kind=kint):: kk
3269  integer(kind=kint):: ks
3270  integer(kind=kint):: m
3271  integer(kind=kint):: n
3272 
3273  ks = xlnzr(ic)
3274  ke = xlnzr(ic+1)
3275  s(1:ndeg,1:ndeg) = 0.0d0
3276 
3277  do k = ks, ke - 1
3278  jc = colno(k)
3279  indx(jc) = ic
3280  do jj = xlnzr(jc), xlnzr(jc+1) - 1
3281  j = colno(jj)
3282  if ( indx(j)==ic ) then
3283  do m = 1, ndeg
3284  do n = 1, ndeg
3285  do kk = 1, ndeg
3286  s(n,m) = s(n,m) + temp(n,kk,j)*zln(m,kk,jj)
3287  enddo
3288  enddo
3289  enddo
3290  endif
3291  enddo
3292  do m = 1, ndeg
3293  do n = 1, ndeg
3294  temp(n,m,jc) = zln(n,m,k) - s(n,m)
3295  zln(n,m,k) = temp(n,m,jc)
3296  s(n,m) = 0.0d0
3297  enddo
3298  enddo
3299  enddo
3300  end subroutine sxum1
3301 
3302  !======================================================================!
3304  !======================================================================!
3305  subroutine sxum2(Neqns,Nstop,Xlnzr,Colno,Zln,Diag,Dsln,Temp,Indx,Ndeg,Ndegl)
3306  implicit none
3307  !------
3308  integer(kind=kint), intent(in):: Ndeg
3309  integer(kind=kint), intent(in):: Ndegl
3310  integer(kind=kint), intent(in):: Neqns
3311  integer(kind=kint), intent(in):: Nstop
3312  integer(kind=kint), intent(in):: Xlnzr(:)
3313  integer(kind=kint), intent(in):: Colno(:)
3314  integer(kind=kint), intent(inout):: Indx(:)
3315  real(kind=kreal), intent(inout):: diag(ndegl,*)
3316  real(kind=kreal), intent(inout):: dsln(ndeg,ndeg,*)
3317  real(kind=kreal), intent(inout):: temp(ndeg,ndeg,*)
3318  real(kind=kreal), intent(inout):: zln(ndeg,ndeg,*)
3319  !------
3320  integer(kind=kint):: ic
3321  integer(kind=kint):: j
3322  integer(kind=kint):: j1
3323  integer(kind=kint):: j2
3324  integer(kind=kint):: jc
3325  integer(kind=kint):: jj
3326  integer(kind=kint):: joc
3327  integer(kind=kint):: k
3328  integer(kind=kint):: ke
3329  integer(kind=kint):: kk
3330  integer(kind=kint):: ks
3331  integer(kind=kint):: locd
3332  integer(kind=kint):: m
3333  integer(kind=kint):: n
3334 
3335  joc = 0
3336  do ic = nstop, neqns
3337  ks = xlnzr(ic)
3338  ke = xlnzr(ic+1) - 1
3339  do k = ks, ke
3340  jj = colno(k)
3341  do m = 1, ndeg
3342  temp(1:ndeg,m,jj) = zln(1:ndeg,m,k)
3343  indx(jj) = ic
3344  enddo
3345  enddo
3346  do k = ks, ke
3347  jj = colno(k)
3348  call invxx(zln(1,1,k),temp(1,1,jj),diag(1,jj),ndeg)
3349  enddo
3350 
3351  locd = 0
3352  do n = 1, ndeg
3353  do m = 1, n
3354  locd = locd + 1
3355  do k = ks, ke
3356  jj = colno(k)
3357  do kk = 1, ndeg
3358  diag(locd,ic) = diag(locd,ic) - temp(n,kk,jj)*zln(m,kk,k)
3359  enddo
3360  enddo
3361  enddo
3362  enddo
3363  do jc = nstop, ic - 1
3364  joc = joc + 1
3365  j1 = xlnzr(jc)
3366  j2 = xlnzr(jc+1)
3367  do jj = xlnzr(jc), xlnzr(jc+1) - 1
3368  j = colno(jj)
3369  if ( indx(j)==ic ) then
3370  do m = 1, ndeg
3371  do n = 1, ndeg
3372  do k = 1, ndeg
3373  dsln(n,m,joc) = dsln(n,m,joc) - temp(n,k,j)*zln(m,k,jj)
3374  enddo
3375  enddo
3376  enddo
3377  endif
3378  enddo
3379  enddo
3380  enddo
3381  end subroutine sxum2
3382 
3383  !======================================================================!
3385  !======================================================================!
3386  subroutine sxum3(Nn,Dsln,Diag,Indx,Temp,Ndeg,Ndegl,T)
3387  implicit none
3388  !------
3389  integer(kind=kint), intent(in):: Ndeg
3390  integer(kind=kint), intent(in):: Ndegl
3391  integer(kind=kint), intent(in):: Nn
3392  integer(kind=kint), intent(out):: Indx(:)
3393  real(kind=kreal), intent(inout):: diag(ndegl,*)
3394  real(kind=kreal), intent(inout):: dsln(ndeg,ndeg,*)
3395  real(kind=kreal), intent(out):: t(ndeg,ndeg)
3396  real(kind=kreal), intent(out):: temp(ndeg,ndeg,*)
3397  !------
3398  integer(kind=kint):: i
3399  integer(kind=kint):: ir
3400  integer(kind=kint):: j
3401  integer(kind=kint):: joc
3402  integer(kind=kint):: locd
3403  integer(kind=kint):: m
3404  integer(kind=kint):: n
3405 
3406  if ( nn>0 ) then
3407  indx(1) = 0
3408  joc = 1
3409  call invx(diag(1,1),ndeg,ir)
3410  do i = 2, nn
3411  indx(i) = joc
3412  do j = 1, i - 1
3413  call dxdot(ndeg,t,dsln(1,1,indx(i)),dsln(1,1,indx(j)),j-1)
3414  do m = 1, ndeg
3415  dsln(1:ndeg,m,joc) = dsln(1:ndeg,m,joc) - t(1:ndeg,m)
3416  enddo
3417  joc = joc + 1
3418  enddo
3419  call vxprod(ndeg,ndegl,dsln(1,1,indx(i)),diag,temp,i-1)
3420  call dxdotl(ndeg,t,temp,dsln(1,1,indx(i)),i-1)
3421  locd = 0
3422  do n = 1, ndeg
3423  do m = 1, n
3424  locd = locd + 1
3425  diag(locd,i) = diag(locd,i) - t(n,m)
3426  enddo
3427  enddo
3428  call vcopy(temp,dsln(1,1,indx(i)),ndeg*ndeg*(i-1))
3429  call invx(diag(1,i),ndeg,ir)
3430  enddo
3431  endif
3432  end subroutine sxum3
3433 
3434  !======================================================================!
3436  !======================================================================!
3437  subroutine inv2(Dsln,Ir)
3438  implicit none
3439  !------
3440  integer(kind=kint), intent(out):: Ir
3441  real(kind=kreal), intent(inout):: dsln(3)
3442  !------
3443  real(kind=kreal):: t
3444 
3445  ir = 0
3446  if ( dabs(dsln(1))<rmin ) then
3447  ir = 10
3448  return
3449  endif
3450  dsln(1) = 1.0d0/dsln(1)
3451  t = dsln(2)*dsln(1)
3452  dsln(3) = dsln(3) - t*dsln(2)
3453  dsln(2) = t
3454  if ( dabs(dsln(3))<rmin ) then
3455  ir = 10
3456  return
3457  endif
3458  dsln(3) = 1.0d0/dsln(3)
3459  end subroutine inv2
3460 
3461  !======================================================================!
3463  !======================================================================!
3464  subroutine inv22(Zln,Zz,Diag)
3465  implicit none
3466  !------
3467  real(kind=kreal), intent(in):: diag(3)
3468  real(kind=kreal), intent(in):: zz(4)
3469  real(kind=kreal), intent(out):: zln(4)
3470  !------
3471  zln(3) = zz(3) - zz(1)*diag(2)
3472  zln(1) = zz(1)*diag(1)
3473  zln(3) = zln(3)*diag(3)
3474  zln(1) = zln(1) - zln(3)*diag(2)
3475 
3476  zln(4) = zz(4) - zz(2)*diag(2)
3477  zln(2) = zz(2)*diag(1)
3478  zln(4) = zln(4)*diag(3)
3479  zln(2) = zln(2) - zln(4)*diag(2)
3480  end subroutine inv22
3481 
3482  !======================================================================!
3484  !======================================================================!
3485  subroutine inv3(Dsln,Ir)
3486  implicit none
3487  !------
3488  integer(kind=kint), intent(out):: Ir
3489  real(kind=kreal), intent(inout):: dsln(6)
3490  !------
3491  real(kind=kreal):: t(2)
3492 
3493  ir = 0
3494  do
3495  if ( dabs(dsln(1))<rmin ) exit
3496  dsln(1) = 1.0d0/dsln(1)
3497  t(1) = dsln(2)*dsln(1)
3498  dsln(3) = dsln(3) - t(1)*dsln(2)
3499  dsln(2) = t(1)
3500  if ( dabs(dsln(3))<rmin ) exit
3501  dsln(3) = 1.0d0/dsln(3)
3502  t(1) = dsln(4)*dsln(1)
3503  dsln(5) = dsln(5) - dsln(2)*dsln(4)
3504  t(2) = dsln(5)*dsln(3)
3505  dsln(6) = dsln(6) - t(1)*dsln(4) - t(2)*dsln(5)
3506  dsln(4) = t(1)
3507  dsln(5) = t(2)
3508  if ( dabs(dsln(6))<rmin ) exit
3509  dsln(6) = 1.0d0/dsln(6)
3510  return
3511  enddo
3512 
3513  dsln(1) = 1.0d0
3514  dsln(2) = 0.0d0
3515  dsln(3) = 1.0d0
3516  dsln(4) = 0.0d0
3517  dsln(5) = 0.0d0
3518  dsln(6) = 1.0d0
3519  end subroutine inv3
3520 
3521  !======================================================================!
3523  !======================================================================!
3524  subroutine inv33(Zln,Zz,Diag)
3525  implicit none
3526  !------
3527  real(kind=kreal), intent(in):: diag(6)
3528  real(kind=kreal), intent(in):: zz(9)
3529  real(kind=kreal), intent(out):: zln(9)
3530  !------
3531  zln(4) = zz(4) - zz(1)*diag(2)
3532  zln(7) = zz(7) - zz(1)*diag(4) - zln(4)*diag(5)
3533  zln(1) = zz(1)*diag(1)
3534  zln(4) = zln(4)*diag(3)
3535  zln(7) = zln(7)*diag(6)
3536  zln(4) = zln(4) - zln(7)*diag(5)
3537  zln(1) = zln(1) - zln(4)*diag(2) - zln(7)*diag(4)
3538 
3539  zln(5) = zz(5) - zz(2)*diag(2)
3540  zln(8) = zz(8) - zz(2)*diag(4) - zln(5)*diag(5)
3541  zln(2) = zz(2)*diag(1)
3542  zln(5) = zln(5)*diag(3)
3543  zln(8) = zln(8)*diag(6)
3544  zln(5) = zln(5) - zln(8)*diag(5)
3545  zln(2) = zln(2) - zln(5)*diag(2) - zln(8)*diag(4)
3546 
3547  zln(6) = zz(6) - zz(3)*diag(2)
3548  zln(9) = zz(9) - zz(3)*diag(4) - zln(6)*diag(5)
3549  zln(3) = zz(3)*diag(1)
3550  zln(6) = zln(6)*diag(3)
3551  zln(9) = zln(9)*diag(6)
3552  zln(6) = zln(6) - zln(9)*diag(5)
3553  zln(3) = zln(3) - zln(6)*diag(2) - zln(9)*diag(4)
3554  end subroutine inv33
3555 
3556  !======================================================================!
3558  !======================================================================!
3559  subroutine inv6(Dsln,Ir)
3560  implicit none
3561  !------
3562  integer(kind=kint), intent(out):: Ir
3563  real(kind=kreal), intent(inout):: dsln(21)
3564  !------
3565  real(kind=kreal):: t(5)
3566 
3567  ir = 0
3568  dsln(1) = 1.0d0/dsln(1)
3569  t(1) = dsln(2)*dsln(1)
3570  dsln(3) = 1.0d0/(dsln(3)-t(1)*dsln(2))
3571  dsln(2) = t(1)
3572  dsln(5) = dsln(5) - dsln(4)*dsln(2)
3573  t(1) = dsln(4)*dsln(1)
3574  t(2) = dsln(5)*dsln(3)
3575  dsln(6) = 1.0d0/(dsln(6)-t(1)*dsln(4)-t(2)*dsln(5))
3576  dsln(4) = t(1)
3577  dsln(5) = t(2)
3578  dsln(8) = dsln(8) - dsln(7)*dsln(2)
3579  dsln(9) = dsln(9) - dsln(7)*dsln(4) - dsln(8)*dsln(5)
3580  t(1) = dsln(7)*dsln(1)
3581  t(2) = dsln(8)*dsln(3)
3582  t(3) = dsln(9)*dsln(6)
3583  dsln(10) = 1.0d0/(dsln(10)-t(1)*dsln(7)-t(2)*dsln(8)-t(3)*dsln(9))
3584  dsln(7) = t(1)
3585  dsln(8) = t(2)
3586  dsln(9) = t(3)
3587  dsln(12) = dsln(12) - dsln(11)*dsln(2)
3588  dsln(13) = dsln(13) - dsln(11)*dsln(4) - dsln(12)*dsln(5)
3589  dsln(14) = dsln(14) - dsln(11)*dsln(7) - dsln(12)*dsln(8) - dsln(13)*dsln(9)
3590  t(1) = dsln(11)*dsln(1)
3591  t(2) = dsln(12)*dsln(3)
3592  t(3) = dsln(13)*dsln(6)
3593  t(4) = dsln(14)*dsln(10)
3594  dsln(15) = 1.0d0/(dsln(15)-t(1)*dsln(11)-t(2)*dsln(12)-t(3)*dsln(13)-t(4)*dsln(14))
3595  dsln(11) = t(1)
3596  dsln(12) = t(2)
3597  dsln(13) = t(3)
3598  dsln(14) = t(4)
3599  dsln(17) = dsln(17) - dsln(16)*dsln(2)
3600  dsln(18) = dsln(18) - dsln(16)*dsln(4) - dsln(17)*dsln(5)
3601  dsln(19) = dsln(19) - dsln(16)*dsln(7) - dsln(17)*dsln(8) - dsln(18)*dsln(9)
3602  dsln(20) = dsln(20) - dsln(16)*dsln(11) - dsln(17)*dsln(12) - dsln(18)*dsln(13) - dsln(19)*dsln(14)
3603  t(1) = dsln(16)*dsln(1)
3604  t(2) = dsln(17)*dsln(3)
3605  t(3) = dsln(18)*dsln(6)
3606  t(4) = dsln(19)*dsln(10)
3607  t(5) = dsln(20)*dsln(15)
3608  dsln(21) = 1.0d0/(dsln(21)-t(1)*dsln(16)-t(2)*dsln(17)-t(3) *dsln(18)-t(4)*dsln(19)-t(5)*dsln(20))
3609  dsln(16) = t(1)
3610  dsln(17) = t(2)
3611  dsln(18) = t(3)
3612  dsln(19) = t(4)
3613  dsln(20) = t(5)
3614  end subroutine inv6
3615 
3616  !======================================================================!
3618  !======================================================================!
3619  subroutine inv66(Zln,Zz,Diag)
3620  implicit none
3621  !------
3622  real(kind=kreal), intent(in):: diag(21)
3623  real(kind=kreal), intent(in):: zz(36)
3624  real(kind=kreal), intent(out):: zln(36)
3625  !------
3626  integer(kind=kint):: i
3627 
3628  do i = 0, 5
3629  zln(i+7) = zz(i+7) - zz(i+1)*diag(2)
3630  zln(i+13) = zz(i+13) - zz(i+1)*diag(4) - zln(i+7)*diag(5)
3631  zln(i+19) = zz(i+19) - zz(i+1)*diag(7) - zln(i+7)*diag(8) - zln(i+13)*diag(9)
3632  zln(i+25) = zz(i+25) - zz(i+1)*diag(11) - zln(i+7)*diag(12) - zln(i+13)*diag(13)&
3633  - zln(i+19)*diag(14)
3634  zln(i+31) = zz(i+31) - zz(i+1)*diag(16) - zln(i+7)*diag(17) - zln(i+13)*diag(18)&
3635  - zln(i+19)*diag(19) - zln(i+25)*diag(20)
3636  zln(i+1) = zz(i+1)*diag(1)
3637  zln(i+7) = zln(i+7)*diag(3)
3638  zln(i+13) = zln(i+13)*diag(6)
3639  zln(i+19) = zln(i+19)*diag(10)
3640  zln(i+25) = zln(i+25)*diag(15)
3641  zln(i+31) = zln(i+31)*diag(21)
3642  zln(i+25) = zln(i+25) - zln(i+31)*diag(20)
3643  zln(i+19) = zln(i+19) - zln(i+31)*diag(19) - zln(i+25)*diag(14)
3644  zln(i+13) = zln(i+13) - zln(i+31)*diag(18) - zln(i+25)*diag(13)- zln(i+19)*diag(9)
3645  zln(i+7) = zln(i+7) - zln(i+31)*diag(17) - zln(i+25)*diag(12)- zln(i+19)*diag(8)&
3646  - zln(i+13)*diag(5)
3647  zln(i+1) = zln(i+1) - zln(i+31)*diag(16) - zln(i+25)*diag(11)- zln(i+19)*diag(7)&
3648  - zln(i+13)*diag(4) - zln(i+7)*diag(2)
3649  enddo
3650  end subroutine inv66
3651 
3652  !======================================================================!
3654  !======================================================================!
3655  subroutine invx(Dsln,Ndeg,Ir)
3656  implicit none
3657  !------
3658  integer(kind=kint), intent(in):: Ndeg
3659  integer(kind=kint), intent(out):: Ir
3660  real(kind=kreal), intent(inout):: dsln(*)
3661  !------
3662  integer(kind=kint):: i
3663  integer(kind=kint):: j
3664  integer(kind=kint):: k
3665  integer(kind=kint):: k0
3666  integer(kind=kint):: l
3667  integer(kind=kint):: l0
3668  integer(kind=kint):: ld
3669  integer(kind=kint):: ll
3670  real(kind=kreal):: t
3671  real(kind=kreal):: tem
3672 
3673  ir = 0
3674  l = 1
3675  dsln(1) = 1.0d0/dsln(1)
3676  do i = 2, ndeg
3677  ld = 0
3678  l0 = l
3679  do j = 1, i - 1
3680  l = l + 1
3681  do k = 1, j - 1
3682  ld = ld + 1
3683  dsln(l) = dsln(l) - dsln(l0+k)*dsln(ld)
3684  enddo
3685  ld = ld + 1
3686  enddo
3687  t = 0.0d0
3688  k0 = 0
3689  ll = 0
3690  do k = l - i + 2, l
3691  ll = ll + 1
3692  k0 = k0 + ll
3693  tem = dsln(k)*dsln(k0)
3694  t = t + tem*dsln(k)
3695  dsln(k) = tem
3696  enddo
3697  l = l + 1
3698  dsln(l) = dsln(l) - t
3699  dsln(l) = 1.0d0/dsln(l)
3700  enddo
3701  end subroutine invx
3702 
3703  !======================================================================!
3705  !======================================================================!
3706  subroutine invxx(Zln,Zz,Diag,Ndeg)
3707  implicit none
3708  !------
3709  integer(kind=kint), intent(in):: Ndeg
3710  real(kind=kreal), intent(in):: diag(*)
3711  real(kind=kreal), intent(in):: zz(ndeg,ndeg)
3712  real(kind=kreal), intent(out):: zln(ndeg,ndeg)
3713  !------
3714  integer(kind=kint):: joc
3715  integer(kind=kint):: l
3716  integer(kind=kint):: loc1
3717  integer(kind=kint):: m
3718  integer(kind=kint):: n
3719 
3720  zln = zz
3721  do l = 1, ndeg, 2
3722  joc = 0
3723  do m = 1, ndeg - 1
3724  joc = joc + m
3725  loc1 = joc + m
3726  do n = m + 1, ndeg
3727  zln(l,n) = zln(l,n) - zln(l,m)*diag(loc1)
3728  zln(l+1,n) = zln(l+1,n) - zln(l+1,m)*diag(loc1)
3729  loc1 = loc1 + n
3730  enddo
3731  enddo
3732  joc = 0
3733  do m = 1, ndeg
3734  joc = joc + m
3735  zln(l,m) = zln(l,m)*diag(joc)
3736  zln(l+1,m) = zln(l+1,m)*diag(joc)
3737  enddo
3738  do n = ndeg, 2, -1
3739  joc = joc - 1
3740  do m = n - 1, 1, -1
3741  zln(l,m) = zln(l,m) - zln(l,n)*diag(joc)
3742  zln(l+1,m) = zln(l+1,m) - zln(l+1,n)*diag(joc)
3743  joc = joc - 1
3744  enddo
3745  enddo
3746  enddo
3747  end subroutine invxx
3748 
3749  !======================================================================!
3751  !======================================================================!
3752  real(kind=kreal) function ddot(A,B,N)
3753  implicit none
3754  !------
3755  integer(kind=kint), intent(in):: N
3756  real(kind=kreal), intent(in):: a(n)
3757  real(kind=kreal), intent(in):: b(n)
3758  !------
3759  integer(kind=kint):: i
3760  real(kind=kreal):: s
3761 
3762  s = 0.0d0
3763  do i = 1, n
3764  s = s + a(i)*b(i)
3765  enddo
3766  ddot = s
3767  end function ddot
3768 
3769  !======================================================================!
3771  !======================================================================!
3772  subroutine d2dot(T,A,B,N)
3773  implicit none
3774  !------
3775  integer(kind=kint), intent(in):: N
3776  real(kind=kreal), intent(in):: a(4,*)
3777  real(kind=kreal), intent(in):: b(4,*)
3778  real(kind=kreal), intent(out):: t(4)
3779  !------
3780  integer(kind=kint):: jj
3781 
3782  t(1:4) = 0.0d0
3783 
3784  do jj = 1, n
3785  t(1) = t(1) + a(1,jj)*b(1,jj) + a(3,jj)*b(3,jj)
3786  t(2) = t(2) + a(2,jj)*b(1,jj) + a(4,jj)*b(3,jj)
3787  t(3) = t(3) + a(1,jj)*b(2,jj) + a(3,jj)*b(4,jj)
3788  t(4) = t(4) + a(2,jj)*b(2,jj) + a(4,jj)*b(4,jj)
3789  enddo
3790  end subroutine d2dot
3791 
3792  !======================================================================!
3794  !======================================================================!
3795  subroutine d3dot(T,A,B,N)
3796  implicit none
3797  !------
3798  integer(kind=kint), intent(in):: N
3799  real(kind=kreal), intent(in):: a(9,*)
3800  real(kind=kreal), intent(in):: b(9,*)
3801  real(kind=kreal), intent(out):: t(9)
3802  !------
3803  integer(kind=kint):: jj
3804 
3805  t(1:9) = 0.0d0
3806  do jj = 1, n
3807  t(1) = t(1) + a(1,jj)*b(1,jj) + a(4,jj)*b(4,jj) + a(7,jj)*b(7,jj)
3808  t(2) = t(2) + a(2,jj)*b(1,jj) + a(5,jj)*b(4,jj) + a(8,jj)*b(7,jj)
3809  t(3) = t(3) + a(3,jj)*b(1,jj) + a(6,jj)*b(4,jj) + a(9,jj)*b(7,jj)
3810 
3811  t(4) = t(4) + a(1,jj)*b(2,jj) + a(4,jj)*b(5,jj) + a(7,jj)*b(8,jj)
3812  t(5) = t(5) + a(2,jj)*b(2,jj) + a(5,jj)*b(5,jj) + a(8,jj)*b(8,jj)
3813  t(6) = t(6) + a(3,jj)*b(2,jj) + a(6,jj)*b(5,jj) + a(9,jj)*b(8,jj)
3814 
3815  t(7) = t(7) + a(1,jj)*b(3,jj) + a(4,jj)*b(6,jj) + a(7,jj)*b(9,jj)
3816  t(8) = t(8) + a(2,jj)*b(3,jj) + a(5,jj)*b(6,jj) + a(8,jj)*b(9,jj)
3817  t(9) = t(9) + a(3,jj)*b(3,jj) + a(6,jj)*b(6,jj) + a(9,jj)*b(9,jj)
3818  enddo
3819  end subroutine d3dot
3820 
3821  !======================================================================!
3823  !======================================================================!
3824  subroutine d3dotl(T,A,B,N)
3825  implicit none
3826  !------
3827  integer(kind=kint), intent(in):: N
3828  real(kind=kreal), intent(in):: a(9,*)
3829  real(kind=kreal), intent(in):: b(9,*)
3830  real(kind=kreal), intent(out):: t(6)
3831  !------
3832  integer(kind=kint):: jj
3833 
3834  t(1:6) = 0.0d0
3835  do jj = 1, n
3836  t(1) = t(1) + a(1,jj)*b(1,jj) + a(4,jj)*b(4,jj) + a(7,jj)*b(7,jj)
3837  t(2) = t(2) + a(2,jj)*b(1,jj) + a(5,jj)*b(4,jj) + a(8,jj)*b(7,jj)
3838 
3839  t(3) = t(3) + a(2,jj)*b(2,jj) + a(5,jj)*b(5,jj) + a(8,jj)*b(8,jj)
3840  t(4) = t(4) + a(3,jj)*b(1,jj) + a(6,jj)*b(4,jj) + a(9,jj)*b(7,jj)
3841 
3842  t(5) = t(5) + a(3,jj)*b(2,jj) + a(6,jj)*b(5,jj) + a(9,jj)*b(8,jj)
3843  t(6) = t(6) + a(3,jj)*b(3,jj) + a(6,jj)*b(6,jj) + a(9,jj)*b(9,jj)
3844  enddo
3845  end subroutine d3dotl
3846 
3847  !======================================================================!
3849  !======================================================================!
3850  subroutine dxdot(Ndeg,T,A,B,L)
3851  implicit none
3852  !------
3853  integer(kind=kint), intent(in):: L
3854  integer(kind=kint), intent(in):: Ndeg
3855  real(kind=kreal), intent(in):: a(ndeg,ndeg,*)
3856  real(kind=kreal), intent(in):: b(ndeg,ndeg,*)
3857  real(kind=kreal), intent(out):: t(ndeg,ndeg)
3858  !------
3859  integer(kind=kint):: jj
3860  integer(kind=kint):: k
3861  integer(kind=kint):: m
3862  integer(kind=kint):: n
3863 
3864  do n = 1, ndeg
3865  do m = 1, ndeg
3866  t(n,m) = 0.0d0
3867  do k = 1, ndeg
3868  do jj = 1, l
3869  t(n,m) = t(n,m) + a(n,k,jj)*b(m,k,jj)
3870  enddo
3871  enddo
3872  enddo
3873  enddo
3874  end subroutine dxdot
3875 
3876  !======================================================================!
3878  !======================================================================!
3879  subroutine dxdotl(Ndeg,T,A,B,L)
3880  implicit none
3881  !------
3882  integer(kind=kint), intent(in):: L
3883  integer(kind=kint), intent(in):: Ndeg
3884  real(kind=kreal), intent(in):: a(ndeg,ndeg,*)
3885  real(kind=kreal), intent(in):: b(ndeg,ndeg,*)
3886  real(kind=kreal), intent(out):: t(ndeg,ndeg)
3887  !------
3888  integer(kind=kint):: jj
3889  integer(kind=kint):: k
3890  integer(kind=kint):: m
3891  integer(kind=kint):: n
3892 
3893  do n = 1, ndeg
3894  do m = 1, n
3895  t(n,m) = 0.0d0
3896  do k = 1, ndeg
3897  do jj = 1, l
3898  t(n,m) = t(n,m) + a(n,k,jj)*b(m,k,jj)
3899  enddo
3900  enddo
3901  enddo
3902  enddo
3903  end subroutine dxdotl
3904 
3905  !======================================================================!
3907  !======================================================================!
3908  subroutine vprod(A,B,C,N)
3909  implicit none
3910  !------
3911  integer(kind=kint), intent(in):: N
3912  real(kind=kreal), intent(in):: a(n)
3913  real(kind=kreal), intent(in):: b(n)
3914  real(kind=kreal), intent(out):: c(n)
3915  !------
3916  c(1:n) = a(1:n)*b(1:n)
3917  end subroutine vprod
3918 
3919  !======================================================================!
3921  !======================================================================!
3922  subroutine v2prod(A,B,C,N)
3923  implicit none
3924  !------
3925  integer(kind=kint), intent(in):: N
3926  real(kind=kreal), intent(in):: a(4,n)
3927  real(kind=kreal), intent(in):: b(3,n)
3928  real(kind=kreal), intent(out):: c(4,n)
3929  !------
3930  integer(kind=kint):: i
3931 
3932  do i = 1, n
3933  c(3,i) = a(3,i) - a(1,i)*b(2,i)
3934  c(1,i) = a(1,i)*b(1,i)
3935  c(3,i) = c(3,i)*b(3,i)
3936  c(1,i) = c(1,i) - c(3,i)*b(2,i)
3937 
3938  c(4,i) = a(4,i) - a(2,i)*b(2,i)
3939  c(2,i) = a(2,i)*b(1,i)
3940  c(4,i) = c(4,i)*b(3,i)
3941  c(2,i) = c(2,i) - c(4,i)*b(2,i)
3942  enddo
3943  end subroutine v2prod
3944 
3945  !======================================================================!
3947  !======================================================================!
3948  subroutine v3prod(Zln,Diag,Zz,N)
3949  implicit none
3950  !------
3951  integer(kind=kint), intent(in):: N
3952  real(kind=kreal), intent(in):: diag(6,n)
3953  real(kind=kreal), intent(in):: zln(9,n)
3954  real(kind=kreal), intent(out):: zz(9,n)
3955  !------
3956  integer(kind=kint):: i
3957 
3958  do i = 1, n
3959  zz(4,i) = zln(4,i) - zln(1,i)*diag(2,i)
3960  zz(7,i) = zln(7,i) - zln(1,i)*diag(4,i) - zz(4,i)*diag(5,i)
3961  zz(1,i) = zln(1,i)*diag(1,i)
3962  zz(4,i) = zz(4,i)*diag(3,i)
3963  zz(7,i) = zz(7,i)*diag(6,i)
3964  zz(4,i) = zz(4,i) - zz(7,i)*diag(5,i)
3965  zz(1,i) = zz(1,i) - zz(4,i)*diag(2,i) - zz(7,i)*diag(4,i)
3966 
3967  zz(5,i) = zln(5,i) - zln(2,i)*diag(2,i)
3968  zz(8,i) = zln(8,i) - zln(2,i)*diag(4,i) - zz(5,i)*diag(5,i)
3969  zz(2,i) = zln(2,i)*diag(1,i)
3970  zz(5,i) = zz(5,i)*diag(3,i)
3971  zz(8,i) = zz(8,i)*diag(6,i)
3972  zz(5,i) = zz(5,i) - zz(8,i)*diag(5,i)
3973  zz(2,i) = zz(2,i) - zz(5,i)*diag(2,i) - zz(8,i)*diag(4,i)
3974 
3975  zz(6,i) = zln(6,i) - zln(3,i)*diag(2,i)
3976  zz(9,i) = zln(9,i) - zln(3,i)*diag(4,i) - zz(6,i)*diag(5,i)
3977  zz(3,i) = zln(3,i)*diag(1,i)
3978  zz(6,i) = zz(6,i)*diag(3,i)
3979  zz(9,i) = zz(9,i)*diag(6,i)
3980  zz(6,i) = zz(6,i) - zz(9,i)*diag(5,i)
3981  zz(3,i) = zz(3,i) - zz(6,i)*diag(2,i) - zz(9,i)*diag(4,i)
3982  enddo
3983  end subroutine v3prod
3984 
3985  !======================================================================!
3987  !======================================================================!
3988  subroutine vxprod(Ndeg,Ndegl,Zln,Diag,Zz,N)
3989  implicit none
3990  !------
3991  integer(kind=kint), intent(in):: Ndeg
3992  integer(kind=kint), intent(in):: Ndegl
3993  real(kind=kreal), intent(in):: diag(ndegl,n)
3994  real(kind=kreal), intent(in):: zln(ndeg*ndeg,n)
3995  real(kind=kreal), intent(out):: zz(ndeg*ndeg,n)
3996  !------
3997  integer(kind=kint):: i
3998  integer(kind=kint):: N
3999 
4000  do i = 1, n
4001  call invxx(zz(1,i),zln(1,i),diag(1,i),ndeg)
4002  enddo
4003  end subroutine vxprod
4004 
4005  !======================================================================!
4007  !======================================================================!
4008  subroutine vcopy(A,C,N)
4009  implicit none
4010  !------
4011  integer(kind=kint), intent(in):: N
4012  real(kind=kreal), intent(in):: a(n)
4013  real(kind=kreal), intent(out):: c(n)
4014  !------
4015  c = a
4016  end subroutine vcopy
4017 
4018  !======================================================================!
4020  ! (i/o)
4021  ! r_h_s on entry right hand side vector
4022  ! on exit solution vector
4023  ! iv communication array
4024  ! updates STAge
4025  !======================================================================!
4026  subroutine nusol0(R_h_s,FCT,Ir)
4027  implicit none
4028  !------
4029  real(kind=kreal), intent(inout):: r_h_s(:)
4030  type (cholesky_factor), intent(inout):: FCT
4031  integer(kind=kint), intent(out):: Ir
4032  !------
4033  integer(kind=kint):: ndegl
4034  integer(kind=kint):: ierror
4035  real(kind=kreal), pointer :: wk(:)
4036 
4037  if ( fct%STAge/=30 .and. fct%STAge/=40 ) then
4038  ir = 50
4039  return
4040  else
4041  ir = 0
4042  endif
4043 
4044  allocate (wk(fct%NDEg*fct%NEQns),stat=ierror)
4045  if ( ierror/=0 ) then
4046  write (*,*) "##Error: not enough memory"
4047  call hecmw_abort(hecmw_comm_get_comm())
4048  endif
4049  !rmiv
4050  ndegl = fct%NDEg*(fct%NDEg+1)
4051  ndegl = ndegl/2
4052 
4053  select case( fct%NDEg )
4054  case (1)
4055  call nusol1(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,fct%IPErm,r_h_s,wk,fct%NEQns,fct%NSTop)
4056  case (2)
4057  call nusol2(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,fct%IPErm,r_h_s,wk,fct%NEQns,fct%NSTop)
4058  case (3)
4059  call nusol3(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,fct%IPErm,r_h_s,wk,fct%NEQns,fct%NSTop)
4060  case (6)
4061  call nusolx(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,fct%IPErm,r_h_s,wk,fct%NEQns,fct%NSTop,fct%NDEg,ndegl)
4062  case default
4063  call nusolx(fct%XLNzr,fct%COLno,fct%DSLn,fct%ZLN,fct%DIAg,fct%IPErm,r_h_s,wk,fct%NEQns,fct%NSTop,fct%NDEg,ndegl)
4064  endselect
4065 
4066  fct%STAge = 40
4067  deallocate (wk)
4068  end subroutine nusol0
4069 
4070  !======================================================================!
4072  !======================================================================!
4073  subroutine nusol1(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
4074  implicit none
4075  !------
4076  integer(kind=kint), intent(in):: Neqns
4077  integer(kind=kint), intent(in):: Nstop
4078  integer(kind=kint), intent(in):: Xlnzr(:)
4079  integer(kind=kint), intent(in):: Colno(:)
4080  integer(kind=kint), intent(in):: Iperm(:)
4081  real(kind=kreal), intent(in):: zln(:)
4082  real(kind=kreal), intent(in):: diag(:)
4083  real(kind=kreal), intent(in):: dsln(:)
4084  real(kind=kreal), intent(inout):: b(:)
4085  real(kind=kreal), intent(out):: wk(:)
4086  !------
4087  integer(kind=kint):: i
4088  integer(kind=kint):: j
4089  integer(kind=kint):: joc
4090  integer(kind=kint):: k
4091  integer(kind=kint):: ke
4092  integer(kind=kint):: ks
4093 
4094  ! forward
4095  do i = 1, neqns
4096  wk(i) = b(iperm(i))
4097  enddo
4098  joc = 1
4099  do i = 1, neqns
4100  ks = xlnzr(i)
4101  ke = xlnzr(i+1) - 1
4102  if ( ke>=ks ) wk(i) = wk(i) - spdot2(wk,zln,colno,ks,ke)
4103  if ( i>nstop ) then
4104  wk(i) = wk(i) - ddot(wk(nstop:),dsln(joc:),i-nstop)
4105  joc = joc + i - nstop
4106  endif
4107  enddo
4108  do i = 1, neqns
4109  wk(i) = wk(i)*diag(i)
4110  enddo
4111  ! back ward
4112  do i = neqns, 1, -1
4113  if ( i>=nstop ) then
4114  do j = i - 1, nstop, -1
4115  joc = joc - 1
4116  wk(j) = wk(j) - wk(i)*dsln(joc)
4117  enddo
4118  endif
4119  ks = xlnzr(i)
4120  ke = xlnzr(i+1) - 1
4121  if ( ke>=ks ) then
4122  do k = ks, ke
4123  j = colno(k)
4124  wk(j) = wk(j) - wk(i)*zln(k)
4125  enddo
4126  endif
4127  enddo
4128  ! permutation
4129  do i = 1, neqns
4130  b(iperm(i)) = wk(i)
4131  enddo
4132  end subroutine nusol1
4133 
4134  !======================================================================!
4136  !======================================================================!
4137  subroutine nusol2(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
4138  implicit none
4139  !------
4140  integer(kind=kint), intent(in):: Neqns
4141  integer(kind=kint), intent(in):: Nstop
4142  integer(kind=kint), intent(in):: Xlnzr(:)
4143  integer(kind=kint), intent(in):: Colno(:)
4144  integer(kind=kint), intent(in):: Iperm(:)
4145  real(kind=kreal), intent(in):: zln(4,*)
4146  real(kind=kreal), intent(in):: diag(3,*)
4147  real(kind=kreal), intent(in):: dsln(4,*)
4148  real(kind=kreal), intent(inout):: b(2,*)
4149  real(kind=kreal), intent(out):: wk(2,*)
4150  !------
4151  integer(kind=kint):: i
4152  integer(kind=kint):: j
4153  integer(kind=kint):: joc
4154  integer(kind=kint):: k
4155  integer(kind=kint):: ke
4156  integer(kind=kint):: ks
4157 
4158  ! forward
4159  do i = 1, neqns
4160  wk(1,i) = b(1,iperm(i))
4161  wk(2,i) = b(2,iperm(i))
4162  enddo
4163  joc = 1
4164  do i = 1, neqns
4165  ks = xlnzr(i)
4166  ke = xlnzr(i+1) - 1
4167  if ( ke>=ks ) call s2pdot(wk(1,i),wk,zln,colno,ks,ke)
4168  if ( i>nstop ) then
4169  call d2sdot(wk(1,i),wk(1,nstop),dsln(1,joc),i-nstop)
4170  joc = joc + i - nstop
4171  endif
4172  enddo
4173  do i = 1, neqns
4174  wk(2,i) = wk(2,i) - wk(1,i)*diag(2,i)
4175  wk(1,i) = wk(1,i)*diag(1,i)
4176  wk(2,i) = wk(2,i)*diag(3,i)
4177  wk(1,i) = wk(1,i) - wk(2,i)*diag(2,i)
4178  enddo
4179  ! back ward
4180  do i = neqns, 1, -1
4181  if ( i>=nstop ) then
4182  do j = i - 1, nstop, -1
4183  joc = joc - 1
4184  wk(1,j) = wk(1,j) - wk(1,i)*dsln(1,joc) - wk(2,i)*dsln(2,joc)
4185  wk(2,j) = wk(2,j) - wk(1,i)*dsln(3,joc) - wk(2,i)*dsln(4,joc)
4186  enddo
4187  endif
4188  ks = xlnzr(i)
4189  ke = xlnzr(i+1) - 1
4190  if ( ke>=ks ) then
4191  do k = ks, ke
4192  j = colno(k)
4193  wk(1,j) = wk(1,j) - wk(1,i)*zln(1,k) - wk(2,i)*zln(2,k)
4194  wk(2,j) = wk(2,j) - wk(1,i)*zln(3,k) - wk(2,i)*zln(4,k)
4195  enddo
4196  endif
4197  enddo
4198  ! permutation
4199  do i = 1, neqns
4200  b(1,iperm(i)) = wk(1,i)
4201  b(2,iperm(i)) = wk(2,i)
4202  enddo
4203  end subroutine nusol2
4204 
4205  !======================================================================!
4207  !======================================================================!
4208  subroutine nusol3(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop)
4209  implicit none
4210  !------
4211  integer(kind=kint), intent(in):: Neqns
4212  integer(kind=kint), intent(in):: Nstop
4213  integer(kind=kint), intent(in):: Xlnzr(:)
4214  integer(kind=kint), intent(in):: Colno(:)
4215  integer(kind=kint), intent(in):: Iperm(:)
4216  real(kind=kreal), intent(in):: zln(9,*)
4217  real(kind=kreal), intent(in):: diag(6,*)
4218  real(kind=kreal), intent(in):: dsln(9,*)
4219  real(kind=kreal), intent(inout):: b(3,*)
4220  real(kind=kreal), intent(out):: wk(3,*)
4221  !------
4222  integer(kind=kint):: i
4223  integer(kind=kint):: j
4224  integer(kind=kint):: joc
4225  integer(kind=kint):: k
4226  integer(kind=kint):: ke
4227  integer(kind=kint):: ks
4228 
4229  ! forward
4230  do i = 1, neqns
4231  wk(1,i) = b(1,iperm(i))
4232  wk(2,i) = b(2,iperm(i))
4233  wk(3,i) = b(3,iperm(i))
4234  enddo
4235  joc = 1
4236  do i = 1, neqns
4237  ks = xlnzr(i)
4238  ke = xlnzr(i+1) - 1
4239  if ( ke>=ks ) call s3pdot(wk(1,i),wk,zln,colno,ks,ke)
4240  if ( i>nstop ) then
4241  call d3sdot(wk(1,i),wk(1,nstop),dsln(1,joc),i-nstop)
4242  joc = joc + i - nstop
4243  endif
4244  enddo
4245  do i = 1, neqns
4246  wk(2,i) = wk(2,i) - wk(1,i)*diag(2,i)
4247  wk(3,i) = wk(3,i) - wk(1,i)*diag(4,i) - wk(2,i)*diag(5,i)
4248  wk(1,i) = wk(1,i)*diag(1,i)
4249  wk(2,i) = wk(2,i)*diag(3,i)
4250  wk(3,i) = wk(3,i)*diag(6,i)
4251  wk(2,i) = wk(2,i) - wk(3,i)*diag(5,i)
4252  wk(1,i) = wk(1,i) - wk(2,i)*diag(2,i) - wk(3,i)*diag(4,i)
4253  enddo
4254  ! back ward
4255  do i = neqns, 1, -1
4256  if ( i>=nstop ) then
4257  do j = i - 1, nstop, -1
4258  joc = joc - 1
4259  wk(1,j) = wk(1,j) - wk(1,i)*dsln(1,joc) - wk(2,i)*dsln(2,joc) - wk(3,i)*dsln(3,joc)
4260  wk(2,j) = wk(2,j) - wk(1,i)*dsln(4,joc) - wk(2,i)*dsln(5,joc) - wk(3,i)*dsln(6,joc)
4261  wk(3,j) = wk(3,j) - wk(1,i)*dsln(7,joc) - wk(2,i)*dsln(8,joc) - wk(3,i)*dsln(9,joc)
4262  enddo
4263  endif
4264  ks = xlnzr(i)
4265  ke = xlnzr(i+1) - 1
4266  if ( ke>=ks ) then
4267  do k = ks, ke
4268  j = colno(k)
4269  wk(1,j) = wk(1,j) - wk(1,i)*zln(1,k) - wk(2,i)*zln(2,k) - wk(3,i)*zln(3,k)
4270  wk(2,j) = wk(2,j) - wk(1,i)*zln(4,k) - wk(2,i)*zln(5,k) - wk(3,i)*zln(6,k)
4271  wk(3,j) = wk(3,j) - wk(1,i)*zln(7,k) - wk(2,i)*zln(8,k) - wk(3,i)*zln(9,k)
4272  enddo
4273  endif
4274  enddo
4275  ! permutation
4276  do i = 1, neqns
4277  b(1,iperm(i)) = wk(1,i)
4278  b(2,iperm(i)) = wk(2,i)
4279  b(3,iperm(i)) = wk(3,i)
4280  enddo
4281  end subroutine nusol3
4282 
4283  !======================================================================!
4285  !======================================================================!
4286  subroutine nusolx(Xlnzr,Colno,Dsln,Zln,Diag,Iperm,B,Wk,Neqns,Nstop,Ndeg,Ndegl)
4287  implicit none
4288  !------
4289  integer(kind=kint), intent(in):: Ndeg
4290  integer(kind=kint), intent(in):: Ndegl
4291  integer(kind=kint), intent(in):: Neqns
4292  integer(kind=kint), intent(in):: Nstop
4293  integer(kind=kint), intent(in):: Xlnzr(:)
4294  integer(kind=kint), intent(in):: Colno(:)
4295  integer(kind=kint), intent(in):: Iperm(:)
4296  real(kind=kreal), intent(in):: zln(ndeg,ndeg,*)
4297  real(kind=kreal), intent(in):: diag(ndegl,*)
4298  real(kind=kreal), intent(in):: dsln(ndeg,ndeg,*)
4299  real(kind=kreal), intent(inout):: b(ndeg,*)
4300  real(kind=kreal), intent(out):: wk(ndeg,*)
4301  !------
4302  integer(kind=kint):: i
4303  integer(kind=kint):: j
4304  integer(kind=kint):: joc
4305  integer(kind=kint):: joc1
4306  integer(kind=kint):: k
4307  integer(kind=kint):: ke
4308  integer(kind=kint):: ks
4309  integer(kind=kint):: l
4310  integer(kind=kint):: loc1
4311  integer(kind=kint):: locd
4312  integer(kind=kint):: m
4313  integer(kind=kint):: n
4314 
4315  ! forward
4316  do l = 1, ndeg
4317  wk(l,1:neqns) = b(l,iperm(1:neqns))
4318  enddo
4319  joc = 1
4320  do i = 1, neqns
4321  ks = xlnzr(i)
4322  ke = xlnzr(i+1) - 1
4323  if ( ke>=ks ) call sxpdot(ndeg,wk(1,i),wk,zln,colno,ks,ke)
4324  if ( i>nstop ) then
4325  joc1 = i - nstop
4326  call dxsdot(ndeg,wk(1,i),wk(1,nstop),dsln(1,1,joc),joc1)
4327  joc = joc + joc1
4328  endif
4329  enddo
4330  do i = 1, neqns
4331  locd = 0
4332  do m = 1, ndeg - 1
4333  locd = locd + m
4334  loc1 = locd + m
4335  do n = m + 1, ndeg
4336  wk(n,i) = wk(n,i) - wk(m,i)*diag(loc1,i)
4337  loc1 = loc1 + n
4338  enddo
4339  enddo
4340  locd = 0
4341  do m = 1, ndeg
4342  locd = locd + m
4343  wk(m,i) = wk(m,i)*diag(locd,i)
4344  enddo
4345  do n = ndeg, 2, -1
4346  locd = locd - 1
4347  do m = n - 1, 1, -1
4348  wk(m,i) = wk(m,i) - wk(n,i)*diag(locd,i)
4349  locd = locd - 1
4350  enddo
4351  enddo
4352  enddo
4353  ! back ward
4354  do i = neqns, 1, -1
4355  if ( i>=nstop ) then
4356  do j = i - 1, nstop, -1
4357  joc = joc - 1
4358  do m = 1, ndeg
4359  do n = 1, ndeg
4360  wk(m,j) = wk(m,j) - wk(n,i)*dsln(n,m,joc)
4361  enddo
4362  enddo
4363  enddo
4364  endif
4365  ks = xlnzr(i)
4366  ke = xlnzr(i+1) - 1
4367  if ( ke>=ks ) then
4368  do k = ks, ke
4369  j = colno(k)
4370  do m = 1, ndeg
4371  do n = 1, ndeg
4372  wk(m,j) = wk(m,j) - wk(n,i)*zln(n,m,k)
4373  enddo
4374  enddo
4375  enddo
4376  endif
4377  enddo
4378  ! permutation
4379  do l = 1, ndeg
4380  b(l,iperm(1:neqns)) = wk(l,1:neqns)
4381  enddo
4382  end subroutine nusolx
4383 
4384  !======================================================================!
4386  !======================================================================!
4387  real(kind=kreal) function spdot2(B,Zln,Colno,Ks,Ke)
4388  implicit none
4389  !------
4390  integer(kind=kint), intent(in):: Ke
4391  integer(kind=kint), intent(in):: Ks
4392  integer(kind=kint), intent(in):: Colno(:)
4393  real(kind=kreal), intent(in):: zln(:)
4394  real(kind=kreal), intent(in):: b(:)
4395  !------
4396  integer(kind=kint):: j
4397  integer(kind=kint):: jj
4398  real(kind=kreal):: s
4399 
4400  s = 0.0d0
4401  do jj = ks, ke
4402  j = colno(jj)
4403  s = s + zln(jj)*b(j)
4404  enddo
4405  spdot2 = s
4406  end function spdot2
4407 
4408  !======================================================================!
4410  !======================================================================!
4411  subroutine s2pdot(Bi,B,Zln,Colno,Ks,Ke)
4412  implicit none
4413  !------
4414  integer(kind=kint), intent(in):: Ke
4415  integer(kind=kint), intent(in):: Ks
4416  integer(kind=kint), intent(in):: Colno(:)
4417  real(kind=kreal), intent(in):: zln(4,*)
4418  real(kind=kreal), intent(in):: b(2,*)
4419  real(kind=kreal), intent(inout):: bi(2)
4420  !------
4421  integer(kind=kint):: j
4422  integer(kind=kint):: jj
4423 
4424  do jj = ks, ke
4425  j = colno(jj)
4426  bi(1) = bi(1) - zln(1,jj)*b(1,j) - zln(3,jj)*b(2,j)
4427  bi(2) = bi(2) - zln(2,jj)*b(1,j) - zln(4,jj)*b(2,j)
4428  enddo
4429  end subroutine s2pdot
4430 
4431  !======================================================================!
4433  !======================================================================!
4434  subroutine s3pdot(Bi,B,Zln,Colno,Ks,Ke)
4435  implicit none
4436  !------
4437  integer(kind=kint), intent(in):: Ke
4438  integer(kind=kint), intent(in):: Ks
4439  integer(kind=kint), intent(in):: Colno(:)
4440  real(kind=kreal), intent(in):: zln(9,*)
4441  real(kind=kreal), intent(in):: b(3,*)
4442  real(kind=kreal), intent(inout):: bi(3)
4443  !------
4444  integer(kind=kint):: j
4445  integer(kind=kint):: jj
4446 
4447  do jj = ks, ke
4448  j = colno(jj)
4449  bi(1) = bi(1) - zln(1,jj)*b(1,j) - zln(4,jj)*b(2,j) - zln(7,jj)*b(3,j)
4450  bi(2) = bi(2) - zln(2,jj)*b(1,j) - zln(5,jj)*b(2,j) - zln(8,jj)*b(3,j)
4451  bi(3) = bi(3) - zln(3,jj)*b(1,j) - zln(6,jj)*b(2,j) - zln(9,jj)*b(3,j)
4452  enddo
4453  end subroutine s3pdot
4454 
4455  !======================================================================!
4457  !======================================================================!
4458  subroutine s6pdot(Bi,B,Zln,Colno,Ks,Ke)
4459  implicit none
4460  !------
4461  integer(kind=kint), intent(in):: Ke
4462  integer(kind=kint), intent(in):: Ks
4463  integer(kind=kint), intent(in):: Colno(:)
4464  real(kind=kreal), intent(in):: zln(36,*)
4465  real(kind=kreal), intent(in):: b(6,*)
4466  real(kind=kreal), intent(inout):: bi(6)
4467  !------
4468  integer(kind=kint):: j
4469  integer(kind=kint):: jj
4470 
4471  do jj = ks, ke
4472  j = colno(jj)
4473  bi(1) = bi(1) - zln(1,jj)*b(1,j) - zln(7,jj)*b(2,j) - zln(13,jj)*b(3,j)&
4474  - zln(19,jj)*b(4,j) - zln(25,jj)*b(5,j) - zln(31,jj)*b(6,j)
4475  bi(2) = bi(2) - zln(2,jj)*b(1,j) - zln(8,jj)*b(2,j) - zln(14,jj)*b(3,j)&
4476  - zln(20,jj)*b(4,j) - zln(26,jj)*b(5,j) - zln(32,jj)*b(6,j)
4477  bi(3) = bi(3) - zln(3,jj)*b(1,j) - zln(9,jj)*b(2,j) - zln(15,jj)*b(3,j)&
4478  - zln(21,jj)*b(4,j) - zln(27,jj)*b(5,j) - zln(33,jj)*b(6,j)
4479  bi(4) = bi(4) - zln(4,jj)*b(1,j) - zln(10,jj)*b(2,j) - zln(16,jj)*b(3,j)&
4480  - zln(22,jj)*b(4,j) - zln(28,jj)*b(5,j) - zln(34,jj)*b(6,j)
4481  bi(5) = bi(5) - zln(5,jj)*b(1,j) - zln(11,jj)*b(2,j) - zln(17,jj)*b(3,j)&
4482  - zln(23,jj)*b(4,j) - zln(29,jj)*b(5,j) - zln(35,jj)*b(6,j)
4483  bi(6) = bi(6) - zln(6,jj)*b(1,j) - zln(12,jj)*b(2,j) - zln(18,jj)*b(3,j)&
4484  - zln(25,jj)*b(4,j) - zln(30,jj)*b(5,j) - zln(36,jj)*b(6,j)
4485  enddo
4486  end subroutine s6pdot
4487 
4488  !======================================================================!
4490  !======================================================================!
4491  subroutine sxpdot(Ndeg,Bi,B,Zln,Colno,Ks,Ke)
4492  implicit none
4493  !------
4494  integer(kind=kint), intent(in):: Ke
4495  integer(kind=kint), intent(in):: Ks
4496  integer(kind=kint), intent(in):: Ndeg
4497  integer(kind=kint), intent(in):: Colno(:)
4498  real(kind=kreal), intent(in):: zln(ndeg,ndeg,*)
4499  real(kind=kreal), intent(in):: b(ndeg,*)
4500  real(kind=kreal), intent(inout):: bi(ndeg)
4501  !------
4502  integer(kind=kint):: j
4503  integer(kind=kint):: jj
4504  integer(kind=kint):: m
4505  integer(kind=kint):: n
4506 
4507  do jj = ks, ke
4508  j = colno(jj)
4509  do m = 1, ndeg
4510  do n = 1, ndeg
4511  bi(n) = bi(n) - zln(n,m,jj)*b(m,j)
4512  enddo
4513  enddo
4514  enddo
4515  end subroutine sxpdot
4516 
4517  !======================================================================!
4519  !======================================================================!
4520  subroutine d2sdot(Wi,A,B,N)
4521  implicit none
4522  !------
4523  integer(kind=kint), intent(in):: N
4524  real(kind=kreal), intent(in):: a(2,*)
4525  real(kind=kreal), intent(in):: b(4,*)
4526  real(kind=kreal), intent(inout):: wi(2)
4527  !------
4528  integer(kind=kint):: jj
4529 
4530  do jj = 1, n
4531  wi(1) = wi(1) - a(1,jj)*b(1,jj) - a(2,jj)*b(3,jj)
4532  wi(2) = wi(2) - a(1,jj)*b(2,jj) - a(2,jj)*b(4,jj)
4533  enddo
4534  end subroutine d2sdot
4535 
4536  !======================================================================!
4538  !======================================================================!
4539  subroutine d3sdot(Wi,A,B,N)
4540  implicit none
4541  !------
4542  integer(kind=kint), intent(in):: N
4543  real(kind=kreal), intent(in):: a(3,*)
4544  real(kind=kreal), intent(in):: b(9,*)
4545  real(kind=kreal), intent(inout):: wi(3)
4546  !------
4547  integer(kind=kint):: jj
4548 
4549  do jj = 1, n
4550  wi(1) = wi(1) - a(1,jj)*b(1,jj) - a(2,jj)*b(4,jj) - a(3,jj)*b(7,jj)
4551  wi(2) = wi(2) - a(1,jj)*b(2,jj) - a(2,jj)*b(5,jj) - a(3,jj)*b(8,jj)
4552  wi(3) = wi(3) - a(1,jj)*b(3,jj) - a(2,jj)*b(6,jj) - a(3,jj)*b(9,jj)
4553  enddo
4554  end subroutine d3sdot
4555 
4556  !======================================================================!
4558  !======================================================================!
4559  subroutine dxsdot(Ndeg,Wi,A,B,N)
4560  implicit none
4561  !------
4562  integer(kind=kint), intent(in):: Ndeg
4563  real(kind=kreal), intent(in):: a(ndeg,*)
4564  real(kind=kreal), intent(in):: b(ndeg,ndeg,*)
4565  real(kind=kreal), intent(inout):: wi(ndeg)
4566  integer(kind=kint), intent(inout):: N
4567  !------
4568  integer(kind=kint):: jj
4569  integer(kind=kint):: m
4570 
4571  do jj = 1, n
4572  do m = 1, ndeg
4573  do n = 1, ndeg
4574  wi(n) = wi(n) - b(n,m,jj)*a(m,jj)
4575  enddo
4576  enddo
4577  enddo
4578  end subroutine dxsdot
4579 
4580 end module hecmw_solver_direct
hecmw_matrix_dump::hecmw_mat_dump_solution
subroutine, public hecmw_mat_dump_solution(hecMAT)
Definition: hecmw_matrix_dump.f90:340
hecmw_ordering
HECMW_ORDERING is a program for fill-reducing ordering.
Definition: hecmw_ordering.F90:10
hecmw_util::hecmw_wtime
real(kind=kreal) function hecmw_wtime()
Definition: hecmw_util_f.F90:549
hecmw_util::hecmw_abort
subroutine hecmw_abort(comm)
Definition: hecmw_util_f.F90:534
hecmw_solver_direct::zpivot
subroutine zpivot(Neqns, Neqnsz, Nttbr, Jcol, Irow, Zpiv, Ir)
ZPIVOT.
Definition: hecmw_solver_direct.f90:432
hecmw_solver_direct::ptime
subroutine ptime(Cputim)
PTIME.
Definition: hecmw_solver_direct.f90:186
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::hecmwst_local_mesh
Definition: hecmw_util_f.F90:234
hecmw_solver_direct::hecmw_solve_direct
subroutine, public hecmw_solve_direct(hecMESH, hecMAT, Ifmsg)
HECMW_SOLVE_DIRECT is a program for the matrix solver.
Definition: hecmw_solver_direct.f90:61
hecmw_solver_direct
HECMW_SOLVE_DIRECT is a program for the matrix direct solver.
Definition: hecmw_solver_direct.f90:10
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
hecmw_matrix_ass
Definition: hecmw_mat_ass.f90:6
hecmw_matrix_dump
Definition: hecmw_matrix_dump.f90:6
hecmw_matrix_dump::hecmw_mat_dump
subroutine, public hecmw_mat_dump(hecMAT, hecMESH)
Definition: hecmw_matrix_dump.f90:32
hecmw_ordering::hecmw_ordering_gen
subroutine, public hecmw_ordering_gen(Neqns, Nttbr, Xadj, Adj0, Perm, Invp, opt, loglevel)
hecmw_ordering_gen
Definition: hecmw_ordering.F90:31
hecmw_util::hecmw_comm_get_comm
integer(kind=kint) function hecmw_comm_get_comm()
Definition: hecmw_util_f.F90:571
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444