FrontISTR  5.7.0
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 100 ic=2,nstop-1
909  call sum(ic,xlnzr,colno,zln(1,:),diag(1,:),nch,parent,neqns)
910  100 continue
911 
912  !
913  ! phase II
914  ! LDU decompose of C (nstop..neqnsA+neqnsd)
915  !
916  do 200 ic=nstop,neqns
917  call sum1(ic,xlnzr,colno,zln(1,:),diag(1,:),parent,neqns)
918  200 continue
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 100 l=1,neqns
1056  zpiv(l)=1
1057  100 continue
1058 
1059  do 200 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  200 continue
1071 
1072  do 310 i=neqns,1,-1
1073  if(zpiv(i).eq.0) then
1074  neqnsz=i
1075  goto 320
1076  endif
1077  310 continue
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 10 i=1,2*nttbr
1099  jcpt(i)=0
1100  jcolno(i)=0
1101  10 continue
1102  do 20 i=1,neqns
1103  jcpt(i)=i+neqns
1104  jcolno(i+neqns)=i
1105  20 continue
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 100 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  100 continue
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 100 i=1,neqns
1219  invp(i)=i
1220  iperm(i)=i
1221  100 continue
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 10 i=1,xadj(neqns+1)-1
1242  adjncy(i)=adj0(i)
1243  10 continue
1244  do 100 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  100 continue
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 100 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  100 continue
1336  do 200 i=1,neqns
1337  if(parent(i).eq.0) parent(i)=neqns+1
1338  200 continue
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 10 i=1,neqns+1
1361  btree(1,i)=0
1362  btree(2,i)=0
1363  10 continue
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 200 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  200 continue
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 210 i=1,neqns
1452  adjt(i)=0
1453  210 continue
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 240 k=xadj(iperm(loc)),xadj(iperm(loc)+1)-1
1461  adjt(invp(adjncy(k)))=1
1462  240 continue
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 260 ll=l+1,nanc
1470  if(adjt(anc(ll)).eq.0) then
1471  l=l+1
1472  goto 200
1473  endif
1474  260 continue
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 300 i=1,neqns
1483  adjt(i)=0
1484  300 continue
1485  do 310 ll=l,nanc
1486  adjt(anc(ll))=1
1487  310 continue
1488  k=0
1489  do 320 ll=1,neqns
1490  if(adjt(ll).eq.0) then
1491  k=k+1
1492  invp(iperm(ll))=k
1493  endif
1494  320 continue
1495  ! (2) followed by nodes in Ancestor(iy)-Adj(T(iy))
1496  330 continue
1497  do 340 i=1,neqns
1498  adjt(i)=0
1499  340 continue
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 370 kk=xadj(iperm(loc)),xadj(iperm(loc)+1)-1
1507  adjt(invp(adjncy(kk)))=1
1508  370 continue
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 390 ll=l,nanc
1516  if(adjt(anc(ll)).eq.0) then
1517  k=k+1
1518  invp(iperm(anc(ll)))=k
1519  endif
1520  390 continue
1521  ! (3) and finally number the node in Adj(t(iy))
1522  do 400 ll=l,nanc
1523  if(adjt(anc(ll)).ne.0) then
1524  k=k+1
1525  invp(iperm(anc(ll)))=k
1526  endif
1527  400 continue
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 610 i=1,neqns
1545  iperm(invp(i))=i
1546  610 continue
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 200 i=1,neqns
1594  if(invp(iperm(i)).ne.i) goto 210
1595  if(iperm(invp(i)).ne.i) goto 210
1596  200 continue
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 5 i=1,neqns
1621  mch(i)=0
1622  pordr(i)=0
1623  5 continue
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 100 i=1,neqns
1644  ipinv=pordr(invp(i))
1645  invp(i)=ipinv
1646  iperm(ipinv)=i
1647  iw(pordr(i))=i
1648  100 continue
1649  do 110 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  110 continue
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 105 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  105 continue
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 130 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  130 continue
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 200 k=1,neqns
1837  write(idbg,6100) k
1838  write(idbg,6000) (colno(i),i=xlnzr(k),xlnzr(k+1)-1)
1839  200 continue
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  400 if(marker(node).ne.0) goto 500
1885  rchsze=rchsze+1
1886  rchset(rchsze)=node
1887  marker(node)=1
1888  500 continue
1889  600 continue
1890  return
1891  end subroutine qmdrch
1892 
1893  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1894 
1895  subroutine qmdupd(xadj,adjncy,nlist,list,deg,qsize,qlink,marker,rchset,nbrhd)
1896 
1897  implicit none
1898 
1899  integer(kind=kint), intent(in) :: adjncy(:),list(:),xadj(:)
1900  integer(kind=kint), intent(out) :: marker(:),nbrhd(:),rchset(:),deg(:),qsize(:),qlink(:)
1901  integer(kind=kint), intent(in) :: nlist
1902 
1903  integer(kind=kint) :: i,j,k,l, deg0,deg1,il,inhd,inode,irch,jstrt,jstop,mark,nabor,nhdsze,node,rchsze
1904 
1905  if(nlist.le.0) return
1906  deg0=0
1907  nhdsze=0
1908  do 200 il=1,nlist
1909  node=list(il)
1910  deg0=deg0+qsize(node)
1911  jstrt=xadj(node)
1912  jstop=xadj(node+1)-1
1913 
1914  do 100 j=jstrt,jstop
1915  nabor=adjncy(j)
1916  if(marker(nabor).ne.0.or.deg(nabor).ge.0) goto 100
1917  marker(nabor)=-1
1918  nhdsze=nhdsze+1
1919  nbrhd(nhdsze)=nabor
1920  100 continue
1921  200 continue
1922 
1923  if(nhdsze.gt.0) call qmdmrg(xadj,adjncy,deg,qsize,qlink,marker,deg0,nhdsze,nbrhd,rchset,nbrhd(nhdsze+1:))
1924  do 600 il=1,nlist
1925  node=list(il)
1926  mark=marker(node)
1927  if(mark.gt.1.or.mark.lt.0) goto 600
1928  call qmdrch(node,xadj,adjncy,deg,marker,rchsze,rchset,nhdsze,nbrhd)
1929  deg1=deg0
1930  if(rchsze.le.0) goto 400
1931  do 300 irch=1,rchsze
1932  inode=rchset(irch)
1933  deg1=deg1+qsize(inode)
1934  marker(inode)=0
1935  300 continue
1936  400 deg(node)=deg1-1
1937  if(nhdsze.le.0) goto 600
1938  do 500 inhd=1,nhdsze
1939  inode=nbrhd(inhd)
1940  marker(inode)=0
1941  500 continue
1942  600 continue
1943  return
1944  end subroutine qmdupd
1945 
1946  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1947 
1948  subroutine qmdot(root,xadj,adjncy,marker,rchsze,rchset,nbrhd)
1949 
1950  implicit none
1951 
1952  integer(kind=kint), intent(in) :: marker(:),rchset(:),nbrhd(:),xadj(:)
1953  integer(kind=kint), intent(out) :: adjncy(:)
1954  integer(kind=kint), intent(in) :: rchsze,root
1955 
1956  integer(kind=kint) :: i,j,k,l,irch,inhd,node,jstrt,jstop,link,nabor
1957 
1958  irch=0
1959  inhd=0
1960  node=root
1961  100 jstrt=xadj(node)
1962  jstop=xadj(node+1)-2
1963  if(jstop.lt.jstrt) goto 300
1964  do 200 j=jstrt,jstop
1965  irch=irch+1
1966  adjncy(j)=rchset(irch)
1967  if(irch.ge.rchsze) goto 400
1968  200 continue
1969  300 link=adjncy(jstop+1)
1970  node=-link
1971  if(link.lt.0) goto 100
1972  inhd=inhd+1
1973  node=nbrhd(inhd)
1974  adjncy(jstop+1)=-node
1975  goto 100
1976  400 adjncy(j+1)=0
1977  do 600 irch=1,rchsze
1978  node=rchset(irch)
1979  if(marker(node).lt.0) goto 600
1980  jstrt=xadj(node)
1981  jstop=xadj(node+1)-1
1982  do 500 j=jstrt,jstop
1983  nabor=adjncy(j)
1984  if(marker(nabor).ge.0) goto 500
1985  adjncy(j)=root
1986  goto 600
1987  500 continue
1988  600 continue
1989  return
1990  end subroutine qmdot
1991 
1992  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1993 
1994  subroutine qmdmrg(xadj,adjncy,deg,qsize,qlink,marker,deg0,nhdsze,nbrhd,rchset,ovrlp)
1995 
1996  implicit none
1997 
1998  integer(kind=kint), intent(in) :: adjncy(:),nbrhd(:),xadj(:)
1999  integer(kind=kint), intent(out) :: deg(:),marker(:),rchset(:),ovrlp(:),qsize(:),qlink(:)
2000  integer(kind=kint), intent(in) :: nhdsze
2001 
2002  integer(kind=kint) :: i,j,k,l, deg0,deg1,head,inhd,iov,irch,jstrt,jstop,link,lnode,mark,mrgsze,nabor,node,novrlp,rchsze,root
2003 
2004 
2005  if(nhdsze.le.0) return
2006  do 100 inhd=1,nhdsze
2007  root=nbrhd(inhd)
2008  marker(root)=0
2009  100 continue
2010  do 1400 inhd=1,nhdsze
2011  root=nbrhd(inhd)
2012  marker(root)=-1
2013  rchsze=0
2014  novrlp=0
2015  deg1=0
2016  200 jstrt=xadj(root)
2017  jstop=xadj(root+1)-1
2018  do 600 j=jstrt,jstop
2019  nabor=adjncy(j)
2020  root=-nabor
2021  if(nabor) 200,700,300
2022  300 mark=marker(nabor)
2023  if(mark)600,400,500
2024  400 rchsze=rchsze+1
2025  rchset(rchsze)=nabor
2026  deg1=deg1+qsize(nabor)
2027  marker(nabor)=1
2028  goto 600
2029  500 if(mark.gt.1) goto 600
2030  novrlp=novrlp+1
2031  ovrlp(novrlp)=nabor
2032  marker(nabor)=2
2033  600 continue
2034  700 head=0
2035  mrgsze=0
2036  do 1100 iov=1,novrlp
2037  node=ovrlp(iov)
2038  jstrt=xadj(node)
2039  jstop=xadj(node+1)-1
2040  do 800 j=jstrt,jstop
2041  nabor=adjncy(j)
2042  if(marker(nabor).ne.0) goto 800
2043  marker(node)=1
2044  goto 1100
2045  800 continue
2046  mrgsze=mrgsze+qsize(node)
2047  marker(node)=-1
2048  lnode=node
2049  900 link=qlink(lnode)
2050  if(link.le.0) goto 1000
2051  lnode=link
2052  goto 900
2053  1000 qlink(lnode)=head
2054  head=node
2055  1100 continue
2056  if(head.le.0) goto 1200
2057  qsize(head)=mrgsze
2058  deg(head)=deg0+deg1-1
2059  marker(head)=2
2060  1200 root=nbrhd(inhd)
2061  marker(root)=0
2062  if(rchsze.le.0) goto 1400
2063  do 1300 irch=1,rchsze
2064  node=rchset(irch)
2065  marker(node)=0
2066  1300 continue
2067  1400 continue
2068  return
2069  end subroutine qmdmrg
2070 
2071  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2072 
2073  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)
2074  ! find fill-in position in C which placed under A and set it in xlnzr_c, colno_c
2075 
2076  use m_crs_matrix_lag
2077  implicit none
2078  !input
2079  integer(kind=kint), intent(in) :: xlnzr_a(:)
2080  integer(kind=kint), intent(in) :: colno_a(:)
2081  integer(kind=kint), intent(in) :: iperm_a(:)
2082  integer(kind=kint), intent(in) :: invp_a(:)
2083  integer(kind=kint), intent(in) :: ndeg
2084  integer(kind=kint), intent(in) :: nttbr_c
2085  integer(kind=kint), intent(in) :: irow_c(:)
2086  integer(kind=kint), intent(inout) :: jcol_c(:)
2087  integer(kind=kint), intent(in) :: ncol
2088  integer(kind=kint), intent(in) :: nrow
2089 
2090  !output
2091  integer(kind=kint), pointer :: xlnzr_c(:)
2092  integer(kind=kint), pointer :: colno_c(:)
2093  integer(kind=kint), intent(out) :: lncol_c
2094 
2095  ! internal
2096  integer(kind=kint) :: i,j,k,l,m,n
2097  integer(kind=kint) :: ks, ke, ipass, ierr
2098  logical, allocatable :: cnz(:)
2099  type(crs_matrix) :: crs_c
2100 
2101  !permtate column in C for crs_c
2102  do i=1,nttbr_c
2103  jcol_c(i)=invp_a(jcol_c(i))
2104  end do
2105 
2106  ! make Compact Column Storoge using symbolic information.
2107  call symbolicirjctocrs(ndeg, nttbr_c, irow_c, jcol_c, ncol, nrow, crs_c)
2108 
2109  ! symbolic LDU factorization for C matrix
2110  allocate(cnz(ncol), stat=ierr)
2111  if(ierr .ne. 0) then
2112  call errtrp('stop due to allocation error.')
2113  end if
2114  do ipass = 1,2
2115  lncol_c = 0
2116  do k=1,nrow
2117  ! set cnz as non-zero pattern of C
2118  cnz = .false.
2119  ks = crs_c%ia(k)
2120  ke = crs_c%ia(k+1)-1
2121  if (ke .lt. ks) then
2122  if (ipass .eq. 2) then
2123  xlnzr_c(k+1)=lncol_c+1
2124  end if
2125  cycle ! in case of zero vector, no need to check dot product. not cycle?
2126  end if
2127 
2128  do i=ks,ke
2129  cnz(crs_c%ja(i)) = .true.
2130  end do
2131 
2132  ! check for non-zero dot product and update cnz for each point of cnz
2133  do i=2,ncol
2134  ks = xlnzr_a(i)
2135  ke = xlnzr_a(i+1)-1
2136  if (ke .lt. ks) then ! in case of column of A is zero vector.
2137  cycle
2138  end if
2139  do j=ks,ke
2140  if (cnz(colno_a(j))) then
2141  cnz(i) = .true.
2142  exit
2143  end if
2144  end do
2145  end do
2146 
2147  do i=1,ncol
2148  if (cnz(i)) then
2149  lncol_c = lncol_c + 1
2150  if (ipass .eq. 2) then
2151  colno_c(lncol_c) = i
2152  end if
2153  end if
2154  end do
2155  if (ipass .eq. 2) then
2156  xlnzr_c(k+1)=lncol_c + 1
2157  end if
2158  end do
2159 
2160  if (ipass .eq. 1) then
2161  allocate(xlnzr_c(nrow+1),colno_c(lncol_c), stat=ierr)
2162  if(ierr .ne. 0) then
2163  call errtrp('stop due to allocation error.')
2164  end if
2165  xlnzr_c(1)=1
2166  end if
2167  end do
2168 
2169  ! restore order of C column.
2170  do i=1,nttbr_c
2171  jcol_c(i)=iperm_a(jcol_c(i))
2172  end do
2173 
2174  return
2175 
2176  end subroutine ldudecomposec
2177 
2178  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2179 
2180  subroutine qqsort(iw,ik)
2181 
2182  implicit none
2183 
2184  integer(kind=kint), intent(out) :: iw(:)
2185  integer(kind=kint), intent(in) :: ik
2186 
2187  integer(kind=kint) :: l,m,itemp
2188 
2189  !----------------------------------------------------------------------
2190  ! sort in increasing order up to i
2191  !
2192  ! iw array
2193  ! ik number of input/output
2194  ! i deal with numbers less than this numberi
2195  !
2196  !----------------------------------------------------------------------
2197 
2198  if(ik.le.1) return
2199  do 100 l=1,ik-1
2200  do 110 m=l+1,ik
2201  if(iw(l).lt.iw(m)) goto 110
2202  itemp=iw(l)
2203  iw(l)=iw(m)
2204  iw(m)=itemp
2205  110 continue
2206  100 continue
2207  200 continue
2208  return
2209  end subroutine qqsort
2210 
2211 
2212  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2213 
2214  subroutine staij1(isw,i,j,aij,dsi,ir)
2215 
2216  implicit none
2217 
2218  !----------------------------------------------------------------------
2219  !
2220  ! this routine sets an non-zero entry of the matrix.
2221  ! (symmetric version)
2222  !
2223  ! (i)
2224  ! isw =0 set the value
2225  ! =1 add the value
2226  ! i row entry
2227  ! j column entry
2228  ! aij value
2229  !
2230  ! (o)
2231  ! iv communication array
2232  !
2233  ! #coded by t.arakawa
2234  !
2235  !----------------------------------------------------------------------
2236  !
2237  type(dsinfo) :: dsi
2238  real(kind=kreal), intent(out) :: aij(:) ! ndeg*ndeg
2239  integer(kind=kint), intent(in) :: isw, i, j
2240  integer(kind=kint), intent(out) :: ir
2241 
2242  integer(kind=kint) :: ndeg, neqns, nstop, ndeg2, ndeg2l, ierr
2243  ndeg=dsi%ndeg
2244  neqns=dsi%neqns
2245  nstop=dsi%nstop
2246  ndeg2=ndeg*ndeg
2247  ndeg2l=ndeg*(ndeg+1)/2
2248 
2249  ir=0
2250  ierr=0
2251 
2252 
2253  ! array allocation
2254  if(dsi%stage.ne.20) then
2255  if(dsi%stage.eq.30) write(ilog,*) 'Warning a matrix was build up but never solved.'
2256  !
2257  ! for diagonal
2258  !
2259  allocate(dsi%diag(ndeg2l,neqns), stat=ierr)
2260  if(ierr .ne. 0) then
2261  call errtrp('stop due to allocation error.')
2262  end if
2263  dsi%diag=0
2264  !
2265  ! for lower triangle
2266  !
2267  allocate(dsi%zln(ndeg2,dsi%lncol), stat=ierr)
2268  if(ierr .ne. 0) then
2269  call errtrp('stop due to allocation error.')
2270  end if
2271  dsi%zln=0
2272  !
2273  ! for dense window !TODO delete this and corresponding line in addr3()
2274  !
2275  ! allocate(dsi%dsln(ndeg2,dsi%lndsln))! because of there is no dense window
2276  ! dsi%dsln=0
2277 
2278  dsi%stage=20
2279  endif
2280 
2281  ! set value
2282  if(ndeg.le.2) then
2283  call addr0(isw,i,j,aij,dsi%invp,dsi%xlnzr,dsi%colno,dsi%diag,dsi%zln,dsi%dsln,nstop,dsi%ndeg,ir)
2284  elseif(ndeg.eq.3) then
2285  write(idbg,*) 'ndeg=1 only'
2286  stop
2287  else
2288  write(idbg,*) 'ndeg=1 only'
2289  stop
2290  endif
2291  1000 continue
2292  return
2293  end subroutine staij1
2294 
2295  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2296  !
2297  ! After here, routines specialized for ndeg = 1
2298  !
2299  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2300 
2301  ! LDU decompose of A (1..nstop-1) region
2302  subroutine sum(ic,xlnzr,colno,zln,diag,nch,par,neqns)
2303 
2304  implicit none
2305 
2306  integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
2307  integer(kind=kint), intent(in) :: ic, neqns
2308  real(kind=kreal), intent(inout) :: zln(:),diag(:)
2309  integer(kind=kint), intent(out) :: nch(:)
2310 
2311  real(kind=kreal) :: s, t, zz, piv
2312  integer(kind=kint) :: ks, ke, kk, k, jc, jj, j, ierr
2313  integer(kind=kint) :: isem
2314  real(kind=kreal),allocatable :: temp(:)
2315  integer(kind=kint),allocatable :: indx(:)
2316  allocate(temp(neqns),indx(neqns), stat=ierr)
2317  if(ierr .ne. 0) then
2318  call errtrp('stop due to allocation error.')
2319  end if
2320 
2321  2 continue
2322  ks=xlnzr(ic)
2323  ke=xlnzr(ic+1)
2324  t=0.0d0
2325  ! do 100 i=1,ic
2326  ! temp(i)=0.0d0
2327  ! 100 continue
2328  do 200 k=ks,ke-1
2329  jc=colno(k)
2330  indx(jc)=ic
2331  s=0.0d0
2332  do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
2333  j=colno(jj)
2334  if(indx(j).eq.ic) then
2335  s=s+temp(j)*zln(jj)
2336  endif
2337  310 continue
2338  ! j1=xlnzr(jc)
2339  ! jj=xlnzr(jc+1)-j1
2340  ! ss=ddoti(jj,zln(j1),colno(j1),temp)
2341  ! zz=zln(k)-ddoti(jj,zln(j1),colno(j1),temp)
2342  zz=zln(k)-s
2343  zln(k)=zz*diag(jc)
2344  temp(jc)=zz
2345  t=t+zz*zln(k)
2346  200 continue
2347  piv=diag(ic)-t
2348  if(dabs(piv).gt.rmin) then
2349  diag(ic)=1.0d0/piv
2350  endif
2351  1 continue
2352  ! if(isem.eq.1) then !DBG
2353  isem=0
2354  nch(ic)=-1
2355  kk=par(ic)
2356  nch(kk)=nch(kk)-1
2357  isem=1
2358  ! else !DBG
2359  ! goto 1 !DBG
2360  ! endi !DBG
2361  return
2362  end subroutine sum
2363 
2364  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2365 
2366  ! LDU decompose of C (nstop..neqnsA+neqnsd) region
2367  subroutine sum1(ic,xlnzr,colno,zln,diag,par,neqns)
2368 
2369  implicit none
2370 
2371  integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
2372  integer(kind=kint), intent(in) :: ic, neqns
2373  real(kind=kreal), intent(inout) :: zln(:),diag(:)
2374 
2375  real(kind=kreal) :: s, t, zz
2376  integer(kind=kint) :: ks, ke, k, jc, j, jj, ierr
2377  real(kind=kreal),allocatable :: temp(:)
2378  integer(kind=kint),allocatable :: indx(:)
2379  integer(kind=kint) :: i
2380 
2381  ierr=0
2382 
2383  allocate(temp(neqns),indx(neqns), stat=ierr)
2384  if(ierr .ne. 0) then
2385  call errtrp('stop due to allocation error.')
2386  end if
2387 
2388  do i=1,neqns
2389  temp(i)=0
2390  end do
2391 
2392  ks=xlnzr(ic)
2393  ke=xlnzr(ic+1)
2394  t=0.0d0
2395  ! do 100 i=1,ic
2396  ! temp(i)=0.0d0
2397  ! 100 continue
2398  do 200 k=ks,ke-1
2399  jc=colno(k)
2400  indx(jc)=ic
2401  s=0.0d0
2402  do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
2403  j=colno(jj)
2404  if(indx(j).eq.ic) then
2405  s=s+temp(j)*zln(jj)
2406  endif
2407  310 continue
2408  zz=zln(k)-s
2409  ! j1=xlnzr(jc)
2410  ! jj=xlnzr(jc+1)-j1
2411  ! zz=zln(k)-ddoti(jj,zln(j1),colno(j1),temp)
2412  zln(k)=zz
2413  temp(jc)=zz
2414  ! t=t+zz*zz*diag(jc)
2415  200 continue
2416  ! diag(ic)=diag(ic)-t
2417  return
2418  end subroutine sum1
2419 
2420  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2421 
2422  ! LDU decompose and Update D region.
2423  subroutine sum2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln)
2424 
2425  implicit none
2426 
2427  integer(kind=kint), intent(in) :: neqns, nstop
2428  integer(kind=kint), intent(in) :: xlnzr(:),colno(:)
2429  real(kind=kreal), intent(inout) :: zln(:),diag(:)
2430  integer(kind=kint), pointer :: spdslnidx(:)
2431  real(kind=kreal), pointer :: spdslnval(:,:)
2432  integer(kind=kint), intent(out) :: nspdsln
2433 
2434  real(kind=kreal) :: s, t
2435  integer(kind=kint) :: ks, ke, kk, k, jc, jj, j, j1,j2
2436  integer(kind=kint) :: ic, i, loc, ierr
2437  integer(kind=kint) :: ispdsln
2438  logical :: ftflag
2439  real(kind=kreal),allocatable :: temp(:)
2440  integer(kind=kint),allocatable :: indx(:)
2441  ierr=0
2442  allocate(temp(neqns),indx(neqns), stat=ierr)
2443  if(ierr .ne. 0) then
2444  call errtrp('stop due to allocation error.')
2445  end if
2446  temp=0
2447 
2448  nspdsln=0
2449  do ic=nstop,neqns
2450  ks=xlnzr(ic)
2451  ke=xlnzr(ic+1)-1
2452  do k=ks,ke
2453  jj=colno(k)
2454  indx(jj)=ic
2455  end do
2456  do jc=nstop,ic-1
2457  j1=xlnzr(jc)
2458  j2=xlnzr(jc+1)
2459  do jj=xlnzr(jc),xlnzr(jc+1)-1
2460  j=colno(jj)
2461  if(indx(j).eq.ic) then
2462  nspdsln=nspdsln+1
2463  exit
2464  endif
2465  end do
2466  end do
2467  end do
2468  allocate(spdslnidx(nspdsln),spdslnval(1,nspdsln), stat=ierr)
2469  if(ierr .ne. 0) then
2470  call errtrp('stop due to allocation error.')
2471  end if
2472 
2473  loc=0
2474  ispdsln=0
2475  spdslnval=0
2476  ftflag = .true.
2477  do 100 ic=nstop,neqns
2478  do 105 i=1,nstop
2479  temp(i)=0.0d0
2480  105 continue
2481  ks=xlnzr(ic)
2482  ke=xlnzr(ic+1)-1
2483  do 110 k=ks,ke
2484  jj=colno(k)
2485  temp(jj)=zln(k)
2486  zln(k)=temp(jj)*diag(jj)
2487  indx(jj)=ic
2488  diag(ic)=diag(ic)-temp(jj)*zln(k)
2489  110 continue
2490  do 120 jc=nstop,ic-1
2491  loc=loc+1
2492  do 220 jj=xlnzr(jc),xlnzr(jc+1)-1
2493  j=colno(jj)
2494  if(indx(j).eq.ic) then
2495  if (ftflag) then
2496  ispdsln=ispdsln+1
2497  ftflag=.false.
2498  end if
2499  spdslnidx(ispdsln)=loc
2500  spdslnval(1,ispdsln)=spdslnval(1,ispdsln)-temp(j)*zln(jj)
2501  endif
2502  220 continue
2503  ftflag = .true.
2504  ! j1=xlnzr(jc)
2505  ! jj=xlnzr(jc+1)-j1
2506  ! dsln(loc)=dsln(loc)-ddoti(jj,zln(j1),colno(j1),temp)
2507  120 continue
2508  100 continue
2509  return
2510  end subroutine sum2_child
2511 
2512  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2513 
2514  subroutine sum3(n,dsln,diag)
2515 
2516  implicit none
2517 
2518  real(kind=kreal), intent(inout) :: dsln(:),diag(:)
2519  integer(kind=kint), intent(in) :: n
2520 
2521  integer(kind=kint) :: i, j, loc, ierr
2522  real(kind=kreal),allocatable :: temp(:)
2523  integer(kind=kint),allocatable :: indx(:)
2524  allocate(temp(n),indx(n), stat=ierr)
2525  if(ierr .ne. 0) then
2526  call errtrp('stop due to allocation error.')
2527  end if
2528 
2529  if(n.le.0) goto 1000
2530  indx(1)=0
2531  loc=1
2532  diag(1)=1.0d0/diag(1)
2533  do 100 i=2,n
2534  indx(i)=loc
2535  do 110 j=1,i-1
2536  dsln(loc)=dsln(loc)-dot_product(dsln(indx(i):indx(i)+j-2),dsln(indx(j):indx(j)+j-2))
2537  loc=loc+1
2538  110 continue
2539  temp(1:i-1)=dsln(indx(i):indx(i)+i-2)*diag(1:i-1)
2540  diag(i)=diag(i)-dot_product(temp(1:i-1),dsln(indx(i):indx(i)+i-2))
2541  dsln(indx(i):indx(i)+i-2)=temp(1:i-1)
2542  diag(i)=1.0d0/diag(i)
2543  100 continue
2544  1000 continue
2545 
2546  return
2547  end subroutine sum3
2548 
2549  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2550 
2551  real(kind=kreal) function spdot2(b,zln,colno,ks,ke)
2552 
2553  implicit none
2554 
2555  integer(kind=kint), intent(in) :: colno(:)
2556  integer(kind=kint), intent(in) :: ks,ke
2557  real(kind=kreal), intent(in) :: zln(:),b(:)
2558 
2559  integer(kind=kint) :: j,jj
2560  real(kind=kreal) :: s
2561 
2562  !----------------------------------------------------------------------
2563  !
2564  ! spdot1 performs inner product of sparse vectors
2565  !
2566  !
2567  ! #coded by t.arakawa
2568  !
2569  !----------------------------------------------------------------------
2570  !
2571  s=0.0d0
2572  do 100 jj=ks,ke
2573  j=colno(jj)
2574  s=s+zln(jj)*b(j)
2575  100 continue
2576  spdot2=s
2577  end function spdot2
2578 
2579  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2580 
2581  real(kind=kreal) function ddot(a,b,n)
2582 
2583  implicit none
2584 
2585  real(kind=kreal), intent(in) :: a(n),b(n)
2586  integer(kind=kint), intent(in) :: n
2587 
2588  real(kind=kreal) :: s
2589  integer(kind=kint) :: i
2590 
2591  s=0.0d0
2592  do 100 i=1,n
2593  s=s+a(i)*b(i)
2594  100 continue
2595  ddot=s
2596  return
2597  end function ddot
2598 
2599  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2600 
2601  subroutine addr0(isw,i,j,aij,invp,xlnzr,colno,diag,zln,dsln,nstop,ndeg,ir)
2602 
2603  implicit none
2604 
2605  integer(kind=kint), intent(in) :: isw ! 0: renew diag, dsln, zln other: add to diag, dsln, zln
2606  integer(kind=kint), intent(in) :: i,j,nstop, ndeg, invp(:),xlnzr(:),colno(:)
2607  real(kind=kreal), intent(inout) :: zln(:,:),diag(:,:),dsln(:,:),aij(:)
2608  integer(kind=kint), intent(out) :: ir
2609 
2610  integer(kind=kint) :: ndeg2, ii, jj, itrans, k, i0, j0, l, ks, ke
2611  integer(kind=kint), parameter :: idbg=0
2612  ndeg2=ndeg*ndeg
2613 
2614  ir=0
2615  ii=invp(i)
2616  jj=invp(j)
2617  if(idbg.ne.0) write(idbg,*) 'addr0',ii,jj,aij
2618  if(ii.eq.jj) then
2619  if(ndeg2.eq.1) then
2620  if(isw.eq.0) then
2621  diag(1,ii)=aij(1)
2622  else
2623  diag(1,ii)=diag(1,ii)+aij(1)
2624  endif
2625  elseif(ndeg2.eq.4) then
2626  if(isw.eq.0) then
2627  diag(1,ii)=aij(1)
2628  diag(2,ii)=aij(2)
2629  diag(3,ii)=aij(4)
2630  else
2631  diag(1,ii)=diag(1,ii)+aij(1)
2632  diag(2,ii)=diag(2,ii)+aij(2)
2633  diag(3,ii)=diag(3,ii)+aij(4)
2634  endif
2635  endif
2636  goto 1000
2637  endif
2638  itrans=0
2639  if(jj.gt.ii) then
2640  k=jj
2641  jj=ii
2642  ii=k
2643  itrans=1
2644  endif
2645  if(jj.ge.nstop) then
2646  i0=ii-nstop
2647  j0=jj-nstop+1
2648  k=i0*(i0-1)/2+j0
2649  if(ndeg2.eq.1) then
2650  dsln(1,k)=aij(1)
2651  goto 1000
2652  elseif(ndeg2.eq.4) then
2653  if(itrans.eq.0) then
2654  do 3 l=1,ndeg2
2655  dsln(l,k)=aij(l)
2656  3 continue
2657  goto 1000
2658  else
2659  dsln(1,k)=aij(1)
2660  dsln(2,k)=aij(3)
2661  dsln(3,k)=aij(2)
2662  dsln(4,k)=aij(4)
2663  goto 1000
2664  endif
2665  endif
2666  endif
2667  ks=xlnzr(ii)
2668  ke=xlnzr(ii+1)-1
2669  do 100 k=ks,ke
2670  if(colno(k).eq.jj) then
2671  if(isw.eq.0) then
2672  if(ndeg2.eq.1) then
2673  zln(1,k)=aij(1)
2674  elseif(ndeg2.eq.4) then
2675  if(itrans.eq.0) then
2676  do 4 l=1,ndeg2
2677  zln(l,k)=aij(l)
2678  4 continue
2679  else
2680  zln(1,k)=aij(1)
2681  zln(2,k)=aij(3)
2682  zln(3,k)=aij(2)
2683  zln(4,k)=aij(4)
2684  endif
2685  endif
2686  else
2687  if(ndeg2.eq.1) then
2688  zln(1,k)=zln(1,k)+aij(1)
2689  elseif(ndeg2.eq.4) then
2690  if(itrans.eq.0) then
2691  do 5 l=1,ndeg2
2692  zln(l,k)=zln(l,k)+aij(l)
2693  5 continue
2694  else
2695  zln(1,k)=zln(1,k)+aij(1)
2696  zln(2,k)=zln(2,k)+aij(3)
2697  zln(3,k)=zln(3,k)+aij(2)
2698  zln(4,k)=zln(4,k)+aij(4)
2699  endif
2700  endif
2701  endif
2702  goto 1000
2703  endif
2704  100 continue
2705  ir=20
2706  1000 continue
2707  return
2708  end subroutine addr0
2709 
2710  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2711  subroutine vcopy(a,c,n)
2712  implicit none
2713 
2714  integer(kind=kint) :: n
2715  real(kind=kreal) :: a(n),c(n)
2716  ! do 100 i=1,n
2717  ! c(i)=a(i)
2718  ! 100 continue
2719  c=a
2720  return
2721  end subroutine vcopy
2722 
2723  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2724 
2725  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)
2726 
2727  implicit none
2728 
2729  integer(kind=kint), intent(in) :: ndeg
2730 
2731  integer(kind=kint), intent(in) :: irow_a0(:),jcol_a0(:)
2732  integer(kind=kint), intent(in) :: neqns_a0,nttbr_a0
2733  real(kind=kreal), intent(in) :: val_a0(:,:)
2734  integer(kind=kint), intent(in) :: irow_l(:),jcol_l(:)
2735  integer(kind=kint), intent(in) :: neqns_l,nttbr_l
2736  real(kind=kreal), intent(in) :: val_l(:,:)
2737 
2738  real(kind=kreal), intent(in) :: x(:,:)
2739  real(kind=kreal), intent(out) :: rhs(:,:)
2740 
2741  integer(kind=kint) :: i,j,k,l,m
2742  real(kind=kreal) :: rel,err
2743  !
2744  !----------------------------------------------------------------------
2745  !
2746  ! verify the solution(symmetric matrix)
2747  !
2748  ! ndeg=1 only
2749  !
2750  ! include lagrange elements as L region.
2751  !
2752  ! A0 | +
2753  ! A0 | |
2754  ! A0 | | neqns_a0
2755  ! A0 | |
2756  ! A0| +
2757  ! ----------+--- -
2758  ! |0 +
2759  ! | 0 | neqns_l
2760  ! laglange | 0 +
2761  !
2762  !
2763  !----------------------------------------------------------------------
2764  !
2765  rel=0.0d0
2766  do 10 i=1,neqns_a0+neqns_l
2767  do 10 l=1,ndeg
2768  rel=rel+dabs(rhs(l,i))
2769  10 continue
2770 
2771  ! A0 region
2772  do k=1,nttbr_a0
2773  i=irow_a0(k)
2774  j=jcol_a0(k)
2775  do l=1,ndeg
2776  do m=1,ndeg
2777  rhs(l,i)=rhs(l,i)-val_a0(1,k)*x(m,j)
2778  if(i.ne.j) rhs(l,j)=rhs(l,j)-val_a0(1,k)*x(m,i)
2779  end do
2780  end do
2781  end do
2782 
2783  ! lagrange region
2784  do k=1,nttbr_l
2785  i=irow_l(k)+neqns_a0
2786  j=jcol_l(k)
2787  do l=1,ndeg
2788  do m=1,ndeg
2789  rhs(l,i)=rhs(l,i)-val_l(1,k)*x(m,j)
2790  if(i.ne.j) rhs(l,j)=rhs(l,j)-val_l(1,k)*x(m,i)
2791  end do
2792  end do
2793  end do
2794 
2795  err=0.0d0
2796  do 200 i=1,neqns_a0 + neqns_l
2797  do 200 l=1,ndeg
2798  err=err+dabs(rhs(l,i))
2799  200 continue
2800 
2801  write(imsg,6000) err,rel,err/rel
2802  6000 format(' ***verification***(symmetric)'/&
2803  & 'norm(Ax-b) = ',1pd20.10/&
2804  & 'norm(b) = ',1pd20.10/&
2805  & 'norm(Ax-b)/norm(b) = ',1pd20.10)
2806  6010 format(1p4d15.7)
2807  return
2808  end subroutine verif0
2809 
2810 end module hecmw_solver_direct_serial_lag
hecmw_solver_direct_serial_lag::qqsort
subroutine qqsort(iw, ik)
Definition: hecmw_solver_direct_serial_lag.f90:2181
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
m_child_matrix_lag
Definition: m_child_matrix_lag.f90:5
my_hecmw_util_lag
Definition: my_hecmw_util_f_lag.f90:6
verif0
subroutine verif0(Neqns, Ndeg, Nttbr, Irow, Jcol, Val, Rhs, X)
Definition: _unused_code.f90:301
m_crs_matrix_lag::symbolicirjctocrs
subroutine, public symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
Definition: m_crs_matrix_lag.f90:29
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
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
m_crs_matrix_lag
Definition: m_crs_matrix_lag.f90:5
m_child_matrix_lag::child_matrix
Definition: m_child_matrix_lag.f90:10
idntty
subroutine idntty(Neqns, Invp, Iperm)
Definition: _unused_code.f90:98
hecmw_solver_direct_serial_lag
Definition: hecmw_solver_direct_serial_lag.f90:6
m_irjc_matrix_lag
Definition: m_irjc_matrix_lag.f90:5