FrontISTR  5.7.1
Large-scale structural analysis program with finit element method
hecmw_solver_direct_serial_lag.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 ! module for Parallel Direct Solver
6 module hecmw_solver_direct_serial_lag
7 
8  use m_child_matrix_lag
9  use m_irjc_matrix_lag
10  use my_hecmw_util_lag
11 
12  ! access control !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
13 
14  private ! default
15 
16  public hecmw_solve_direct_serial_lag ! only entry point of Parallel Direct Solver is public
17 
18  ! internal type definition !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
19 
20 
21  type dsinfo ! direct solver information
22  integer(kind=kint) :: ndeg ! dimension of small matrix
23  integer(kind=kint) :: neqns ! number of equations
24  integer(kind=kint) :: nstop ! beginning point of C
25  integer(kind=kint) :: stage ! calculation stage
26  integer(kind=kint) :: lncol ! length of col
27  integer(kind=kint) :: lndsln ! length of dsln
28 
29  integer(kind=kint), pointer :: zpiv(:) ! in zpivot()
30  integer(kind=kint), pointer :: iperm(:) ! permtation vector
31  integer(kind=kint), pointer :: invp(:) ! inverse permtation of iperm
32  integer(kind=kint), pointer :: parent(:) !
33  integer(kind=kint), pointer :: nch(:) !
34  integer(kind=kint), pointer :: xlnzr(:) ! ia index of whole sparse matrix. (neqns_t + 1)
35  integer(kind=kint), pointer :: colno(:) ! ja index of whole sparse matrix.
36 
37  real(kind=kreal), pointer :: diag(:,:) ! diagonal element
38  real(kind=kreal), pointer :: zln(:,:) ! non diagonal sparse
39  real(kind=kreal), pointer :: dsln(:,:) ! non diagonal dens
40  end type dsinfo
41 
42  ! internal global variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43 
44  real(kind=kreal), parameter :: rmin = 1.00d-200 ! for inv3() pivot
45 
46  integer :: imsg ! output file handler
47 
48  integer, parameter :: ilog = 16 ! according to FSTR
49  logical, parameter :: ldbg = .false.
50  !integer, parameter :: idbg = 52 ! according to FSTR
51  integer :: idbg = 10 ! debug output fort.*
52  logical :: lelap = .false. ! debug out
53 
54 contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
55 
56  subroutine hecmw_solve_direct_serial_lag(nrows, ilag_sta, nttbr, pointers, indices, values, b)
57  ! wrapper for parallel direct solver hecmw_solve_direct_parallel_internal()
58 
59  implicit none
60 
61 
62  ! arguments !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
63 
64  ! input
65  ! CRS style matrix
66  integer(kind=kint), intent(in) :: nrows ! number of rows of whole square matrix, including lagrange elements
67  integer(kind=kint), intent(in) :: ilag_sta ! row index of start point of lagrange elements
68  integer(kind=kint), intent(in) :: nttbr ! number of none zero elements of whole square matrix. include lagrange elements and both lower, upper matrix
69  integer(kind=kint), intent(in) :: pointers(:) ! whole matrix in CRS format. size is (0:nrows)
70  integer(kind=kint), intent(in) :: indices(:) ! whole matrix in CRS format. size is (nttbr)
71  real(kind=kreal), intent(in) :: values(:) ! whole matrix in CRS format. size is (nttbr)
72 
73  !input/output
74  real(kind=kreal), intent(inout) :: b(:) ! right hand side vector. result will return via this array.
75 
76 
77  ! internal valuables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
78 
79  ! stiffness matrix
80  type(irjc_square_matrix), target :: a0
81 
82  ! lagrange elements
83  type(irjc_mn_matrix), target :: lag
84 
85  ! degree of freedom is fixed as 1
86  integer(kind=kint), parameter :: ndeg_prm = 1
87 
88  ! for ASTOM matrix
89  type(irjc_square_matrix), target :: a0tmp
90  type(irjc_mn_matrix), target :: lagtmp
91 
92  ! misc
93  integer(kind=kint) :: i,j,k,l,ii,jj,kk,ll
94 
95  ! change CRS style matrix to irow jcol style !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
96 
97  ! set up a0
98  a0%neqns = ilag_sta - 1 ! non lagrange region only
99  a0%ndeg = ndeg_prm
100 
101  ! count ordinally A0 elements in upper triangle matrix from ASTOM
102  a0%nttbr = 0
103  do i= 1, nrows
104  do l= pointers(i), pointers(i+1)-1
105  j= indices(l)
106  if (j .le. a0%neqns) then
107  a0%nttbr = a0%nttbr+1
108  end if
109  enddo
110  enddo
111  !write (idbg,*) 'a0%nttbr', a0%nttbr
112  allocate( a0%irow(a0%nttbr), a0%jcol(a0%nttbr), a0%val(ndeg_prm,a0%nttbr))
113  allocate( a0tmp%irow(a0%nttbr), a0tmp%jcol(a0%nttbr), a0tmp%val(ndeg_prm,a0%nttbr))
114 
115  kk = 0
116  do i= 1, nrows
117  do l= pointers(i), pointers(i+1)-1
118  j= indices(l)
119  if (j .le. a0%neqns) then
120  !*Lower
121  kk = kk + 1
122  a0tmp%irow(kk) = i
123  a0tmp%jcol(kk) = j
124  a0tmp%val(1,kk)=values(l)
125  end if
126  enddo
127  enddo
128 
129  ! exchange irow and jcol, reordering
130  kk=0
131  do i=1,a0%neqns
132  do k=1,a0%nttbr
133  if (a0tmp%jcol(k) == i) then
134  kk=kk+1
135  a0%irow(kk)=a0tmp%jcol(k)
136  a0%jcol(kk)=a0tmp%irow(k)
137  a0%val(1,kk)=a0tmp%val(1,k)
138  end if
139  end do
140  end do
141 
142  ! set up lagrange elements !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143  lag%nrows = nrows - ilag_sta + 1 ! number of rows of lagrange region
144  lag%ncols = ilag_sta - 1 ! none lagrange region only
145  lag%ndeg = ndeg_prm
146 
147  ! count and allocate lagrange matrix
148  lag%nttbr = 0
149  do i= 1, nrows ! loop for whole matrix
150  do l= pointers(i), pointers(i+1)-1
151  j= indices(l)
152  if ((j .ge. ilag_sta) .and. (i .lt. ilag_sta) ) then
153  lag%nttbr = lag%nttbr+1
154  end if
155  enddo
156  enddo
157  !write (idbg,*) 'lag%nttbr', lag%nttbr
158  allocate( lag%irow(lag%nttbr), lag%jcol(lag%nttbr), lag%val(ndeg_prm,lag%nttbr))
159  allocate( lagtmp%irow(lag%nttbr), lagtmp%jcol(lag%nttbr), lagtmp%val(ndeg_prm,lag%nttbr))
160 
161  kk = 0
162  do i= 1, nrows
163  do l= pointers(i), pointers(i+1)-1
164  j= indices(l)
165  if ((j .ge. ilag_sta) .and. (i .lt. ilag_sta) ) then
166  kk = kk + 1
167  lagtmp%irow(kk) = i
168  lagtmp%jcol(kk) = j - ilag_sta +1
169  lagtmp%val(1,kk)=values(l)
170  end if
171  enddo
172  enddo
173 
174  ! exchange irow and jcol, reordering
175  kk=0
176  do i=1,lag%nrows
177  do k=1,lag%nttbr
178  if (lagtmp%jcol(k) == i) then
179  kk=kk+1
180  lag%irow(kk)=lagtmp%jcol(k)
181  lag%jcol(kk)=lagtmp%irow(k)
182  lag%val(1,kk)=lagtmp%val(1,k)
183  end if
184  end do
185  end do
186  call hecmw_solve_direct_serial_lag_in(a0,lag, b)
187 
188  return
189  end subroutine hecmw_solve_direct_serial_lag
190 
191  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
192 
193  subroutine hecmw_solve_direct_serial_lag_in(a0,lag, b)
194 
195  implicit none
196 
197 
198  type (irjc_square_matrix), intent(inout) :: a0 ! given left side matrix assembled from sp_matrix
199  type (irjc_mn_matrix), intent(inout) :: lag ! lagrange elements
200  real(kind=kreal), intent(inout) :: b(:) ! (a0%neqns) right hand side value vector include both original stifness matrix and followed by lagrange right hand side value.
201 
202  logical, save :: first_time = .true.
203 
204  integer(kind=kint) :: ierr
205 
206 
207  ! start !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
208 
209  imsg=99 ! set message file
210 
211  call sp_direct_parent(a0,lag,b)
212 
213  return
214  end subroutine hecmw_solve_direct_serial_lag_in
215 
216  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
217 
218  subroutine sp_direct_parent(a0, lag, b_in)
219 
220  implicit none
221 
222 
223  !I/O !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
224  ! given A0 x = b0
225  type (irjc_square_matrix), intent(inout) :: a0 ! given left side matrix assembled from sp_matrix
226  type (irjc_mn_matrix), intent(inout) :: lag ! lagrange elements
227  real(kind=kreal), intent(inout) :: b_in(:) ! (ndeg,neqns) right hand side value vector of equation.
228 
229  !internal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
230 
231 
232  ! for sp_direct_child
233  type (child_matrix) :: cm ! irow, jcol matrix
234  type (dsinfo) :: dsi ! direct solver info
235 
236 
237  integer(kind=kint) :: neqns_a ! number of eqns in A matrix
238  integer(kind=kint) :: neqns_lag ! number of eqns in D matrix
239  real(kind=kreal), pointer :: dsln(:,:) ! non-diagonal elements of dens D matrix
240  real(kind=kreal), pointer :: diag(:,:) ! diagonal elements of dens D matrix
241  type(child_matrix), pointer :: dm(:) !divided matrices
242 
243  real(kind=kreal), allocatable :: bd(:,:) ! for right hand side value
244 
245  ! internal use
246  real(kind=kreal), allocatable :: b(:,:)
247  real(kind=kreal), allocatable :: oldb(:,:)
248  real(kind=kreal), allocatable :: b_a(:)
249  real(kind=kreal), allocatable :: b_lag(:)
250 
251 
252 
253  logical, save :: nusol_ready = .false.
254  integer(kind=kint), save :: ndeg, nndeg, ndegt
255  integer(kind=kint), save :: neqns_c, iofst_a2, iofst_c, ndm
256 
257  ! misc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
258  integer(kind=kint) :: ierr
259  integer(kind=kint) :: i,j,k,l,m,n
260 
261 
262  real(kind=kreal), pointer :: spdslnval(:,:), bdbuf(:,:)
263  integer(kind=kint), pointer :: spdslnidx(:)
264  integer(kind=kint) :: nspdsln
265 
266 
267  !! temporary
268  integer(kind=kint), pointer :: iperm_all_inc_lag(:), part_all_inc_lag(:), iperm_rev_inc_lag(:)
269  integer(kind=kint) :: child_lag_nrows, child_lag_ncols, child_lag_nttbr, offset_irow
270  integer(kind=kint) :: ii, jj
271  real(kind=kreal), pointer :: dsln_lag(:,:)
272  real(kind=kreal), pointer :: diag_lag(:,:)
273  real(kind=kreal), allocatable :: wk(:), wk_d(:)
274  integer(kind=kint) :: ks, ke
275 
276 
277  ierr=0
278 
279 
280  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281  !
282  ! STEP01: get a0 from FEM data format hecMAT
283  !
284 
285  ndeg=a0%ndeg
286  nndeg=ndeg*ndeg
287  ndegt = (ndeg+1)*ndeg/2 !triangle element in diag
288 
289 
290  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291  !
292  ! STEP1X: LDU decompose of given A0
293  !
294  ! STEP11: build up matrix.
295  ! A is given a0
296  ! C is lagrange region
297  ! D is kept as dsln (off-diagonal), diag (diagonal).
298 
299 
300  neqns_a = a0%neqns
301  neqns_lag = lag%nrows
302  allocate(dsln_lag(1,neqns_lag*(neqns_lag - 1)/2)) ! size of lagrange
303  allocate(diag_lag(1,neqns_lag)) ! size of lagrange
304  dsln_lag=0.0d0 ! because of lagrange
305  diag_lag=0.0d0 ! because of lagrange
306 
307  ! set matrix for LDU decomposition
308 
309 
310  ! A matrix
311  cm%a%ndeg = a0%ndeg
312  cm%a%neqns = a0%neqns
313  cm%a%nttbr = a0%nttbr
314  allocate(cm%a%irow(cm%a%nttbr), cm%a%jcol(cm%a%nttbr), cm%a%val(1, cm%a%nttbr))
315  cm%a%irow(:) = a0%irow(:)
316  cm%a%jcol(:) = a0%jcol(:)
317  cm%a%val(1,:) = a0%val(1,:)
318 
319  ! C matrix (lagrange region)
320  cm%c%ndeg = lag%ndeg
321  cm%c%nttbr = lag%nttbr
322  cm%c%nrows = lag%nrows
323  cm%c%ncols = lag%ncols
324  allocate(cm%c%irow(cm%c%nttbr), cm%c%jcol(cm%c%nttbr), cm%c%val(1, cm%c%nttbr))
325  cm%c%irow(:) = lag%irow(:)
326  cm%c%jcol(:) = lag%jcol(:)
327  cm%c%val(1,:) = lag%val(1,:)
328 
329  cm%ndeg = cm%a%ndeg
330  cm%ista_c = cm%a%neqns+1
331  cm%neqns_t = cm%a%neqns + cm%c%nrows
332 
333 
334  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
335  !
336  ! STEP1x: LDU decompose for given matrix
337  !
338 
339  ! set up dsi for allocate matrix array for fill-in
340  call matini_para(cm, dsi, ierr)
341 
342  ! set real8 value
343  do i=1,cm%a%nttbr
344  call staij1(0, cm%a%irow(i), cm%a%jcol(i), cm%a%val(:,i), dsi, ierr)
345  end do
346  do i=1,cm%c%nttbr
347  ! call staij1(0, cm%c%irow(i)+cm%a%neqns, dsi%iperm(cm%c%jcol(i)), cm%c%val(:,i), dsi, ierr)
348  call staij1(0, cm%c%irow(i)+cm%a%neqns, cm%c%jcol(i), cm%c%val(:,i), dsi, ierr)
349  end do
350 
351  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
352  !
353  ! following STEP12-15 will be done in nufct0_child()
354  ! and return D region
355  !
356  ! STEP12: LDU decompose of A (1..nstop-1)
357  ! STEP13: LDU decompose of C (nstop..neqnsA+neqnsd)
358  ! STEP14: update D region.
359 
360  call nufct0_child(dsi, ierr, nspdsln, spdslnidx, spdslnval, diag_lag)
361 
362  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
363  !
364  ! STEP13: Receive D region and update D matrix as D' = D - D1' - D2' ...
365  !
366  ! D is receive as dens matrix, which format is given in s3um2() in serial solver.
367  ! to decompose this dens D matrix, use s3um3() on parent.
368 
369  ! off diagonal
370  do i=1,nspdsln
371  dsln_lag(:,spdslnidx(i)) = dsln_lag(:,spdslnidx(i)) + spdslnval(:,i) ! because of child process dsln is already substructed in s3um2()
372  end do
373 
374  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
375  !
376  ! STEP13: ! LDU decompose dens D
377  !
378  call nufct0_parent(dsln_lag, diag_lag, neqns_lag, ndeg)
379 
380  nusol_ready = .true.
381 
382  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
383  !
384  ! STEP2X: Solve Ax=b0
385  !
386 
387  ! forward substitution for A region
388 
389  ! set right hand side vector (b)
390  allocate(b(ndeg,a0%neqns+lag%nrows), stat=ierr)
391  if(ierr .ne. 0) then
392  call errtrp('stop due to allocation error.')
393  end if
394  do i=1,a0%neqns+lag%nrows
395  b(1,i)=b_in(i) !for ndeg=1
396  end do
397 
398  ! for verify
399  allocate(oldb(ndeg,a0%neqns+lag%nrows), stat=ierr)
400  if(ierr .ne. 0) then
401  call errtrp('stop due to allocation error.')
402  end if
403  do i=1, a0%neqns
404  oldb(ndeg,i)=b(ndeg,i)
405  end do
406  do i=a0%neqns+1, a0%neqns+lag%nrows
407  oldb(ndeg,i)=b(ndeg,i)
408  end do
409 
410  allocate(b_a(neqns_a))
411  do i=1,neqns_a
412  b_a(i)=b_in(i)
413  end do
414 
415  allocate(b_lag(neqns_lag))
416  do i=1, neqns_lag
417  b_lag(i)=b_in(neqns_a + i)
418  end do
419 
420  ! old nusol0_child !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
421  ! STEP22: forward substitution for A
422  allocate(wk(neqns_a+neqns_lag), stat=ierr)
423  if(ierr .ne. 0) then
424  call errtrp('stop due to allocation error.')
425  end if
426  wk = 0
427 
428  do i=1,neqns_a
429  wk(i)=b_a(dsi%iperm(i)) ! it sholud be permtated
430  end do
431 
432  ! STEP22: forward substitution for A
433  do i=1,neqns_a
434  ks=dsi%xlnzr(i)
435  ke=dsi%xlnzr(i+1)-1
436  if(ke.lt.ks) then ! logic inverted
437  cycle
438  end if
439  wk(i)=wk(i)-spdot2(wk,dsi%zln(1,:),dsi%colno,ks,ke)
440  end do
441 
442  ! STEP23: forward substitution for C and send it (yi) to parent
443  allocate(wk_d(dsi%nstop:dsi%neqns), stat=ierr)
444  if(ierr .ne. 0) then
445  call errtrp('stop due to allocation error.')
446  end if
447  wk_d=0
448 
449  do i=dsi%nstop,dsi%neqns
450  ks=dsi%xlnzr(i)
451  ke=dsi%xlnzr(i+1)-1
452  if(ke.lt.ks) then ! logic inverted
453  cycle
454  end if
455  wk_d(i)=wk_d(i)-spdot2(wk,dsi%zln(1,:),dsi%colno,ks,ke)
456  end do
457  ! Now wk_d is given. it should be used in parent.
458 
459 
460  ! end old nusol0_child !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
461 
462  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
463  !
464  ! STEP22 forward substitution for D region. and get Y
465  !
466  ! b1, b2... are sended to child processes.
467  ! bd is substituted to D in locally, and results from child processes
468  ! (C1-Y1 substitute, C2-Y2 substitute...) are receive from child processes.
469  ! these value also add for Yd
470 
471  ! update b_lag
472  do i=1, neqns_lag
473  b_lag(i)=b_lag(i) + wk_d(neqns_a + i)
474  end do
475 
476  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
477  !
478  ! STEP22 solve Ax=b for dens matrix, using updated bd
479  !
480  call nusol1_parent(dsln_lag(1,:), diag_lag(1,:), b_lag, neqns_lag)
481  !now b_lag is correct answer
482 
483 
484  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
485  !
486  ! STEP23 solve A region
487  !
488 
489  ! prepare Z
490  do i=1,neqns_a
491  wk(i)=wk(i)*dsi%diag(1,i)
492  end do
493 
494  ! D region is already solved
495  do i=1,neqns_lag
496  wk(neqns_a + i)=b_lag(i)
497  end do
498 
499  do i=dsi%neqns,1,-1
500  ks=dsi%xlnzr(i)
501  ke=dsi%xlnzr(i+1)-1
502  if(ke.lt.ks) then
503  cycle
504  end if
505  do k=ks,ke
506  j=dsi%colno(k)
507  wk(j)=wk(j)-wk(i)*dsi%zln(1,k)
508  end do
509  end do
510 
511  do i=1, neqns_a
512  b(1,dsi%iperm(i))=wk(i)
513  end do
514  do i=1, neqns_lag
515  b(1,neqns_a +i )=b_lag(i)
516  end do
517 
518 
519  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
520  !
521  ! STEP25 restore final result
522  !
523 
524  ! verify result
525  call verif0(ndeg, a0%neqns, a0%nttbr, a0%irow, a0%jcol, a0%val, lag%nrows, lag%nttbr, lag%irow, lag%jcol, lag%val, oldb, b) !verify result oldb will be broken.
526 
527  do i=1,a0%neqns+lag%nrows
528  b_in(i)=b(1,i)
529  end do
530 
531 
532  deallocate(spdslnidx, spdslnval)
533  deallocate(dsln_lag, diag_lag)
534 
535  deallocate(b, oldb, b_a, b_lag)
536 
537  deallocate(cm%a%irow, cm%a%jcol, cm%a%val)
538  deallocate(cm%c%irow, cm%c%jcol, cm%c%val)
539  deallocate(dsi%zpiv)
540  deallocate(dsi%iperm)
541  deallocate(dsi%invp)
542  deallocate(dsi%parent)
543  deallocate(dsi%nch)
544  deallocate(dsi%xlnzr)
545  deallocate(dsi%colno)
546  deallocate(dsi%diag)
547  deallocate(dsi%zln)
548  !deallocate(dsi%dsln) ! not allocated
549 
550  return
551  end subroutine sp_direct_parent
552 
553 
554  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
555  subroutine errtrp(mes)
556  character(*) mes
557  write(ilog,*) mes
558 
559  stop
560  end subroutine errtrp
561 
562  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
563 
564  subroutine matini_para(cm,dsi,ir)
565 
566  !----------------------------------------------------------------------
567  !
568  ! matini initializes storage for sparse matrix solver.
569  ! this routine is used for both symmetric matrices
570  ! and must be called once at the beginning
571  !
572  ! (i)
573  ! neqns number of unknowns
574  ! nttbr number of non0s, pattern of non-zero elements are
575  ! given like following.
576  ! nonz(A)={(i,j);i=irow(l),j=jcol(l); 1<= l <= nttbr}
577  ! irow
578  ! jcol to define non-zero pattern
579  ! lenv length of the array v (iv)
580  !
581  ! (o)
582  ! dsi matrix information
583  ! ir return code
584  ! =0 normal
585  ! =-1 non positive index
586  ! =1 too big index
587  ! =10 insufficient storage
588  !
589  !
590  ! stage 10 after initialization
591  ! 20 building up matrix
592  ! 30 after LU decomposition
593  ! 40 after solving
594  !
595  ! # coded by t.arakawa
596  ! # reviced by t.kitayama of Univ. Tokyo on 20071120
597  !
598  !----------------------------------------------------------------------
599 
600  implicit none
601 
602  type(child_matrix), intent(in) :: cm
603  type(dsinfo), intent(out) :: dsi
604  integer(kind=kint), intent(out) :: ir
605 
606  integer(kind=kint), pointer :: irow_a(:), jcol_a(:)
607  integer(kind=kint), pointer :: irow_c(:), jcol_c(:)
608 
609  integer(kind=kint), pointer :: ia(:) ! in stiaja() neqns+2
610  integer(kind=kint), pointer :: ja(:) ! in stiaja() 2*nttbr
611  integer(kind=kint), pointer :: jcpt(:) ! in stsmat() 2*nttbr
612  integer(kind=kint), pointer :: jcolno(:) ! in stsmat() 2*nttbr
613 
614  integer(kind=kint), pointer :: iperm_a(:)
615  integer(kind=kint), pointer :: invp_a(:)
616 
617  integer(kind=kint), pointer :: xlnzr_a(:)
618  integer(kind=kint), pointer :: colno_a(:)
619 
620  integer(kind=kint), pointer :: xlnzr_c(:)
621  integer(kind=kint), pointer :: colno_c(:)
622 
623 
624  integer(kind=kint), pointer :: adjncy(:) ! in genqmd() 2*nttbr
625  integer(kind=kint), pointer :: qlink(:) ! in genqmd() neqne+2
626  integer(kind=kint), pointer :: qsize(:) ! in genqmd() neqne+2
627  integer(kind=kint), pointer :: nbrhd(:) ! in genqmd() neqne+2
628  integer(kind=kint), pointer :: rchset(:) ! in genqmd() neqne+2
629 
630  integer(kind=kint), pointer :: cstr(:)
631 
632  integer(kind=kint), pointer :: adjt(:) ! in rotate() neqne+2
633  integer(kind=kint), pointer :: anc(:) ! in rotate() neqne+2
634 
635  integer(kind=kint), pointer :: lwk3arr(:)
636  integer(kind=kint), pointer :: lwk2arr(:)
637  integer(kind=kint), pointer :: lwk1arr(:)
638  integer(kind=kint), pointer :: lbtreearr(:,:) ! genbtq() (2,neqns+1)
639  integer(kind=kint), pointer :: lleafarr(:)
640  integer(kind=kint), pointer :: lxleafarr(:)
641  integer(kind=kint), pointer :: ladparr(:)
642  integer(kind=kint), pointer :: lpordrarr(:)
643 
644  integer(kind=kint) :: neqns_a, nttbr_a, neqns_a1, nstop, neqns_t, neqns_d, nttbr_c, ndeg
645  integer(kind=kint) :: lncol_a, lncol_c
646  integer(kind=kint) :: neqnsz, nofsub, izz, izz0, lnleaf ! dummy variables
647  integer(kind=kint) :: ir1
648  integer(kind=kint) :: i, j, k , ipass, ks, ke, ierr
649 
650  ndeg = cm%ndeg
651 
652  neqns_t = cm%neqns_t
653  neqns_d = cm%c%nrows
654 
655  neqns_a = cm%a%neqns
656  nttbr_a = cm%a%nttbr
657  irow_a => cm%a%irow
658  jcol_a => cm%a%jcol
659 
660  nttbr_c = cm%c%nttbr
661  irow_c => cm%c%irow
662  jcol_c => cm%c%jcol
663 
664  dsi%neqns=neqns_t ! because direct solver treat A + C as one matrix.
665  dsi%ndeg=ndeg
666 
667  neqns_a1=neqns_a+2
668  ir=0
669  ierr=0
670  izz0=0
671  !
672  ! set z pivot
673  !
674  allocate(dsi%zpiv(neqns_a), stat=ierr)
675  if(ierr .ne. 0) then
676  call errtrp('stop due to allocation error.')
677  end if
678  call zpivot(neqns_a,neqnsz,nttbr_a,jcol_a,irow_a,dsi%zpiv,ir1)
679  if(ir1.ne.0) then
680  ir=ir1
681  goto 1000
682  endif
683  !
684  ! build jcpt,jcolno
685  !
686  allocate(jcpt(2*nttbr_a), jcolno(2*nttbr_a), stat=ierr)
687  if(ierr .ne. 0) then
688  call errtrp('stop due to allocation error.')
689  end if
690  call stsmat(neqns_a,nttbr_a,irow_a,jcol_a,jcpt,jcolno)
691  !
692  ! build ia,ja
693  !
694  allocate(ia(neqns_a1), ja(2*nttbr_a), stat=ierr)
695  if(ierr .ne. 0) then
696  call errtrp('stop due to allocation error.')
697  end if
698  call stiaja(neqns_a, neqns_a,ia,ja,jcpt,jcolno)
699  !
700  ! get permutation vector iperm,invp
701  !
702 
703  ! setup identity permtation for C matrix
704  allocate(iperm_a(neqns_a), invp_a(neqns_a), stat=ierr)
705  if(ierr .ne. 0) then
706  call errtrp('stop due to allocation error.')
707  end if
708  call idntty(neqns_a,invp_a,iperm_a)
709 
710  ! reorder A matrix
711  allocate(adjncy(2*nttbr_a),qlink(neqns_a1),qsize(neqns_a1),nbrhd(neqns_a1),rchset(neqns_a1), stat=ierr)
712  if(ierr .ne. 0) then
713  call errtrp('stop due to allocation error.')
714  end if
715  allocate(lwk2arr(neqns_a1),lwk1arr(neqns_a1), stat=ierr)
716  if(ierr .ne. 0) then
717  call errtrp('stop due to allocation error.')
718  end if
719  call genqmd(neqns_a,ia,ja,iperm_a,invp_a,lwk1arr,lwk2arr,rchset,nbrhd,qsize,qlink,nofsub,adjncy)
720  deallocate(adjncy, qlink, qsize, nbrhd, rchset)
721 
722  ! set ia, ja
723  call stiaja(neqns_a, neqns_a, ia, ja, jcpt,jcolno)
724 
725  ! build up the parent vector parent vector will be saved in
726  ! work2 for a while
727  allocate(cstr(neqns_a1),adjt(neqns_a1), stat=ierr)
728  if(ierr .ne. 0) then
729  call errtrp('stop due to allocation error.')
730  end if
731 10 continue
732  call genpaq(ia,ja,invp_a,iperm_a,lwk2arr,neqns_a,cstr)
733 
734  ! build up the binary tree
735  allocate (lbtreearr(2,neqns_a1), stat=ierr)
736  if(ierr .ne. 0) then
737  call errtrp('stop due to allocation error.')
738  end if
739  call genbtq(ia, ja, invp_a, iperm_a,lwk2arr,lbtreearr,dsi%zpiv,izz,neqns_a)
740 
741  ! rotate the binary tree to avoid a zero pivot
742  if(izz.eq.0) goto 20
743  if(izz0.eq.0) izz0=izz
744  if(izz0.ne.izz) goto 30
745  call rotate(ia, ja, invp_a, iperm_a, lwk2arr,lbtreearr,izz,neqns_a,anc,adjt,ir1)
746  goto 10
747 30 continue
748  call bringu(dsi%zpiv,iperm_a, invp_a, lwk2arr,izz,neqns_a,ir1)
749  goto 10
750 
751  ! post ordering
752 20 continue
753  allocate(lwk3arr(0:neqns_a1),lpordrarr(neqns_a1),dsi%parent(neqns_a1), dsi%nch(neqns_a1), stat=ierr)
754  if(ierr .ne. 0) then
755  call errtrp('stop due to allocation error.')
756  end if
757  call posord(dsi%parent,lbtreearr,invp_a,iperm_a,lpordrarr,dsi%nch,neqns_a,lwk1arr,lwk2arr,lwk3arr)
758 
759  ! generate skeleton graph
760  allocate(lleafarr(nttbr_a),lxleafarr(neqns_a1),ladparr(neqns_a1), stat=ierr)
761  if(ierr .ne. 0) then
762  call errtrp('stop due to allocation error.')
763  end if
764  call gnleaf(ia, ja, invp_a, iperm_a, lpordrarr,dsi%nch,ladparr,lxleafarr,lleafarr,neqns_a,lnleaf)
765 
766  ! build up xlnzr,colno (this is the symbolic fct.)
767  nstop = cm%ista_c
768  call countclno(dsi%parent, lxleafarr, lleafarr, neqns_a, nstop, lncol_a, ir1) ! only for A
769  allocate(colno_a(lncol_a),xlnzr_a(neqns_a1), stat=ierr)
770  if(ierr .ne. 0) then
771  call errtrp('stop due to allocation error.')
772  end if
773  call gnclno(dsi%parent,lpordrarr,lxleafarr,lleafarr,xlnzr_a, colno_a, neqns_a, nstop,lncol_a,ir1) ! only for A
774  if(ir1.ne.0) then
775  ir=10
776  goto 1000
777  endif
778 
779  ! do symbolic LDU decomposition for C region.
780  call ldudecomposec(xlnzr_a,colno_a,invp_a,iperm_a, ndeg, nttbr_c, irow_c, &
781  jcol_c, cm%c%ncols, cm%c%nrows, xlnzr_c, colno_c, lncol_c)
782 
783  ! set calculated information to dsi.
784  allocate(dsi%xlnzr(neqns_t + 1), stat=ierr)
785  if(ierr .ne. 0) then
786  call errtrp('stop due to allocation error.')
787  end if
788  dsi%xlnzr(1:neqns_a)=xlnzr_a(:)
789  dsi%xlnzr(neqns_a+1:neqns_t+1)=xlnzr_c(:)+xlnzr_a(neqns_a+1)-1
790 
791  dsi%lncol=lncol_a + lncol_c
792  allocate(dsi%colno(lncol_a + lncol_c), stat=ierr)
793  if(ierr .ne. 0) then
794  call errtrp('stop due to allocation error.')
795  end if
796 
797  do i=1,lncol_a
798  dsi%colno(i)=colno_a(i)
799  end do
800  do i=1, lncol_c
801  dsi%colno(lncol_a + i)=colno_c(i)
802  end do
803 
804  allocate(dsi%invp(neqns_t), dsi%iperm(neqns_t), stat=ierr)
805  if(ierr .ne. 0) then
806  call errtrp('stop due to allocation error.')
807  end if
808  dsi%invp(1:neqns_a)=invp_a(1:neqns_a)
809  dsi%iperm(1:neqns_a)=iperm_a(1:neqns_a)
810  do i=neqns_a+1,neqns_t
811  dsi%invp(i)=i
812  dsi%iperm(i)=i
813  end do
814 
815  deallocate(xlnzr_a, colno_a, xlnzr_c, colno_c, invp_a, iperm_a)
816 
817  dsi%nstop=nstop
818  dsi%stage=10
819 1000 continue
820 
821  return
822  end subroutine matini_para
823 
824  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
825 
826  subroutine nufct0_child(dsi,ir, nspdsln, spdslnidx, spdslnval, diag_lag)
827 
828  implicit none
829  type(dsinfo), intent(inout) :: dsi
830  integer(kind=kint), intent(out) :: ir
831 
832  integer(kind=kint), intent(inout) :: nspdsln
833  real(kind=kreal), pointer :: spdslnval(:,:), bdbuf(:,:)
834  integer(kind=kint), pointer :: spdslnidx(:)
835 
836  real(kind=kreal), intent(inout) :: diag_lag(:,:)
837 
838  !
839  ! this performs Cholesky factorization
840  !
841 
842  if(dsi%stage.ne.20) then
843  ir=40
844  goto 1000
845  else
846  ir=0
847  endif
848  if(dsi%ndeg.eq.1) then
849  call nufct1_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%neqns,dsi%parent,dsi%nch,dsi%nstop,ir, &
850  nspdsln, spdslnidx, spdslnval, diag_lag)
851  else if(dsi%ndeg.eq.2) then
852  write(idbg,*) 'ndeg=1 only'
853  stop
854  else if(dsi%ndeg.eq.3) then
855  write(idbg,*) 'ndeg=1 only'
856  stop
857  else if(dsi%ndeg.eq.6) then
858  write(idbg,*) 'ndeg=1 only'
859  stop
860  else
861  write(idbg,*) 'ndeg=1 only'
862  stop
863  end if
864 
865  dsi%stage=30
866 1000 continue
867  return
868  end subroutine nufct0_child
869 
870  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
871 
872  subroutine nufct1_child(xlnzr,colno,zln,diag,neqns,parent,nch,nstop,ir, nspdsln, spdslnidx, spdslnval, diag_lag)
873 
874  implicit none
875  integer(kind=kint), intent(in) :: xlnzr(:),colno(:),parent(:)
876  integer(kind=kint), intent(in) :: neqns, nstop, ir
877  integer(kind=kint), intent(out) :: nch(:)
878  real(kind=kreal), intent(out) :: zln(:,:),diag(:,:) !zln(1,:), diag(1,:)
879 
880  integer(kind=kint) :: neqns_c
881  integer(kind=kint) :: i,j,k,l, ic,ierr,imp
882  integer(kind=kint) :: nspdsln
883  integer(kind=kint), pointer :: spdslnidx(:)
884  real(kind=kreal), pointer :: spdslnval(:,:)
885  real(kind=kreal), intent(out) :: diag_lag(:,:) ! diag(1,:)
886 
887  !----------------------------------------------------------------------
888  !
889  ! nufct1 performs cholesky factorization in row order for ndeg=1
890  !
891  ! (i) xlnzr,colno,zln,diag
892  ! symbolicaly factorized
893  !
894  ! (o) zln,diag,dsln
895  !
896  ! #coded by t.arakawa
897  !
898  !----------------------------------------------------------------------
899 
900  !
901  ! phase I
902  ! LDU decompose of A (1..nstop-1)
903  !
904  diag(1,1)=1.0d0/diag(1,1)
905  l=parent(1)
906  nch(l)=nch(l)-1
907  nch(1)=-1
908  do ic=2,nstop-1
909  call sum(ic,xlnzr,colno,zln(1,:),diag(1,:),nch,parent,neqns)
910  enddo
911 
912  !
913  ! phase II
914  ! LDU decompose of C (nstop..neqnsA+neqnsd)
915  !
916  do ic=nstop,neqns
917  call sum1(ic,xlnzr,colno,zln(1,:),diag(1,:),parent,neqns)
918  enddo
919 
920  !
921  ! phase III
922  ! Update D region.
923  !
924 
925  ! clear dummy diagonal value for D region
926  do i=nstop,neqns
927  diag(:,i)=0.0
928  end do
929 
930  neqns_c = neqns - nstop + 1
931  call sum2_child(neqns,nstop,xlnzr,colno,zln(1,:),diag(1,:),spdslnidx,spdslnval,nspdsln)
932  ! send D region to parent
933  ! imp = m_pds_procinfo%imp
934  ! call MPI_SEND(nspdsln, 1,MPI_INTEGER,IMP,1,MPI_COMM_WORLD,ierr)
935  ! call MPI_SEND(spdslnidx, nspdsln,MPI_INTEGER,IMP,1,MPI_COMM_WORLD,ierr)
936  ! call MPI_SEND(spdslnval, nspdsln,MPI_REAL8,IMP,1,MPI_COMM_WORLD,ierr)
937  ! call MPI_SEND(diag(1,nstop), neqns_c,MPI_REAL8,IMP,1,MPI_COMM_WORLD,ierr)
938 
939  do i=1, neqns_c
940  diag_lag(1,i) = diag(1,nstop+i-1) ! lagrange region
941  end do
942 
943  return
944  end subroutine nufct1_child
945 
946 
947 
948  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
949 
950  subroutine nufct0_parent(dsln, diag, neqns, ndeg)
951  ! select LDU decomposer for dens matrix according to ndeg
952 
953  implicit none
954 
955  real(kind=kreal), intent(inout) :: dsln(:,:)
956  real(kind=kreal), intent(inout) :: diag(:,:)
957  integer(kind=kint), intent(in) :: neqns, ndeg
958 
959  integer(kind=kint) :: ndegl
960 
961  if (ndeg .eq. 1) then
962  call sum3(neqns, dsln(1,:), diag(1,:))
963  else if (ndeg .eq. 3) then
964  write(idbg,*) 'ndeg=1 only'
965  stop
966  else
967  write(idbg,*) 'ndeg=1 only'
968  stop
969  end if
970 
971  return
972  end subroutine nufct0_parent
973 
974  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
975 
976  subroutine nusol0_parent(dsln, diag, b, neqns, ndeg)
977  ! select solvers according to ndeg
978 
979  implicit none
980 
981  real(kind=kreal), intent(in) :: dsln(:,:)
982  real(kind=kreal), intent(in) :: diag(:,:)
983  real(kind=kreal), intent(inout) :: b(:,:)
984 
985  integer(kind=kint), intent(in) :: neqns, ndeg
986 
987  if (ndeg .eq. 1) then
988  call nusol1_parent(dsln(1,:), diag(1,:), b(1,:), neqns)
989  else if (ndeg .eq. 3) then
990  write(idbg,*) 'ndeg=1 only'
991  stop
992  else
993  write(idbg,*) 'ndeg=1 only'
994  stop
995  end if
996 
997  return
998  end subroutine nusol0_parent
999 
1000  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1001 
1002  subroutine nusol1_parent(dsln, diag, b, neqns)
1003  ! solve Ax=b for dens matrix with ndeg=1
1004  ! require dsln, diag is already LDU decomposed.
1005  ! currently not tested 20071121
1006 
1007  implicit none
1008 
1009  real(kind=kreal), intent(in) :: dsln(:) !((neqns+1)*neqns/2)
1010  real(kind=kreal), intent(in) :: diag(:) !(neqns)
1011  real(kind=kreal), intent(inout) :: b(:) !(3,neqns)
1012  integer(kind=kint), intent(in) :: neqns
1013 
1014  integer(kind=kint) :: i,j,k,l,loc
1015 
1016  ! forward substitution
1017  do i=2,neqns
1018  k=(i-1)*(i-2)/2 + 1 ! first element of i'th row.
1019  b(i)=b(i)-dot_product(b(1:i-1),dsln(k:k+i-2))
1020  end do
1021 
1022  ! divide by D (because of diag is already inverted (1/Dii))
1023  b(:)=b(:)*diag(:)
1024 
1025  ! Backward substitution.
1026  ! Substitute Zi into D and get Xd results.
1027  loc=(neqns-1)*neqns/2
1028  do i=neqns,1,-1
1029  do j=i-1,1,-1
1030  b(j)=b(j)-b(i)*dsln(loc)
1031  loc=loc-1
1032  end do
1033  end do
1034 
1035  return
1036  end subroutine nusol1_parent
1037 
1038  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1039 
1040 
1041  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1042 
1043  subroutine zpivot(neqns,neqnsz,nttbr,jcol,irow,zpiv,ir)
1044 
1045  implicit none
1046 
1047  integer(kind=kint), intent(in) :: jcol(:),irow(:)
1048  integer(kind=kint), intent(out) :: zpiv(:)
1049  integer(kind=kint), intent(in) :: neqns,nttbr
1050  integer(kind=kint), intent(out) :: neqnsz,ir
1051 
1052  integer(kind=kint) :: i,j,k,l
1053 
1054  ir=0
1055  do l=1,neqns
1056  zpiv(l)=1
1057  enddo
1058 
1059  do l=1,nttbr
1060  i=irow(l)
1061  j=jcol(l)
1062  if(i.le.0.or.j.le.0) then
1063  ir=-1
1064  goto 1000
1065  elseif(i.gt.neqns.or.j.gt.neqns) then
1066  ir=1
1067  goto 1000
1068  endif
1069  if(i.eq.j) zpiv(i)=0
1070  enddo
1071 
1072  do i=neqns,1,-1
1073  if(zpiv(i).eq.0) then
1074  neqnsz=i
1075  goto 320
1076  endif
1077  enddo
1078 320 continue
1079 1000 continue
1080  if(ldbg) write(idbg,*) '# zpivot ########################'
1081  if(ldbg) write(idbg,60) (zpiv(i),i=1,neqns)
1082 60 format(20i3)
1083  return
1084  end subroutine zpivot
1085 
1086  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1087 
1088  subroutine stsmat(neqns,nttbr,irow,jcol,jcpt,jcolno)
1089 
1090  implicit none
1091 
1092  integer(kind=kint), intent(in) :: irow(:), jcol(:)
1093  integer(kind=kint), intent(out) :: jcpt(:), jcolno(:)
1094  integer(kind=kint), intent(in) :: neqns, nttbr
1095 
1096  integer(kind=kint) :: i,j,k,l,loc,locr
1097 
1098  do i=1,2*nttbr
1099  jcpt(i)=0
1100  jcolno(i)=0
1101  enddo
1102  do i=1,neqns
1103  jcpt(i)=i+neqns
1104  jcolno(i+neqns)=i
1105  enddo
1106 
1107  k=2*neqns
1108  do 100 l=1,nttbr
1109  i=irow(l)
1110  j=jcol(l)
1111  if(i.eq.j) goto 100
1112  loc=jcpt(i)
1113  locr=i
1114 110 continue
1115  if(loc.eq.0) goto 120
1116  if(jcolno(loc).eq.j) then
1117  goto 100
1118  elseif(jcolno(loc).gt.j) then
1119  goto 130
1120  endif
1121  locr=loc
1122  loc=jcpt(loc)
1123  goto 110
1124 120 continue
1125  k=k+1
1126  jcpt(locr)=k
1127  jcolno(k)=j
1128  goto 150
1129 130 continue
1130  k=k+1
1131  jcpt(locr)=k
1132  jcpt(k)=loc
1133  jcolno(k)=j
1134 150 continue
1135  loc=jcpt(j)
1136  locr=j
1137 160 continue
1138  if(loc.eq.0) goto 170
1139  if(jcolno(loc).eq.i) then
1140  goto 100
1141  elseif(jcolno(loc).gt.i) then
1142  goto 180
1143  endif
1144  locr=loc
1145  loc=jcpt(loc)
1146  goto 160
1147 170 continue
1148  k=k+1
1149  jcpt(locr)=k
1150  jcolno(k)=i
1151  goto 100
1152 180 continue
1153  k=k+1
1154  jcpt(locr)=k
1155  jcpt(k)=loc
1156  jcolno(k)=i
1157 100 continue
1158  if(ldbg) then
1159  write(idbg,*) 'jcolno'
1160  write(idbg,60) (jcolno(i),i=1,k)
1161  write(idbg,*) 'jcpt'
1162  write(idbg,60) (jcpt(i),i=1,k)
1163 60 format(10i7)
1164  endif
1165  return
1166  end subroutine stsmat
1167 
1168  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1169  subroutine stiaja(neqns,neqnsz,ia,ja,jcpt,jcolno)
1170 
1171  implicit none
1172  !
1173  ! coded by t.arakawa
1174  !
1175  integer(kind=kint), intent(in) :: jcpt(:),jcolno(:)
1176  integer(kind=kint), intent(out) :: ia(:),ja(:)
1177  integer(kind=kint), intent(in) :: neqns, neqnsz
1178 
1179  integer(kind=kint) :: i,j,k,l,ii,loc
1180  !
1181 
1182  ia(1)=1
1183  l=0
1184  do k=1,neqns
1185  loc=jcpt(k)
1186 110 continue
1187  if(loc.eq.0) goto 120
1188  ii=jcolno(loc)
1189  if(ii.eq.k.or.ii.gt.neqnsz) goto 130
1190  l=l+1
1191  ja(l)=ii
1192 130 continue
1193  loc=jcpt(loc)
1194  goto 110
1195 120 ia(k+1)=l+1
1196  enddo
1197  if(ldbg) then
1198  write(idbg,*) 'stiaja(): ia '
1199  write(idbg,60) (ia(i),i=1,neqns+1)
1200  write(idbg,*) 'stiaja(): ja '
1201  write(idbg,60) (ja(i),i=1,ia(neqns+1))
1202  endif
1203 60 format(10i7)
1204  return
1205  end subroutine stiaja
1206 
1207  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1208 
1209  subroutine idntty(neqns,invp,iperm)
1210 
1211  implicit none
1212 
1213  integer(kind=kint), intent(out) :: invp(:),iperm(:)
1214  integer(kind=kint), intent(in) :: neqns
1215 
1216  integer(kind=kint) :: i
1217 
1218  do i=1,neqns
1219  invp(i)=i
1220  iperm(i)=i
1221  enddo
1222  return
1223  end subroutine idntty
1224 
1225  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1226 
1227  subroutine genqmd(neqns,xadj,adj0,perm,invp,deg,marker,rchset,nbrhd,qsize,qlink,nofsub,adjncy)
1228 
1229  implicit none
1230 
1231  integer(kind=kint), intent(in) :: adj0(:),xadj(:)
1232  integer(kind=kint), intent(out) :: rchset(:),nbrhd(:),adjncy(:),perm(:),invp(:),deg(:),marker(:),qsize(:),qlink(:)
1233  integer(kind=kint), intent(in) :: neqns
1234  integer(kind=kint), intent(out) :: nofsub
1235 
1236  integer(kind=kint) :: inode,ip,irch,mindeg,nhdsze,node,np,num,nump1,nxnode,rchsze,search,thresh,ndeg
1237  integer(kind=kint) :: i,j,k,l
1238 
1239  mindeg=neqns
1240  nofsub=0
1241  do i=1,xadj(neqns+1)-1
1242  adjncy(i)=adj0(i)
1243  enddo
1244  do node=1,neqns
1245  perm(node)=node
1246  invp(node)=node
1247  marker(node)=0
1248  qsize(node)=1
1249  qlink(node)=0
1250  ndeg=xadj(node+1)-xadj(node)
1251  deg(node)=ndeg
1252  if(ndeg.lt.mindeg) mindeg=ndeg
1253  enddo
1254 
1255  num=0
1256 200 search=1
1257  thresh=mindeg
1258  mindeg=neqns
1259 300 nump1=num+1
1260  if(nump1.gt.search) search=nump1
1261  do 400 j=search,neqns
1262  node=perm(j)
1263  if(marker(node).lt.0) goto 400
1264  ndeg=deg(node)
1265  if(ndeg.le.thresh) goto 500
1266  if(ndeg.lt.mindeg) mindeg=ndeg
1267 400 continue
1268  goto 200
1269 
1270 500 search=j
1271  nofsub=nofsub+deg(node)
1272  marker(node)=1
1273  call qmdrch(node,xadj,adjncy,deg,marker,rchsze,rchset,nhdsze,nbrhd)
1274  nxnode=node
1275 600 num=num+1
1276  np=invp(nxnode)
1277  ip=perm(num)
1278  perm(np)=ip
1279  invp(ip)=np
1280  perm(num)=nxnode
1281  invp(nxnode)=num
1282  deg(nxnode)=-1
1283  nxnode=qlink(nxnode)
1284  if(nxnode.gt.0) goto 600
1285  if(rchsze.le.0) goto 800
1286  !
1287  call qmdupd(xadj,adjncy,rchsze,rchset,deg,qsize,qlink,marker,rchset(rchsze+1:),nbrhd(nhdsze+1:))
1288  marker(node)=0
1289  do 700 irch=1,rchsze
1290  inode=rchset(irch)
1291  if(marker(inode).lt.0) goto 700
1292  marker(inode)=0
1293  ndeg=deg(inode)
1294  if(ndeg.lt.mindeg) mindeg=ndeg
1295  if(ndeg.gt.thresh) goto 700
1296  mindeg=thresh
1297  thresh=ndeg
1298  search=invp(inode)
1299 700 continue
1300  if(nhdsze.gt.0) call qmdot(node,xadj,adjncy,marker,rchsze,rchset,nbrhd)
1301 800 if(num.lt.neqns) goto 300
1302  return
1303  end subroutine genqmd
1304 
1305  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1306 
1307  subroutine genpaq(xadj,adjncy,invp,iperm,parent,neqns,ancstr)
1308 
1309  implicit none
1310 
1311  integer(kind=kint), intent(in) :: xadj(:),adjncy(:),invp(:),iperm(:)
1312  integer(kind=kint), intent(out) :: parent(:),ancstr(:)
1313  integer(kind=kint), intent(in) :: neqns
1314 
1315  integer(kind=kint) :: i,j,k,l,ip,it
1316 
1317  do i=1,neqns
1318  parent(i)=0
1319  ancstr(i)=0
1320  ip=iperm(i)
1321  do 110 k=xadj(ip),xadj(ip+1)-1
1322  l=invp(adjncy(k))
1323  if(l.ge.i) goto 110
1324 112 continue
1325  if(ancstr(l).eq.0) goto 111
1326  if(ancstr(l).eq.i) goto 110
1327  it=ancstr(l)
1328  ancstr(l)=i
1329  l=it
1330  goto 112
1331 111 continue
1332  ancstr(l)=i
1333  parent(l)=i
1334 110 continue
1335  enddo
1336  do i=1,neqns
1337  if(parent(i).eq.0) parent(i)=neqns+1
1338  enddo
1339  parent(neqns+1)=0
1340  if(ldbg) write(idbg,6010)
1341  if(ldbg) write(idbg,6000) (i,parent(i),i=1,neqns)
1342 6000 format(2i6)
1343 6010 format(' parent')
1344  return
1345  end subroutine genpaq
1346 
1347  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1348 
1349  subroutine genbtq(xadj,adjncy,invp,iperm,parent,btree,zpiv,izz,neqns)
1350 
1351  implicit none
1352 
1353  integer(kind=kint), intent(in) :: xadj(:),adjncy(:),parent(:),invp(:),iperm(:),zpiv(:)
1354  integer(kind=kint), intent(out) :: btree(:,:) ! btree is (2,:)
1355  integer(kind=kint), intent(in) :: neqns
1356  integer(kind=kint), intent(out) :: izz
1357 
1358  integer(kind=kint) :: i,j,k,l,ip,ib,inext
1359 
1360  do i=1,neqns+1
1361  btree(1,i)=0
1362  btree(2,i)=0
1363  enddo
1364  do 100 i=1,neqns+1
1365  ip=parent(i)
1366  if(ip.le.0) goto 100
1367  ib=btree(1,ip)
1368  if(ib.eq.0) then
1369  btree(1,ip)=i
1370  else
1371 101 continue
1372  inext=btree(2,ib)
1373  if(inext.eq.0) then
1374  btree(2,ib)=i
1375  else
1376  ib=inext
1377  goto 101
1378  endif
1379  endif
1380 100 continue
1381  !
1382  ! find zeropivot
1383  !
1384  do i=1,neqns
1385  if(zpiv(i).ne.0) then
1386  if(btree(1,invp(i)).eq.0) then
1387  izz=i
1388  goto 210
1389  endif
1390  endif
1391  enddo
1392  izz=0
1393 210 continue
1394  if(ldbg) write(idbg,6010)
1395  if(ldbg) write(idbg,6000) (i,btree(1,i),btree(2,i),i=1,neqns)
1396  if(ldbg) write(idbg,6020) izz
1397  ! if(idbg1.ge.2) write(10,6100) neqns
1398  ! if(idbg1.ge.2) write(10,6100) (btree(1,i),btree(2,i),i=1,neqns)
1399 6000 format(i6,'(',2i6,')')
1400 6010 format(' binary tree')
1401 6020 format(' the first zero pivot is ',i4)
1402 6100 format(2i8)
1403  return
1404  end subroutine genbtq
1405 
1406  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1407 
1408  subroutine rotate(xadj,adjncy,invp,iperm,parent,btree,izz,neqns,anc,adjt,irr)
1409 
1410  implicit none
1411 
1412  integer(kind=kint), intent(in) :: xadj(:),adjncy(:),parent(:),btree(:,:)
1413  integer(kind=kint), intent(out) :: anc(:),adjt(:),invp(:),iperm(:)
1414  integer(kind=kint), intent(in) :: neqns,izz
1415  integer(kind=kint), intent(out) :: irr
1416 
1417  integer(kind=kint) :: i,j,k,l,izzz,nanc,loc,locc,ll,kk,iy
1418 
1419  !----------------------------------------------------------------------
1420  ! irr return code irr=0 node izz is not a bottom node
1421  ! irr=1 is a bottom node then rotation is
1422  ! performed
1423  !
1424  !----------------------------------------------------------------------
1425  if(izz.eq.0) then
1426  irr=0
1427  return
1428  endif
1429  izzz=invp(izz)
1430  if(btree(1,izzz).ne.0) then
1431  irr=0
1432  ! return
1433  endif
1434  irr=1
1435  !
1436  ! ancestors of izzz
1437  !
1438  nanc=0
1439  loc=izzz
1440 100 continue
1441  nanc=nanc+1
1442  anc(nanc)=loc
1443  loc=parent(loc)
1444  if(loc.ne.0) goto 100
1445  !
1446  ! to find the eligible node from ancestors of izz
1447  !
1448  ! adjt = Adj(Tree(y))
1449  l=1
1450 200 continue
1451  do i=1,neqns
1452  adjt(i)=0
1453  enddo
1454  locc=anc(l)
1455 220 continue
1456  loc=locc
1457  locc=btree(1,loc)
1458  if(locc.ne.0) goto 220
1459 230 continue
1460  do k=xadj(iperm(loc)),xadj(iperm(loc)+1)-1
1461  adjt(invp(adjncy(k)))=1
1462  enddo
1463  if(loc.ge.anc(l)) goto 250
1464  locc=btree(2,loc)
1465  if(locc.ne.0) goto 220
1466  loc=parent(loc)
1467  goto 230
1468 250 continue
1469  do ll=l+1,nanc
1470  if(adjt(anc(ll)).eq.0) then
1471  l=l+1
1472  goto 200
1473  endif
1474  enddo
1475  if(l.eq.1) goto 500
1476 
1477  !
1478  ! anc(l-1) is the eligible node
1479  !
1480  ! (1) number the node not in Ancestor(iy)
1481  iy=anc(l-1)
1482  do i=1,neqns
1483  adjt(i)=0
1484  enddo
1485  do ll=l,nanc
1486  adjt(anc(ll))=1
1487  enddo
1488  k=0
1489  do ll=1,neqns
1490  if(adjt(ll).eq.0) then
1491  k=k+1
1492  invp(iperm(ll))=k
1493  endif
1494  enddo
1495  ! (2) followed by nodes in Ancestor(iy)-Adj(T(iy))
1496 330 continue
1497  do i=1,neqns
1498  adjt(i)=0
1499  enddo
1500  locc=iy
1501 350 continue
1502  loc=locc
1503  locc=btree(1,loc)
1504  if(locc.ne.0) goto 350
1505 360 continue
1506  do kk=xadj(iperm(loc)),xadj(iperm(loc)+1)-1
1507  adjt(invp(adjncy(kk)))=1
1508  enddo
1509  if(loc.ge.iy) goto 380
1510  locc=btree(2,loc)
1511  if(locc.ne.0) goto 350
1512  loc=parent(loc)
1513  goto 360
1514 380 continue
1515  do ll=l,nanc
1516  if(adjt(anc(ll)).eq.0) then
1517  k=k+1
1518  invp(iperm(anc(ll)))=k
1519  endif
1520  enddo
1521  ! (3) and finally number the node in Adj(t(iy))
1522  do ll=l,nanc
1523  if(adjt(anc(ll)).ne.0) then
1524  k=k+1
1525  invp(iperm(anc(ll)))=k
1526  endif
1527  enddo
1528  goto 600
1529  !
1530  ! izz can be numbered last
1531  !
1532 500 continue
1533  k=0
1534  do 510 i=1,neqns
1535  if(i.eq.izzz) goto 510
1536  k=k+1
1537  invp(iperm(i))=k
1538 510 continue
1539  invp(iperm(izzz))=neqns
1540  !
1541  ! set iperm
1542  !
1543 600 continue
1544  do i=1,neqns
1545  iperm(invp(i))=i
1546  enddo
1547  if(ldbg) write(idbg,6000) (invp(i),i=1,neqns)
1548 6000 format(10i6)
1549  return
1550  end subroutine rotate
1551 
1552  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1553 
1554  subroutine bringu(zpiv,iperm,invp,parent,izz,neqns,irr)
1555 
1556  implicit none
1557 
1558  integer(kind=kint), intent(in) :: zpiv(:),parent(:)
1559  integer(kind=kint), intent(out) :: iperm(:),invp(:)
1560  integer(kind=kint), intent(in) :: neqns,izz
1561  integer(kind=kint), intent(out) :: irr
1562 
1563  integer(kind=kint) :: i,j,k,l,ib0,ib,ibp,izzp
1564 
1565  !----------------------------------------------------------------------
1566  !
1567  ! bringu brings up zero pivots from bottom of the elimination tree
1568  ! to higher nodes
1569  !
1570  ! irr = 0 complete
1571  ! = 1 impossible
1572  !
1573  ! #coded by t.arakawa
1574  !
1575  !----------------------------------------------------------------------
1576 
1577  irr=0
1578  ib0=invp(izz)
1579  ib=ib0
1580 100 continue
1581  if(ib.le.0) goto 1000
1582  ibp=parent(ib)
1583  izzp=iperm(ibp)
1584  if(zpiv(izzp).eq.0) goto 110
1585  ib=ibp
1586  goto 100
1587 110 continue
1588  invp(izz)=ibp
1589  invp(izzp)=ib0
1590  iperm(ibp)=izz
1591  iperm(ib0)=izzp
1592  if(ldbg) then
1593  do i=1,neqns
1594  if(invp(iperm(i)).ne.i) goto 210
1595  if(iperm(invp(i)).ne.i) goto 210
1596  enddo
1597  goto 220
1598 210 continue
1599  write(20,*) 'permutation error'
1600  stop
1601  endif
1602 220 continue
1603  return
1604 1000 continue
1605  irr=1
1606  end subroutine bringu
1607 
1608  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1609 
1610  subroutine posord(parent,btree,invp,iperm,pordr,nch,neqns,iw,qarent,mch)
1611 
1612  implicit none
1613 
1614  integer(kind=kint), intent(in) :: btree(:,:),qarent(:)
1615  integer(kind=kint), intent(out) :: pordr(:),invp(:),iperm(:),nch(:),iw(:),parent(:),mch(0:neqns+1)
1616  integer(kind=kint), intent(in) :: neqns
1617 
1618  integer(kind=kint) :: i,j,k,l,locc,loc,locp,invpos,ipinv,ii
1619 
1620  do i=1,neqns
1621  mch(i)=0
1622  pordr(i)=0
1623  enddo
1624  l=1
1625  locc=neqns+1
1626 10 continue
1627  loc=locc
1628  locc=btree(1,loc)
1629  if(locc.ne.0) goto 10
1630  locp=qarent(loc)
1631  mch(locp)=mch(locp)+1
1632 20 continue
1633  pordr(loc)=l
1634  if(l.ge.neqns) goto 1000
1635  l=l+1
1636  locc=btree(2,loc)
1637  if(locc.ne.0) goto 10
1638  loc=qarent(loc)
1639  locp=qarent(loc)
1640  mch(locp)=mch(locp)+mch(loc)+1
1641  goto 20
1642 1000 continue
1643  do i=1,neqns
1644  ipinv=pordr(invp(i))
1645  invp(i)=ipinv
1646  iperm(ipinv)=i
1647  iw(pordr(i))=i
1648  enddo
1649  do i=1,neqns
1650  invpos=iw(i)
1651  nch(i)=mch(invpos)
1652  ii=qarent(invpos)
1653  if(ii.gt.0.and.ii.le.neqns) then
1654  parent(i)=pordr(ii)
1655  else
1656  parent(i)=qarent(invpos)
1657  endif
1658  enddo
1659  if(ldbg) write(idbg,6020)
1660  if(ldbg) write(idbg,6000) (pordr(i),i=1,neqns)
1661  if(ldbg) write(idbg,6030)
1662  if(ldbg) write(idbg,6050)
1663  if(ldbg) write(idbg,6000) (parent(i),i=1,neqns)
1664  if(ldbg) write(idbg,6000) (invp(i),i=1,neqns)
1665  if(ldbg) write(idbg,6040)
1666  if(ldbg) write(idbg,6000) (iperm(i),i=1,neqns)
1667  if(ldbg) write(idbg,6010)
1668  if(ldbg) write(idbg,6000) (nch(i),i=1,neqns)
1669 6000 format(10i6)
1670 6010 format(' nch')
1671 6020 format(' post order')
1672 6030 format(/' invp ')
1673 6040 format(/' iperm ')
1674 6050 format(/' parent')
1675  return
1676  end subroutine posord
1677 
1678  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1679 
1680  subroutine gnleaf(xadj,adjncy,invp,iperm,pordr,nch,adjncp,xleaf,leaf,neqns,lnleaf)
1681 
1682  implicit none
1683 
1684  integer(kind=kint), intent(in) :: xadj(:),adjncy(:),pordr(:),nch(:),invp(:),iperm(:)
1685  integer(kind=kint), intent(out) :: xleaf(:),leaf(:),adjncp(:)
1686  integer(kind=kint), intent(in) :: neqns
1687 
1688  integer(kind=kint) i,j,k,l,m,n,ik,istart,ip,iq,lnleaf,lc1,lc
1689 
1690  l=1
1691  ik=0
1692  istart=0
1693  do 100 i=1,neqns
1694  xleaf(i)=l
1695  ip=iperm(i)
1696  do k=xadj(ip),xadj(ip+1)-1
1697  iq=invp(adjncy(k))
1698  if(iq.lt.i) then
1699  ik=ik+1
1700  adjncp(ik)=iq
1701  endif
1702  enddo
1703  m=ik-istart
1704  if(m.eq.0) goto 131
1705  call qqsort(adjncp(istart+1:),m)
1706  lc1=adjncp(istart+1)
1707  if(lc1.ge.i) goto 100
1708  leaf(l)=lc1
1709  l=l+1
1710  do k=istart+2,ik
1711  lc=adjncp(k)
1712  ! if(lc.ge.i) goto 125
1713  if(lc1.lt.lc-nch(lc)) then
1714  leaf(l)=lc
1715  l=l+1
1716  endif
1717 125 continue
1718  lc1=lc
1719  enddo
1720  ik=1
1721  istart=ik
1722 131 continue
1723 100 continue
1724  xleaf(neqns+1)=l
1725  lnleaf=l-1
1726  if(ldbg) write(idbg,6020)
1727  if(ldbg) write(idbg,6000) (xleaf(i),i=1,neqns+1)
1728  if(ldbg) write(idbg,6010) lnleaf
1729  if(ldbg) write(idbg,6000) (leaf(i),i=1,lnleaf)
1730  return
1731 6000 format(10i6)
1732 6010 format(' leaf (len = ',i6,')')
1733 6020 format(' xleaf')
1734  end subroutine gnleaf
1735 
1736  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1737 
1738  subroutine countclno(parent,xleaf,leaf,neqns,nstop,lncol,ir)
1739 
1740  implicit none
1741 
1742  ! Count total number of non-zero elements
1743  ! which include fill-in.
1744  ! A and C region of given sparse matrix will consider.
1745  ! D region will not consider because of D is treat as
1746  ! dens matrix.
1747  !
1748  integer(kind=kint), intent(in) :: parent(:),xleaf(:),leaf(:)
1749  integer(kind=kint), intent(in) :: neqns, nstop
1750  integer(kind=kint), intent(out) :: lncol, ir
1751 
1752  integer(kind=kint) :: i,j,k,l,nc,ks,ke,nxleaf
1753 
1754  nc=0
1755  ir=0
1756  l=1
1757  do 100 i=1,neqns
1758  ks=xleaf(i)
1759  ke=xleaf(i+1)-1
1760  if(ke.lt.ks) goto 100
1761  nxleaf=leaf(ks)
1762  do 110 k=ks,ke-1
1763  j=nxleaf
1764  nxleaf=leaf(k+1)
1765 105 continue
1766  if(j.ge.nxleaf) goto 110
1767  if(j.ge.nstop) then
1768  goto 100
1769  endif
1770  l=l+1
1771  j=parent(j)
1772  goto 105
1773 110 continue
1774  j=leaf(ke)
1775 115 continue
1776  if(j.ge.nstop) goto 100
1777  if(j.ge.i.or.j.eq.0) goto 100
1778  l=l+1
1779  j=parent(j)
1780  goto 115
1781 100 continue
1782  lncol=l-1
1783  return
1784  end subroutine countclno
1785 
1786  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1787 
1788  subroutine gnclno(parent,pordr,xleaf,leaf,xlnzr,colno,neqns,nstop,lncol,ir)
1789 
1790  implicit none
1791 
1792  integer(kind=kint), intent(in) :: parent(:),pordr(:),xleaf(:),leaf(:)
1793  integer(kind=kint), intent(out) :: colno(:),xlnzr(:)
1794  integer(kind=kint), intent(in) :: neqns, nstop
1795  integer(kind=kint), intent(out) :: lncol,ir
1796 
1797  integer(kind=kint) :: i,j,k,l,nc,ks,ke,nxleaf
1798 
1799  nc=0
1800  ir=0
1801  l=1
1802  do 100 i=1,neqns
1803  xlnzr(i)=l
1804  ks=xleaf(i)
1805  ke=xleaf(i+1)-1
1806  if(ke.lt.ks) goto 100
1807  nxleaf=leaf(ks)
1808  do 110 k=ks,ke-1
1809  j=nxleaf
1810  nxleaf=leaf(k+1)
1811 105 continue
1812  if(j.ge.nxleaf) goto 110
1813  if(j.ge.nstop) then
1814  goto 100
1815  endif
1816  colno(l)=j
1817  l=l+1
1818  j=parent(j)
1819  goto 105
1820 110 continue
1821  j=leaf(ke)
1822 115 continue
1823  if(j.ge.nstop) goto 100
1824  if(j.ge.i.or.j.eq.0) goto 100
1825  colno(l)=j
1826  l=l+1
1827  j=parent(j)
1828  goto 115
1829 100 continue
1830  xlnzr(neqns+1)=l
1831  lncol=l-1
1832  if(ldbg) write(idbg,6010)
1833  ! if(idbg1.ne.0) write(6,6000) (xlnzr(i),i=1,neqns+1)
1834  if(ldbg) write(idbg,6020) lncol
1835  if(ldbg) then
1836  do k=1,neqns
1837  write(idbg,6100) k
1838  write(idbg,6000) (colno(i),i=xlnzr(k),xlnzr(k+1)-1)
1839  enddo
1840  endif
1841 6000 format(10i4)
1842 6010 format(' xlnzr')
1843 6020 format(' colno (lncol =',i10,')')
1844 6100 format(/' row = ',i6)
1845  return
1846  end subroutine gnclno
1847 
1848  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1849 
1850  subroutine qmdrch(root,xadj,adjncy,deg,marker,rchsze,rchset,nhdsze,nbrhd)
1851 
1852  implicit none
1853 
1854  integer(kind=kint), intent(in) :: deg(:),xadj(:),adjncy(:)
1855  integer(kind=kint), intent(out) :: rchset(:),marker(:),nbrhd(:)
1856  integer(kind=kint), intent(in) :: root
1857  integer(kind=kint), intent(out) :: nhdsze,rchsze
1858 
1859  integer(kind=kint) :: i,j,k,l, istrt, istop, jstrt, jstop, nabor, node
1860 
1861  nhdsze=0
1862  rchsze=0
1863  istrt=xadj(root)
1864  istop=xadj(root+1)-1
1865  if(istop.lt.istrt) return
1866  do 600 i=istrt,istop
1867  nabor=adjncy(i)
1868  if(nabor.eq.0) return
1869  if(marker(nabor).ne.0) goto 600
1870  if(deg(nabor).lt.0) goto 200
1871  rchsze=rchsze+1
1872  rchset(rchsze)=nabor
1873  marker(nabor)=1
1874  goto 600
1875 200 marker(nabor)=-1
1876  nhdsze=nhdsze+1
1877  nbrhd(nhdsze)=nabor
1878 300 jstrt=xadj(nabor)
1879  jstop=xadj(nabor+1)-1
1880  do 500 j=jstrt,jstop
1881  node=adjncy(j)
1882  nabor=-node
1883  !if(node) 300,600,400
1884  if(node<0) then
1885  goto 300
1886  elseif(node==0) then
1887  goto 600
1888  endif
1889 400 if(marker(node).ne.0) goto 500
1890  rchsze=rchsze+1
1891  rchset(rchsze)=node
1892  marker(node)=1
1893 500 continue
1894 600 continue
1895  return
1896  end subroutine qmdrch
1897 
1898  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1899 
1900  subroutine qmdupd(xadj,adjncy,nlist,list,deg,qsize,qlink,marker,rchset,nbrhd)
1901 
1902  implicit none
1903 
1904  integer(kind=kint), intent(in) :: adjncy(:),list(:),xadj(:)
1905  integer(kind=kint), intent(out) :: marker(:),nbrhd(:),rchset(:),deg(:),qsize(:),qlink(:)
1906  integer(kind=kint), intent(in) :: nlist
1907 
1908  integer(kind=kint) :: i,j,k,l, deg0,deg1,il,inhd,inode,irch,jstrt,jstop,mark,nabor,nhdsze,node,rchsze
1909 
1910  if(nlist.le.0) return
1911  deg0=0
1912  nhdsze=0
1913  do il=1,nlist
1914  node=list(il)
1915  deg0=deg0+qsize(node)
1916  jstrt=xadj(node)
1917  jstop=xadj(node+1)-1
1918 
1919  do 100 j=jstrt,jstop
1920  nabor=adjncy(j)
1921  if(marker(nabor).ne.0.or.deg(nabor).ge.0) goto 100
1922  marker(nabor)=-1
1923  nhdsze=nhdsze+1
1924  nbrhd(nhdsze)=nabor
1925 100 continue
1926  enddo
1927 
1928  if(nhdsze.gt.0) call qmdmrg(xadj,adjncy,deg,qsize,qlink,marker,deg0,nhdsze,nbrhd,rchset,nbrhd(nhdsze+1:))
1929  do 600 il=1,nlist
1930  node=list(il)
1931  mark=marker(node)
1932  if(mark.gt.1.or.mark.lt.0) goto 600
1933  call qmdrch(node,xadj,adjncy,deg,marker,rchsze,rchset,nhdsze,nbrhd)
1934  deg1=deg0
1935  if(rchsze.le.0) goto 400
1936  do irch=1,rchsze
1937  inode=rchset(irch)
1938  deg1=deg1+qsize(inode)
1939  marker(inode)=0
1940  enddo
1941 400 deg(node)=deg1-1
1942  if(nhdsze.le.0) goto 600
1943  do inhd=1,nhdsze
1944  inode=nbrhd(inhd)
1945  marker(inode)=0
1946  enddo
1947 600 continue
1948  return
1949  end subroutine qmdupd
1950 
1951  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1952 
1953  subroutine qmdot(root,xadj,adjncy,marker,rchsze,rchset,nbrhd)
1954 
1955  implicit none
1956 
1957  integer(kind=kint), intent(in) :: marker(:),rchset(:),nbrhd(:),xadj(:)
1958  integer(kind=kint), intent(out) :: adjncy(:)
1959  integer(kind=kint), intent(in) :: rchsze,root
1960 
1961  integer(kind=kint) :: i,j,k,l,irch,inhd,node,jstrt,jstop,link,nabor
1962 
1963  irch=0
1964  inhd=0
1965  node=root
1966 100 jstrt=xadj(node)
1967  jstop=xadj(node+1)-2
1968  if(jstop.lt.jstrt) goto 300
1969  do j=jstrt,jstop
1970  irch=irch+1
1971  adjncy(j)=rchset(irch)
1972  if(irch.ge.rchsze) goto 400
1973  enddo
1974 300 link=adjncy(jstop+1)
1975  node=-link
1976  if(link.lt.0) goto 100
1977  inhd=inhd+1
1978  node=nbrhd(inhd)
1979  adjncy(jstop+1)=-node
1980  goto 100
1981 400 adjncy(j+1)=0
1982  do 600 irch=1,rchsze
1983  node=rchset(irch)
1984  if(marker(node).lt.0) goto 600
1985  jstrt=xadj(node)
1986  jstop=xadj(node+1)-1
1987  do 500 j=jstrt,jstop
1988  nabor=adjncy(j)
1989  if(marker(nabor).ge.0) goto 500
1990  adjncy(j)=root
1991  goto 600
1992 500 continue
1993 600 continue
1994  return
1995  end subroutine qmdot
1996 
1997  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1998 
1999  subroutine qmdmrg(xadj,adjncy,deg,qsize,qlink,marker,deg0,nhdsze,nbrhd,rchset,ovrlp)
2000 
2001  implicit none
2002 
2003  integer(kind=kint), intent(in) :: adjncy(:),nbrhd(:),xadj(:)
2004  integer(kind=kint), intent(out) :: deg(:),marker(:),rchset(:),ovrlp(:),qsize(:),qlink(:)
2005  integer(kind=kint), intent(in) :: nhdsze
2006 
2007  integer(kind=kint) :: i,j,k,l, deg0,deg1,head,inhd,iov,irch,jstrt,jstop,link,lnode,mark,mrgsze,nabor,node,novrlp,rchsze,root
2008 
2009 
2010  if(nhdsze.le.0) return
2011  do inhd=1,nhdsze
2012  root=nbrhd(inhd)
2013  marker(root)=0
2014  enddo
2015  do 1400 inhd=1,nhdsze
2016  root=nbrhd(inhd)
2017  marker(root)=-1
2018  rchsze=0
2019  novrlp=0
2020  deg1=0
2021 200 jstrt=xadj(root)
2022  jstop=xadj(root+1)-1
2023  do 600 j=jstrt,jstop
2024  nabor=adjncy(j)
2025  root=-nabor
2026  !if(nabor) 200,700,300
2027  if(nabor<0) then
2028  goto 200
2029  elseif(nabor==0) then
2030  goto 700
2031  endif
2032 300 mark=marker(nabor)
2033  !if(mark)600,400,500
2034  if(mark<0) then
2035  goto 600
2036  elseif(mark>0) then
2037  goto 500
2038  endif
2039 400 rchsze=rchsze+1
2040  rchset(rchsze)=nabor
2041  deg1=deg1+qsize(nabor)
2042  marker(nabor)=1
2043  goto 600
2044 500 if(mark.gt.1) goto 600
2045  novrlp=novrlp+1
2046  ovrlp(novrlp)=nabor
2047  marker(nabor)=2
2048 600 continue
2049 700 head=0
2050  mrgsze=0
2051  do 1100 iov=1,novrlp
2052  node=ovrlp(iov)
2053  jstrt=xadj(node)
2054  jstop=xadj(node+1)-1
2055  do 800 j=jstrt,jstop
2056  nabor=adjncy(j)
2057  if(marker(nabor).ne.0) goto 800
2058  marker(node)=1
2059  goto 1100
2060 800 continue
2061  mrgsze=mrgsze+qsize(node)
2062  marker(node)=-1
2063  lnode=node
2064 900 link=qlink(lnode)
2065  if(link.le.0) goto 1000
2066  lnode=link
2067  goto 900
2068 1000 qlink(lnode)=head
2069  head=node
2070 1100 continue
2071  if(head.le.0) goto 1200
2072  qsize(head)=mrgsze
2073  deg(head)=deg0+deg1-1
2074  marker(head)=2
2075 1200 root=nbrhd(inhd)
2076  marker(root)=0
2077  if(rchsze.le.0) goto 1400
2078  do irch=1,rchsze
2079  node=rchset(irch)
2080  marker(node)=0
2081  enddo
2082 1400 continue
2083  return
2084  end subroutine qmdmrg
2085 
2086  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2087 
2088  subroutine ldudecomposec(xlnzr_a,colno_a,invp_a,iperm_a,ndeg,nttbr_c,irow_c,jcol_c,ncol,nrow,xlnzr_c,colno_c,lncol_c)
2089  ! find fill-in position in C which placed under A and set it in xlnzr_c, colno_c
2090 
2091  use m_crs_matrix_lag
2092  implicit none
2093  !input
2094  integer(kind=kint), intent(in) :: xlnzr_a(:)
2095  integer(kind=kint), intent(in) :: colno_a(:)
2096  integer(kind=kint), intent(in) :: iperm_a(:)
2097  integer(kind=kint), intent(in) :: invp_a(:)
2098  integer(kind=kint), intent(in) :: ndeg
2099  integer(kind=kint), intent(in) :: nttbr_c
2100  integer(kind=kint), intent(in) :: irow_c(:)
2101  integer(kind=kint), intent(inout) :: jcol_c(:)
2102  integer(kind=kint), intent(in) :: ncol
2103  integer(kind=kint), intent(in) :: nrow
2104 
2105  !output
2106  integer(kind=kint), pointer :: xlnzr_c(:)
2107  integer(kind=kint), pointer :: colno_c(:)
2108  integer(kind=kint), intent(out) :: lncol_c
2109 
2110  ! internal
2111  integer(kind=kint) :: i,j,k,l,m,n
2112  integer(kind=kint) :: ks, ke, ipass, ierr
2113  logical, allocatable :: cnz(:)
2114  type(crs_matrix) :: crs_c
2115 
2116  !permtate column in C for crs_c
2117  do i=1,nttbr_c
2118  jcol_c(i)=invp_a(jcol_c(i))
2119  end do
2120 
2121  ! make Compact Column Storoge using symbolic information.
2122  call symbolicirjctocrs(ndeg, nttbr_c, irow_c, jcol_c, ncol, nrow, crs_c)
2123 
2124  ! symbolic LDU factorization for C matrix
2125  allocate(cnz(ncol), stat=ierr)
2126  if(ierr .ne. 0) then
2127  call errtrp('stop due to allocation error.')
2128  end if
2129  do ipass = 1,2
2130  lncol_c = 0
2131  do k=1,nrow
2132  ! set cnz as non-zero pattern of C
2133  cnz = .false.
2134  ks = crs_c%ia(k)
2135  ke = crs_c%ia(k+1)-1
2136  if (ke .lt. ks) then
2137  if (ipass .eq. 2) then
2138  xlnzr_c(k+1)=lncol_c+1
2139  end if
2140  cycle ! in case of zero vector, no need to check dot product. not cycle?
2141  end if
2142 
2143  do i=ks,ke
2144  cnz(crs_c%ja(i)) = .true.
2145  end do
2146 
2147  ! check for non-zero dot product and update cnz for each point of cnz
2148  do i=2,ncol
2149  ks = xlnzr_a(i)
2150  ke = xlnzr_a(i+1)-1
2151  if (ke .lt. ks) then ! in case of column of A is zero vector.
2152  cycle
2153  end if
2154  do j=ks,ke
2155  if (cnz(colno_a(j))) then
2156  cnz(i) = .true.
2157  exit
2158  end if
2159  end do
2160  end do
2161 
2162  do i=1,ncol
2163  if (cnz(i)) then
2164  lncol_c = lncol_c + 1
2165  if (ipass .eq. 2) then
2166  colno_c(lncol_c) = i
2167  end if
2168  end if
2169  end do
2170  if (ipass .eq. 2) then
2171  xlnzr_c(k+1)=lncol_c + 1
2172  end if
2173  end do
2174 
2175  if (ipass .eq. 1) then
2176  allocate(xlnzr_c(nrow+1),colno_c(lncol_c), stat=ierr)
2177  if(ierr .ne. 0) then
2178  call errtrp('stop due to allocation error.')
2179  end if
2180  xlnzr_c(1)=1
2181  end if
2182  end do
2183 
2184  ! restore order of C column.
2185  do i=1,nttbr_c
2186  jcol_c(i)=iperm_a(jcol_c(i))
2187  end do
2188 
2189  return
2190 
2191  end subroutine ldudecomposec
2192 
2193  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2194 
2195  subroutine qqsort(iw,ik)
2196 
2197  implicit none
2198 
2199  integer(kind=kint), intent(out) :: iw(:)
2200  integer(kind=kint), intent(in) :: ik
2201 
2202  integer(kind=kint) :: l,m,itemp
2203 
2204  !----------------------------------------------------------------------
2205  ! sort in increasing order up to i
2206  !
2207  ! iw array
2208  ! ik number of input/output
2209  ! i deal with numbers less than this numberi
2210  !
2211  !----------------------------------------------------------------------
2212 
2213  if(ik.le.1) return
2214  do l=1,ik-1
2215  do 110 m=l+1,ik
2216  if(iw(l).lt.iw(m)) goto 110
2217  itemp=iw(l)
2218  iw(l)=iw(m)
2219  iw(m)=itemp
2220 110 continue
2221  enddo
2222 200 continue
2223  return
2224  end subroutine qqsort
2225 
2226 
2227  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2228 
2229  subroutine staij1(isw,i,j,aij,dsi,ir)
2230 
2231  implicit none
2232 
2233  !----------------------------------------------------------------------
2234  !
2235  ! this routine sets an non-zero entry of the matrix.
2236  ! (symmetric version)
2237  !
2238  ! (i)
2239  ! isw =0 set the value
2240  ! =1 add the value
2241  ! i row entry
2242  ! j column entry
2243  ! aij value
2244  !
2245  ! (o)
2246  ! iv communication array
2247  !
2248  ! #coded by t.arakawa
2249  !
2250  !----------------------------------------------------------------------
2251  !
2252  type(dsinfo) :: dsi
2253  real(kind=kreal), intent(out) :: aij(:) ! ndeg*ndeg
2254  integer(kind=kint), intent(in) :: isw, i, j
2255  integer(kind=kint), intent(out) :: ir
2256 
2257  integer(kind=kint) :: ndeg, neqns, nstop, ndeg2, ndeg2l, ierr
2258  ndeg=dsi%ndeg
2259  neqns=dsi%neqns
2260  nstop=dsi%nstop
2261  ndeg2=ndeg*ndeg
2262  ndeg2l=ndeg*(ndeg+1)/2
2263 
2264  ir=0
2265  ierr=0
2266 
2267 
2268  ! array allocation
2269  if(dsi%stage.ne.20) then
2270  if(dsi%stage.eq.30) write(ilog,*) 'Warning a matrix was build up but never solved.'
2271  !
2272  ! for diagonal
2273  !
2274  allocate(dsi%diag(ndeg2l,neqns), stat=ierr)
2275  if(ierr .ne. 0) then
2276  call errtrp('stop due to allocation error.')
2277  end if
2278  dsi%diag=0
2279  !
2280  ! for lower triangle
2281  !
2282  allocate(dsi%zln(ndeg2,dsi%lncol), stat=ierr)
2283  if(ierr .ne. 0) then
2284  call errtrp('stop due to allocation error.')
2285  end if
2286  dsi%zln=0
2287  !
2288  ! for dense window !TODO delete this and corresponding line in addr3()
2289  !
2290  ! allocate(dsi%dsln(ndeg2,dsi%lndsln))! because of there is no dense window
2291  ! dsi%dsln=0
2292 
2293  dsi%stage=20
2294  endif
2295 
2296  ! set value
2297  if(ndeg.le.2) then
2298  call addr0(isw,i,j,aij,dsi%invp,dsi%xlnzr,dsi%colno,dsi%diag,dsi%zln,dsi%dsln,nstop,dsi%ndeg,ir)
2299  elseif(ndeg.eq.3) then
2300  write(idbg,*) 'ndeg=1 only'
2301  stop
2302  else
2303  write(idbg,*) 'ndeg=1 only'
2304  stop
2305  endif
2306 1000 continue
2307  return
2308  end subroutine staij1
2309 
2310  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2311  !
2312  ! After here, routines specialized for ndeg = 1
2313  !
2314  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2315 
2316  ! LDU decompose of A (1..nstop-1) region
2317  subroutine sum(ic,xlnzr,colno,zln,diag,nch,par,neqns)
2318 
2319  implicit none
2320 
2321  integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
2322  integer(kind=kint), intent(in) :: ic, neqns
2323  real(kind=kreal), intent(inout) :: zln(:),diag(:)
2324  integer(kind=kint), intent(out) :: nch(:)
2325 
2326  real(kind=kreal) :: s, t, zz, piv
2327  integer(kind=kint) :: ks, ke, kk, k, jc, jj, j, ierr
2328  integer(kind=kint) :: isem
2329  real(kind=kreal),allocatable :: temp(:)
2330  integer(kind=kint),allocatable :: indx(:)
2331  allocate(temp(neqns),indx(neqns), stat=ierr)
2332  if(ierr .ne. 0) then
2333  call errtrp('stop due to allocation error.')
2334  end if
2335 
2336 2 continue
2337  ks=xlnzr(ic)
2338  ke=xlnzr(ic+1)
2339  t=0.0d0
2340  ! do 100 i=1,ic
2341  ! temp(i)=0.0d0
2342  ! 100 continue
2343  do k=ks,ke-1
2344  jc=colno(k)
2345  indx(jc)=ic
2346  s=0.0d0
2347  do jj=xlnzr(jc),xlnzr(jc+1)-1
2348  j=colno(jj)
2349  if(indx(j).eq.ic) then
2350  s=s+temp(j)*zln(jj)
2351  endif
2352  enddo
2353  ! j1=xlnzr(jc)
2354  ! jj=xlnzr(jc+1)-j1
2355  ! ss=ddoti(jj,zln(j1),colno(j1),temp)
2356  ! zz=zln(k)-ddoti(jj,zln(j1),colno(j1),temp)
2357  zz=zln(k)-s
2358  zln(k)=zz*diag(jc)
2359  temp(jc)=zz
2360  t=t+zz*zln(k)
2361  enddo
2362  piv=diag(ic)-t
2363  if(dabs(piv).gt.rmin) then
2364  diag(ic)=1.0d0/piv
2365  endif
2366 1 continue
2367  ! if(isem.eq.1) then !DBG
2368  isem=0
2369  nch(ic)=-1
2370  kk=par(ic)
2371  nch(kk)=nch(kk)-1
2372  isem=1
2373  ! else !DBG
2374  ! goto 1 !DBG
2375  ! endi !DBG
2376  return
2377  end subroutine sum
2378 
2379  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2380 
2381  ! LDU decompose of C (nstop..neqnsA+neqnsd) region
2382  subroutine sum1(ic,xlnzr,colno,zln,diag,par,neqns)
2383 
2384  implicit none
2385 
2386  integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
2387  integer(kind=kint), intent(in) :: ic, neqns
2388  real(kind=kreal), intent(inout) :: zln(:),diag(:)
2389 
2390  real(kind=kreal) :: s, t, zz
2391  integer(kind=kint) :: ks, ke, k, jc, j, jj, ierr
2392  real(kind=kreal),allocatable :: temp(:)
2393  integer(kind=kint),allocatable :: indx(:)
2394  integer(kind=kint) :: i
2395 
2396  ierr=0
2397 
2398  allocate(temp(neqns),indx(neqns), stat=ierr)
2399  if(ierr .ne. 0) then
2400  call errtrp('stop due to allocation error.')
2401  end if
2402 
2403  do i=1,neqns
2404  temp(i)=0
2405  end do
2406 
2407  ks=xlnzr(ic)
2408  ke=xlnzr(ic+1)
2409  t=0.0d0
2410  ! do 100 i=1,ic
2411  ! temp(i)=0.0d0
2412  ! 100 continue
2413  do k=ks,ke-1
2414  jc=colno(k)
2415  indx(jc)=ic
2416  s=0.0d0
2417  do jj=xlnzr(jc),xlnzr(jc+1)-1
2418  j=colno(jj)
2419  if(indx(j).eq.ic) then
2420  s=s+temp(j)*zln(jj)
2421  endif
2422  enddo
2423  zz=zln(k)-s
2424  ! j1=xlnzr(jc)
2425  ! jj=xlnzr(jc+1)-j1
2426  ! zz=zln(k)-ddoti(jj,zln(j1),colno(j1),temp)
2427  zln(k)=zz
2428  temp(jc)=zz
2429  ! t=t+zz*zz*diag(jc)
2430  enddo
2431  ! diag(ic)=diag(ic)-t
2432  return
2433  end subroutine sum1
2434 
2435  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2436 
2437  ! LDU decompose and Update D region.
2438  subroutine sum2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln)
2439 
2440  implicit none
2441 
2442  integer(kind=kint), intent(in) :: neqns, nstop
2443  integer(kind=kint), intent(in) :: xlnzr(:),colno(:)
2444  real(kind=kreal), intent(inout) :: zln(:),diag(:)
2445  integer(kind=kint), pointer :: spdslnidx(:)
2446  real(kind=kreal), pointer :: spdslnval(:,:)
2447  integer(kind=kint), intent(out) :: nspdsln
2448 
2449  real(kind=kreal) :: s, t
2450  integer(kind=kint) :: ks, ke, kk, k, jc, jj, j, j1,j2
2451  integer(kind=kint) :: ic, i, loc, ierr
2452  integer(kind=kint) :: ispdsln
2453  logical :: ftflag
2454  real(kind=kreal),allocatable :: temp(:)
2455  integer(kind=kint),allocatable :: indx(:)
2456  ierr=0
2457  allocate(temp(neqns),indx(neqns), stat=ierr)
2458  if(ierr .ne. 0) then
2459  call errtrp('stop due to allocation error.')
2460  end if
2461  temp=0
2462 
2463  nspdsln=0
2464  do ic=nstop,neqns
2465  ks=xlnzr(ic)
2466  ke=xlnzr(ic+1)-1
2467  do k=ks,ke
2468  jj=colno(k)
2469  indx(jj)=ic
2470  end do
2471  do jc=nstop,ic-1
2472  j1=xlnzr(jc)
2473  j2=xlnzr(jc+1)
2474  do jj=xlnzr(jc),xlnzr(jc+1)-1
2475  j=colno(jj)
2476  if(indx(j).eq.ic) then
2477  nspdsln=nspdsln+1
2478  exit
2479  endif
2480  end do
2481  end do
2482  end do
2483  allocate(spdslnidx(nspdsln),spdslnval(1,nspdsln), stat=ierr)
2484  if(ierr .ne. 0) then
2485  call errtrp('stop due to allocation error.')
2486  end if
2487 
2488  loc=0
2489  ispdsln=0
2490  spdslnval=0
2491  ftflag = .true.
2492  do ic=nstop,neqns
2493  do i=1,nstop
2494  temp(i)=0.0d0
2495  enddo
2496  ks=xlnzr(ic)
2497  ke=xlnzr(ic+1)-1
2498  do k=ks,ke
2499  jj=colno(k)
2500  temp(jj)=zln(k)
2501  zln(k)=temp(jj)*diag(jj)
2502  indx(jj)=ic
2503  diag(ic)=diag(ic)-temp(jj)*zln(k)
2504  enddo
2505  do jc=nstop,ic-1
2506  loc=loc+1
2507  do jj=xlnzr(jc),xlnzr(jc+1)-1
2508  j=colno(jj)
2509  if(indx(j).eq.ic) then
2510  if (ftflag) then
2511  ispdsln=ispdsln+1
2512  ftflag=.false.
2513  end if
2514  spdslnidx(ispdsln)=loc
2515  spdslnval(1,ispdsln)=spdslnval(1,ispdsln)-temp(j)*zln(jj)
2516  endif
2517  enddo
2518  ftflag = .true.
2519  ! j1=xlnzr(jc)
2520  ! jj=xlnzr(jc+1)-j1
2521  ! dsln(loc)=dsln(loc)-ddoti(jj,zln(j1),colno(j1),temp)
2522  enddo
2523  enddo
2524  return
2525  end subroutine sum2_child
2526 
2527  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2528 
2529  subroutine sum3(n,dsln,diag)
2530 
2531  implicit none
2532 
2533  real(kind=kreal), intent(inout) :: dsln(:),diag(:)
2534  integer(kind=kint), intent(in) :: n
2535 
2536  integer(kind=kint) :: i, j, loc, ierr
2537  real(kind=kreal),allocatable :: temp(:)
2538  integer(kind=kint),allocatable :: indx(:)
2539  allocate(temp(n),indx(n), stat=ierr)
2540  if(ierr .ne. 0) then
2541  call errtrp('stop due to allocation error.')
2542  end if
2543 
2544  if(n.le.0) goto 1000
2545  indx(1)=0
2546  loc=1
2547  diag(1)=1.0d0/diag(1)
2548  do i=2,n
2549  indx(i)=loc
2550  do j=1,i-1
2551  dsln(loc)=dsln(loc)-dot_product(dsln(indx(i):indx(i)+j-2),dsln(indx(j):indx(j)+j-2))
2552  loc=loc+1
2553  enddo
2554  temp(1:i-1)=dsln(indx(i):indx(i)+i-2)*diag(1:i-1)
2555  diag(i)=diag(i)-dot_product(temp(1:i-1),dsln(indx(i):indx(i)+i-2))
2556  dsln(indx(i):indx(i)+i-2)=temp(1:i-1)
2557  diag(i)=1.0d0/diag(i)
2558  enddo
2559 1000 continue
2560 
2561  return
2562  end subroutine sum3
2563 
2564  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2565 
2566  real(kind=kreal) function spdot2(b,zln,colno,ks,ke)
2567 
2568  implicit none
2569 
2570  integer(kind=kint), intent(in) :: colno(:)
2571  integer(kind=kint), intent(in) :: ks,ke
2572  real(kind=kreal), intent(in) :: zln(:),b(:)
2573 
2574  integer(kind=kint) :: j,jj
2575  real(kind=kreal) :: s
2576 
2577  !----------------------------------------------------------------------
2578  !
2579  ! spdot1 performs inner product of sparse vectors
2580  !
2581  !
2582  ! #coded by t.arakawa
2583  !
2584  !----------------------------------------------------------------------
2585  !
2586  s=0.0d0
2587  do jj=ks,ke
2588  j=colno(jj)
2589  s=s+zln(jj)*b(j)
2590  enddo
2591  spdot2=s
2592  end function spdot2
2593 
2594  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2595 
2596  real(kind=kreal) function ddot(a,b,n)
2597 
2598  implicit none
2599 
2600  real(kind=kreal), intent(in) :: a(n),b(n)
2601  integer(kind=kint), intent(in) :: n
2602 
2603  real(kind=kreal) :: s
2604  integer(kind=kint) :: i
2605 
2606  s=0.0d0
2607  do i=1,n
2608  s=s+a(i)*b(i)
2609  enddo
2610  ddot=s
2611  return
2612  end function ddot
2613 
2614  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2615 
2616  subroutine addr0(isw,i,j,aij,invp,xlnzr,colno,diag,zln,dsln,nstop,ndeg,ir)
2617 
2618  implicit none
2619 
2620  integer(kind=kint), intent(in) :: isw ! 0: renew diag, dsln, zln other: add to diag, dsln, zln
2621  integer(kind=kint), intent(in) :: i,j,nstop, ndeg, invp(:),xlnzr(:),colno(:)
2622  real(kind=kreal), intent(inout) :: zln(:,:),diag(:,:),dsln(:,:),aij(:)
2623  integer(kind=kint), intent(out) :: ir
2624 
2625  integer(kind=kint) :: ndeg2, ii, jj, itrans, k, i0, j0, l, ks, ke
2626  integer(kind=kint), parameter :: idbg=0
2627  ndeg2=ndeg*ndeg
2628 
2629  ir=0
2630  ii=invp(i)
2631  jj=invp(j)
2632  if(idbg.ne.0) write(idbg,*) 'addr0',ii,jj,aij
2633  if(ii.eq.jj) then
2634  if(ndeg2.eq.1) then
2635  if(isw.eq.0) then
2636  diag(1,ii)=aij(1)
2637  else
2638  diag(1,ii)=diag(1,ii)+aij(1)
2639  endif
2640  elseif(ndeg2.eq.4) then
2641  if(isw.eq.0) then
2642  diag(1,ii)=aij(1)
2643  diag(2,ii)=aij(2)
2644  diag(3,ii)=aij(4)
2645  else
2646  diag(1,ii)=diag(1,ii)+aij(1)
2647  diag(2,ii)=diag(2,ii)+aij(2)
2648  diag(3,ii)=diag(3,ii)+aij(4)
2649  endif
2650  endif
2651  goto 1000
2652  endif
2653  itrans=0
2654  if(jj.gt.ii) then
2655  k=jj
2656  jj=ii
2657  ii=k
2658  itrans=1
2659  endif
2660  if(jj.ge.nstop) then
2661  i0=ii-nstop
2662  j0=jj-nstop+1
2663  k=i0*(i0-1)/2+j0
2664  if(ndeg2.eq.1) then
2665  dsln(1,k)=aij(1)
2666  goto 1000
2667  elseif(ndeg2.eq.4) then
2668  if(itrans.eq.0) then
2669  do l=1,ndeg2
2670  dsln(l,k)=aij(l)
2671  enddo
2672  goto 1000
2673  else
2674  dsln(1,k)=aij(1)
2675  dsln(2,k)=aij(3)
2676  dsln(3,k)=aij(2)
2677  dsln(4,k)=aij(4)
2678  goto 1000
2679  endif
2680  endif
2681  endif
2682  ks=xlnzr(ii)
2683  ke=xlnzr(ii+1)-1
2684  do k=ks,ke
2685  if(colno(k).eq.jj) then
2686  if(isw.eq.0) then
2687  if(ndeg2.eq.1) then
2688  zln(1,k)=aij(1)
2689  elseif(ndeg2.eq.4) then
2690  if(itrans.eq.0) then
2691  do l=1,ndeg2
2692  zln(l,k)=aij(l)
2693  enddo
2694  else
2695  zln(1,k)=aij(1)
2696  zln(2,k)=aij(3)
2697  zln(3,k)=aij(2)
2698  zln(4,k)=aij(4)
2699  endif
2700  endif
2701  else
2702  if(ndeg2.eq.1) then
2703  zln(1,k)=zln(1,k)+aij(1)
2704  elseif(ndeg2.eq.4) then
2705  if(itrans.eq.0) then
2706  do l=1,ndeg2
2707  zln(l,k)=zln(l,k)+aij(l)
2708  enddo
2709  else
2710  zln(1,k)=zln(1,k)+aij(1)
2711  zln(2,k)=zln(2,k)+aij(3)
2712  zln(3,k)=zln(3,k)+aij(2)
2713  zln(4,k)=zln(4,k)+aij(4)
2714  endif
2715  endif
2716  endif
2717  goto 1000
2718  endif
2719  enddo
2720  ir=20
2721 1000 continue
2722  return
2723  end subroutine addr0
2724 
2725  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2726  subroutine vcopy(a,c,n)
2727  implicit none
2728 
2729  integer(kind=kint) :: n
2730  real(kind=kreal) :: a(n),c(n)
2731  ! do 100 i=1,n
2732  ! c(i)=a(i)
2733  ! 100 continue
2734  c=a
2735  return
2736  end subroutine vcopy
2737 
2738  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2739 
2740  subroutine verif0(ndeg,neqns_a0,nttbr_a0,irow_a0,jcol_a0,val_a0,neqns_l,nttbr_l,irow_l,jcol_l,val_l,rhs,x)
2741 
2742  implicit none
2743 
2744  integer(kind=kint), intent(in) :: ndeg
2745 
2746  integer(kind=kint), intent(in) :: irow_a0(:),jcol_a0(:)
2747  integer(kind=kint), intent(in) :: neqns_a0,nttbr_a0
2748  real(kind=kreal), intent(in) :: val_a0(:,:)
2749  integer(kind=kint), intent(in) :: irow_l(:),jcol_l(:)
2750  integer(kind=kint), intent(in) :: neqns_l,nttbr_l
2751  real(kind=kreal), intent(in) :: val_l(:,:)
2752 
2753  real(kind=kreal), intent(in) :: x(:,:)
2754  real(kind=kreal), intent(out) :: rhs(:,:)
2755 
2756  integer(kind=kint) :: i,j,k,l,m
2757  real(kind=kreal) :: rel,err
2758  !
2759  !----------------------------------------------------------------------
2760  !
2761  ! verify the solution(symmetric matrix)
2762  !
2763  ! ndeg=1 only
2764  !
2765  ! include lagrange elements as L region.
2766  !
2767  ! A0 | +
2768  ! A0 | |
2769  ! A0 | | neqns_a0
2770  ! A0 | |
2771  ! A0| +
2772  ! ----------+--- -
2773  ! |0 +
2774  ! | 0 | neqns_l
2775  ! laglange | 0 +
2776  !
2777  !
2778  !----------------------------------------------------------------------
2779  !
2780  rel=0.0d0
2781  do i=1,neqns_a0+neqns_l
2782  do l=1,ndeg
2783  rel=rel+dabs(rhs(l,i))
2784  enddo
2785  enddo
2786 
2787  ! A0 region
2788  do k=1,nttbr_a0
2789  i=irow_a0(k)
2790  j=jcol_a0(k)
2791  do l=1,ndeg
2792  do m=1,ndeg
2793  rhs(l,i)=rhs(l,i)-val_a0(1,k)*x(m,j)
2794  if(i.ne.j) rhs(l,j)=rhs(l,j)-val_a0(1,k)*x(m,i)
2795  end do
2796  end do
2797  end do
2798 
2799  ! lagrange region
2800  do k=1,nttbr_l
2801  i=irow_l(k)+neqns_a0
2802  j=jcol_l(k)
2803  do l=1,ndeg
2804  do m=1,ndeg
2805  rhs(l,i)=rhs(l,i)-val_l(1,k)*x(m,j)
2806  if(i.ne.j) rhs(l,j)=rhs(l,j)-val_l(1,k)*x(m,i)
2807  end do
2808  end do
2809  end do
2810 
2811  err=0.0d0
2812  do i=1,neqns_a0 + neqns_l
2813  do l=1,ndeg
2814  err=err+dabs(rhs(l,i))
2815  enddo
2816  enddo
2817 
2818  write(imsg,6000) err,rel,err/rel
2819 6000 format(' ***verification***(symmetric)'/&
2820  & 'norm(Ax-b) = ',1pd20.10/&
2821  & 'norm(b) = ',1pd20.10/&
2822  & 'norm(Ax-b)/norm(b) = ',1pd20.10)
2823 6010 format(1p4d15.7)
2824  return
2825  end subroutine verif0
2826 
2827 end module hecmw_solver_direct_serial_lag
verif0
subroutine verif0(Neqns, Ndeg, Nttbr, Irow, Jcol, Val, Rhs, X)
Definition: _unused_code.f90:301
idntty
subroutine idntty(Neqns, Invp, Iperm)
Definition: _unused_code.f90:98
hecmw_solver_direct_serial_lag
Definition: hecmw_solver_direct_serial_lag.f90:6
hecmw_solver_direct_serial_lag::hecmw_solve_direct_serial_lag
subroutine, public hecmw_solve_direct_serial_lag(nrows, ilag_sta, nttbr, pointers, indices, values, b)
Definition: hecmw_solver_direct_serial_lag.f90:57
hecmw_solver_direct_serial_lag::qqsort
subroutine qqsort(iw, ik)
Definition: hecmw_solver_direct_serial_lag.f90:2196
m_child_matrix_lag
Definition: m_child_matrix_lag.f90:5
m_crs_matrix_lag
Definition: m_crs_matrix_lag.f90:5
m_crs_matrix_lag::symbolicirjctocrs
subroutine, public symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
Definition: m_crs_matrix_lag.f90:29
m_irjc_matrix_lag
Definition: m_irjc_matrix_lag.f90:5
my_hecmw_util_lag
Definition: my_hecmw_util_f_lag.f90:6
my_hecmw_util_lag::kreal
integer(kind=4), parameter kreal
Definition: my_hecmw_util_f_lag.f90:10
m_child_matrix_lag::child_matrix
Definition: m_child_matrix_lag.f90:10
m_irjc_matrix_lag::irjc_mn_matrix
Definition: m_irjc_matrix_lag.f90:26
m_irjc_matrix_lag::irjc_square_matrix
Definition: m_irjc_matrix_lag.f90:13