FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
7
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 ! begining 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
54contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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 !! temporaly
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 givem 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 ! bcause of lagrange
305 diag_lag=0.0d0 ! bcause 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 informations
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 skelton 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 informations 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 ! Backword 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
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 specilized 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
subroutine verif0(neqns, ndeg, nttbr, irow, jcol, val, rhs, x)
subroutine idntty(neqns, invp, iperm)
subroutine, public hecmw_solve_direct_serial_lag(nrows, ilag_sta, nttbr, pointers, indices, values, b)
subroutine, public symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
integer(kind=4), parameter kreal