FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_solver_direct_parallel.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
8#ifndef HECMW_SERIAL
12 use m_elap
13#endif
14
15 use hecmw_util
16
17#ifndef HECMW_SERIAL
20#endif
21
22 ! access control !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
24 private ! default
25
26 public hecmw_solve_direct_parallel ! only entry point of Parallel Direct Solver is public
27
28#ifndef HECMW_SERIAL
29
30 ! internal type definition !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
31
32 type procinfo ! process informations on MPI
33 integer(kind=kint) :: myid
34 integer(kind=kint) :: imp ! mother process id
35 logical :: isparent ! true if this process is parent
36 logical :: ischild ! true if this process is child
37
38 integer(kind=kint) :: nchildren ! number of child process
39 integer(kind=kint), pointer :: ichildren(:) ! array of children process number
40
41 integer(kind=kint) :: ndiv ! count for matrix division.
42 end type procinfo
43
44 type dsinfo ! direct solver information
45 integer(kind=kint) :: ndeg ! dimension of small matrix
46 integer(kind=kint) :: neqns ! number of equations
47 integer(kind=kint) :: nstop ! begining point of C
48 integer(kind=kint) :: stage ! calculation stage
49 integer(kind=kint) :: lncol ! length of col
50 integer(kind=kint) :: lndsln ! length of dsln
51
52 integer(kind=kint), pointer :: zpiv(:) ! in zpivot()
53 integer(kind=kint), pointer :: iperm(:) ! permtation vector
54 integer(kind=kint), pointer :: invp(:) ! inverse permtation of iperm
55 integer(kind=kint), pointer :: parent(:) !
56 integer(kind=kint), pointer :: nch(:) !
57 integer(kind=kint), pointer :: xlnzr(:) ! ia index of whole sparse matrix. (neqns_t + 1)
58 integer(kind=kint), pointer :: colno(:) ! ja index of whole sparse matrix.
59
60 real(kind=kreal), pointer :: diag(:,:) ! diagonal element
61 real(kind=kreal), pointer :: zln(:,:) ! non diagonal sparse
62 real(kind=kreal), pointer :: dsln(:,:) ! non diagonal dens
63 end type dsinfo
64
65 ! internal global variables !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
66
67 type(procinfo) :: m_pds_procinfo ! initialized in initproc()
68 real(kind=kreal), parameter :: rmin = 1.00d-200 ! for inv3() pivot
69
70 integer :: imsg ! output file handler
71
72 integer, parameter :: ilog = 16 ! according to FSTR
73 logical, parameter :: ldbg = .false.
74 integer, parameter :: idbg = 52 ! according to FSTR
75 logical :: lelap = .false.
76
77#endif
78
79contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80
81 subroutine hecmw_solve_direct_parallel(hecMESH, hecMAT, ii)
82
83 implicit none
84
85#ifndef HECMW_SERIAL
86 include 'mpif.h'
87#endif
88
89 type (hecmwst_local_mesh), intent(inout) :: hecmesh
90 type (hecmwst_matrix ), intent(inout) :: hecmat
91 integer(kind=kint), intent(in) :: ii ! output file handler
92
93#ifndef HECMW_SERIAL
94
95 logical, save :: first_time = .true.
96
97 integer(kind=kint) :: ierr
98
99 ! start !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
100
101 call hecmw_mat_dump(hecmat, hecmesh)
102
103 imsg=ii ! set message file
104
105 ! set timelog
106 if (hecmat%Iarray(22) .ge. 1) then ! = timelog = kYES (kYES is defined in m_fstr
107 lelap = .true.
108 end if
109
110 ! set location in process tree
111 if (first_time) then
112 call initproc()
113 hecmat%Iarray(97) = 0 ! numeric factorization done flag
114 hecmat%Iarray(98) = 0 ! symbolic factorization done flag
115 first_time = .false.
116 end if
117
118 if ((hecmat%Iarray(97) .ne. 0) .or. (hecmat%Iarray(98) .ne. 0)) then
119 write(ilog,*) 'Error: Recalculation of LDU decompose is currently not surported'
121 end if
122
123 ! set elap time information
124 call initelap(lelap, idbg) !TODO it should be replaced with lelap
125
126 if (m_pds_procinfo%isparent) then
127 call elapout('hecmw_solve_direct_parallel: entering sp_direct_parent') !elap
128 call sp_direct_parent(hecmesh, hecmat)
129 else if (m_pds_procinfo%ischild) then
130 call elapout('hecmw_solve_direct_parallel: entering sp_direct_child') !elap
131 call sp_direct_child()
132 else
133 call elapout('hecmw_solve_direct_parallel: never come here') !elap
135 end if
136
137 call mpi_bcast(hecmat%x, hecmesh%n_dof*hecmat%NP, mpi_real8, m_pds_procinfo%imp, mpi_comm_world, ierr)
138
139 call hecmw_mat_dump_solution(hecmat)
140
141#endif
142
143 return
144 end subroutine hecmw_solve_direct_parallel
145
146
147#ifndef HECMW_SERIAL
148
149 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150
151 subroutine sp_direct_parent(hecMESH, hecMAT)
152
153 implicit none
154
155 include 'mpif.h'
156
157 !I/O !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
158 type (hecmwst_local_mesh), intent(inout) :: hecmesh
159 type (hecmwst_matrix ), intent(inout) :: hecmat
160
161
162 !internal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
163
164 ! given A0 x = b0
165 type (irjc_square_matrix) :: a0 ! given left side matrix assembled from sp_matrix
166 real(kind=kreal), allocatable :: b(:,:) ! (ndeg,neqns) right hand side value vector of equation.
167
168 ! for divided matrixes
169 type(matrix_partition_info) :: pmi
170 integer(kind=kint), pointer, save :: iperm_rev(:)
171 integer(kind=kint), pointer, save :: iofst_dm(:)
172
173 integer(kind=kint), save :: neqns_d ! number of eqns in D matrix
174 real(kind=kreal), pointer, save :: dsln(:,:) ! non-diagonal elements of dens D matrix
175 real(kind=kreal), pointer, save :: diag(:,:) ! diagonal elements of dens D matrix
176 integer(kind=kint), pointer, save :: part_all(:) ! index of corresponding dm of a0 row
177 integer(kind=kint), pointer, save :: iperm_all(:) ! index in partitioned matrix of a0 row
178 type(child_matrix), pointer, save :: dm(:) !divided matrices
179
180 real(kind=kreal), allocatable :: bd(:,:) ! for right hand side value
181
182 ! internal use
183 real(kind=kreal), allocatable :: oldb(:,:)
184 logical, save :: nusol_ready = .false.
185 integer(kind=kint), save :: ndeg, nndeg, ndegt
186 integer(kind=kint), save :: neqns_c, iofst_a2, iofst_c, ndm
187
188 ! misc !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
189 integer(kind=kint) :: ierr
190 integer(kind=kint) :: i,j,k,l,m,n
191
192 ! for MPI
193 integer(kind=kint) :: istatus(mpi_status_size)
194 integer(kind=kint) :: icp
195 real(kind=kreal), allocatable :: spdslnval(:,:), diagbuf(:,:), bdbuf(:,:)
196 integer(kind=kint), allocatable :: spdslnidx(:)
197 integer(kind=kint) :: nspdsln
198
199 ierr=0
200 call elapout('sp_direct_parent: entered') !elap
201
202
203 if (.not. nusol_ready) then
204
205 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
206 !
207 ! STEP01: get a0 from FEM data format hecMESH
208 !
209
210 call geta0(hecmesh, hecmat, a0)
211 ndeg=a0%ndeg
212 nndeg=ndeg*ndeg
213 ndegt = (ndeg+1)*ndeg/2 !triangle element in diag
214
215 call elapout('sp_direct_parent: make a0 done') !elap
216
217 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
218 !
219 ! STEP1X: Parallel LDU decompose of givem A0
220 !
221 ! STEP11: divide given large matrix A into a1 and a2 ... extend to 2^n division.
222 ! C is added to bottom of A1, A2.
223 ! D is kept as dsln (off-diagonal), diag (diagonal).
224
225 call elapout('sp_direct_parent: enter matrix partition') !elap
227 call elapout('sp_direct_parent: end matrix partition') !elap
228 neqns_d = pmi%neqns_d
229 dsln=>pmi%dsln
230 diag=>pmi%diag
231 part_all=>pmi%part_all
232 iperm_all=>pmi%iperm_all
233 dm=>pmi%dm
234 ndm=pmi%ndm
235
236 ! permtation vector for right hand side value
237 allocate(iofst_dm(0:ndm), stat=ierr)
238 if(ierr .ne. 0) then
239 call errtrp('stop due to allocation error.')
240 end if
241 iofst_dm(0)=0
242 iofst_dm(1)=neqns_d
243 do i=2,ndm
244 iofst_dm(i)=iofst_dm(i-1)+dm(i-1)%a%neqns
245 end do
246 do i=1,a0%neqns
247 iperm_all(i)=iperm_all(i)+iofst_dm(part_all(i))
248 end do
249
250 allocate(iperm_rev(a0%neqns), stat=ierr)
251 if(ierr .ne. 0) then
252 call errtrp('stop due to allocation error.')
253 end if
254 do i=1, a0%neqns
255 iperm_rev(iperm_all(i)) = i
256 end do
257
258 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
259 !
260 ! STEP12: Send divided left hand side matrixes to child processes.
261 !
262 ! nstop (separater of A and C) is also send
263 ! A and C will be LDU decomposed in child processes.
264 ! D region update data will be returned.
265
266 call elapout('sp_direct_parent: send divided matrix to children') !elap
267 do i=1,m_pds_procinfo%nchildren
268 icp = m_pds_procinfo%ichildren(i)
269 call mpi_send(dm(i)%a%ndeg, 1,mpi_integer,icp,1,mpi_comm_world,ierr)
270 call mpi_send(dm(i)%a%neqns, 1,mpi_integer,icp,1,mpi_comm_world,ierr)
271 call mpi_send(dm(i)%a%nttbr, 1,mpi_integer,icp,1,mpi_comm_world,ierr)
272
273 call mpi_send(dm(i)%a%irow, dm(i)%a%nttbr, mpi_integer,icp,1,mpi_comm_world,ierr)
274 call mpi_send(dm(i)%a%jcol, dm(i)%a%nttbr, mpi_integer,icp,1,mpi_comm_world,ierr)
275 call mpi_send(dm(i)%a%val, dm(i)%a%nttbr*dm(i)%a%ndeg*dm(i)%a%ndeg, mpi_real8,icp,1,mpi_comm_world,istatus,ierr)
276
277 call mpi_send(dm(i)%c%ndeg, 1,mpi_integer,icp,1,mpi_comm_world,ierr)
278 call mpi_send(dm(i)%c%nttbr, 1,mpi_integer,icp,1,mpi_comm_world,ierr)
279 call mpi_send(dm(i)%c%nrows, 1,mpi_integer,icp,1,mpi_comm_world,ierr)
280 call mpi_send(dm(i)%c%ncols, 1,mpi_integer,icp,1,mpi_comm_world,ierr)
281
282 call mpi_send(dm(i)%c%irow, dm(i)%c%nttbr, mpi_integer,icp,1,mpi_comm_world,ierr)
283 call mpi_send(dm(i)%c%jcol, dm(i)%c%nttbr, mpi_integer,icp,1,mpi_comm_world,ierr)
284 call mpi_send(dm(i)%c%val, dm(i)%c%nttbr*dm(i)%c%ndeg*dm(i)%c%ndeg, mpi_real8,icp,1,mpi_comm_world,istatus,ierr)
285 end do
286
287 !call MPI_BARRIER(MPI_COMM_WORLD, ierr)
288
289 ! clean up
290 do i=1,m_pds_procinfo%nchildren
291 deallocate(dm(i)%a%irow,dm(i)%a%jcol,dm(i)%a%val)
292 deallocate(dm(i)%c%irow,dm(i)%c%jcol,dm(i)%c%val)
293 end do
294
295 call elapout('sp_direct_parent: end send matrix') !elap
296
297 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
298 !
299 ! STEP13: Receive D region and update D matrix as D' = D - D1' - D2' ...
300 !
301 ! D is receive as dens matrix, which format is given in s3um2() in serial solver.
302 ! to decompose this dens D matrix, use s3um3() on parent.
303
304 call elapout('sp_direct_parent: receive D matrix element') !elap
305 allocate(diagbuf(ndegt, neqns_d), stat=ierr)
306 if(ierr .ne. 0) then
307 call errtrp('stop due to allocation error.')
308 end if
309
310 do k=1,m_pds_procinfo%nchildren
311 icp=m_pds_procinfo%ichildren(k)
312 call mpi_recv(nspdsln, 1,mpi_integer,icp,1,mpi_comm_world,istatus,ierr)
313 allocate(spdslnidx(nspdsln),spdslnval(nndeg,nspdsln),stat=ierr)
314 if(ierr .ne. 0) then
315 call errtrp('stop due to allocation error.')
316 end if
317 call mpi_recv(spdslnidx, nspdsln, mpi_integer,icp,1,mpi_comm_world,istatus,ierr)
318 call mpi_recv(spdslnval, nspdsln*nndeg,mpi_real8,icp,1,mpi_comm_world,istatus,ierr)
319 call mpi_recv(diagbuf, neqns_d*ndegt,mpi_real8,icp,1,mpi_comm_world,istatus,ierr)
320
321 ! off diagonal
322 do i=1,nspdsln
323 dsln(:,spdslnidx(i)) = dsln(:,spdslnidx(i)) + spdslnval(:,i) ! because of child process dsln is already substructed in s3um2()
324 end do
325
326 ! diagonal
327 do i=1,neqns_d
328 do j=1,ndegt
329 diag(j,i) = diag(j,i) + diagbuf(j,i)
330 end do
331 end do
332 deallocate(spdslnidx, spdslnval)
333 end do
334 deallocate(diagbuf)
335 call elapout('sp_direct_parent: end receive D matrix element') !elap
336
337 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
338 !
339 ! STEP13: ! LDU decompose dens D
340 !
341 call elapout('sp_direct_parent: LDU decompose of D. entering nufct0_parent') !elap
342 call nufct0_parent(dsln, diag, neqns_d, ndeg)
343 call elapout('sp_direct_parent: exit nufct0_parent') !elap
344
345
346 nusol_ready = .true.
347 end if
348
349 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
350 !
351 ! STEP2X: Solve Ax=b0
352 !
353 ! STEP21: divide right hand side b0 and send it to child in charge.
354 !
355 ! right hand side vector b is reordered as follows:
356 ! b(1:neqns_d)=for parent process
357 ! b(neqns_d+1:number of b in dm(1)) for dm(1)
358 ! b(end of b1:number of b in dm(2)) for dm(2)
359 ! ..
360 ! b(end of b_n-1:neqns_t) for dm(ndm)
361
362 ! set right hand side vector (b)
363 allocate(b(ndeg,a0%neqns), stat=ierr)
364 if(ierr .ne. 0) then
365 call errtrp('stop due to allocation error.')
366 end if
367 do i=1,a0%neqns
368 do j=1,ndeg
369 b(j,i)=hecmat%b(ndeg*(i-1)+j)
370 end do
371 end do
372
373 ! for verify
374 allocate(oldb(ndeg,a0%neqns), stat=ierr)
375 if(ierr .ne. 0) then
376 call errtrp('stop due to allocation error.')
377 end if
378 oldb=b
379
380 call reovec(b, iperm_all)
381 do i=1,m_pds_procinfo%nchildren
382 icp=m_pds_procinfo%ichildren(i)
383 call mpi_send(b(1,iofst_dm(i)+1), dm(i)%ndeg*dm(i)%a%neqns, mpi_real8, icp, 1,mpi_comm_world, ierr)
384 end do
385
386 allocate(bd(ndeg,neqns_d), stat=ierr)
387 if(ierr .ne. 0) then
388 call errtrp('stop due to allocation error.')
389 end if
390 bd(:,1:neqns_d) = b(:,1:neqns_d)
391
392 call elapout('sp_direct_parent: end send b') !elap
393
394 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
395 !
396 ! STEP22 forward substitution for D region. and get Y
397 !
398 ! b1, b2... are sended to child processes.
399 ! bd is substituted to D in locally, and results from child processes
400 ! (C1-Y1 substitute, C2-Y2 substitute...) are receive from child processes.
401 ! these value also add for Yd
402 call elapout('sp_direct_parent: begin receive bd') !elap
403 allocate(bdbuf(ndeg,neqns_d), stat=ierr)
404 if(ierr .ne. 0) then
405 call errtrp('stop due to allocation error.')
406 end if
407 bdbuf=0
408 do k=1,m_pds_procinfo%nchildren
409 icp=m_pds_procinfo%ichildren(k)
410 call mpi_recv(bdbuf, ndeg*neqns_d, mpi_real8, icp, 1,mpi_comm_world, istatus, ierr)
411 do i=1,neqns_d
412 do j=1,ndeg
413 bd(j,i) = bd(j,i) + bdbuf(j,i)
414 end do
415 end do
416 end do
417
418 call elapout('sp_direct_parent: end receive bd') !elap
419
420 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
421 !
422 ! STEP22 solve Ax=b for dens matrix, using updated bd
423 !
424 call elapout('sp_direct_parent: begin solve Ax_d=b_d') !elap
425 call nusol0_parent(dsln, diag, bd, neqns_d, ndeg)
426 call elapout('sp_direct_parent: end solve Ax_d=b_d') !elap
427
428
429 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
430 !
431 ! STEP23 send Xd to children
432 !
433 call elapout('sp_direct_parent: begin send Xd') !elap
434 call mpi_bcast(bd, ndeg*neqns_d, mpi_real8, m_pds_procinfo%imp, mpi_comm_world, ierr)
435 call elapout('sp_direct_parent: end send Xd') !elap
436
437 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
438 !
439 ! STEP24 Receive results Xi from children.
440 !
441 call elapout('sp_direct_parent: begin receive X') !elap
442 do k=1,m_pds_procinfo%nchildren
443 icp=m_pds_procinfo%ichildren(k)
444 call mpi_recv(b(1,iofst_dm(k)+1), dm(k)%ndeg*dm(k)%a%neqns, mpi_real8, icp, 1,mpi_comm_world, istatus, ierr)
445 end do
446 b(:,1:neqns_d)=bd(:,:) ! set xd
447 call elapout('sp_direct_parent: end receive X') !elap
448
449 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
450 !
451 ! STEP25 restore final result
452 !
453 ! permtation divided matrix => initial order
454 call elapout('sp_direct_parent: begin permutate X') !elap
455 call reovec(b, iperm_rev)
456 call elapout('sp_direct_parent: end permutate X') !elap
457
458 ! verify result
459 call verif0(a0%neqns, ndeg, a0%nttbr, a0%irow, a0%jcol, a0%val, oldb, b) !verify result oldb will be broken.
460
461 ! set result to FEM data
462 do i=1,a0%neqns
463 do j=1,ndeg
464 hecmat%x(ndeg*(i-1)+j)=b(j,i)
465 end do
466 end do
467
468 call elapout('sp_direct_parent: end solve Ax=b') !elap
469 deallocate(b, bd, bdbuf, oldb)
470 return
471 end subroutine sp_direct_parent
472
473 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
474
475 subroutine geta0(hecMESH,hecMAT, a0)
476
477 implicit none
478
479 type (hecmwst_local_mesh), intent(in) :: hecmesh
480 type (hecmwst_matrix ), intent(in) :: hecmat
481 type (irjc_square_matrix), intent(out) :: a0
482
483 integer(kind=kint) :: i,j,k,l,ierr,numnp,ndof,kk,ntotal
484 integer(kind=kint) :: iis,iie,kki,kkj,ndof2
485
486 numnp = hecmat%NP
487 ndof = hecmesh%n_dof
488 ntotal = numnp*ndof
489
490 !*NUFACT variables
491 a0%neqns = numnp
492 a0%ndeg = ndof
493 a0%nttbr = hecmat%NP+hecmat%NPL !+hecMAT%NPU if unsymmetric
494
495 !*Allocations
496 allocate(a0%irow(a0%nttbr),stat=ierr)
497 allocate(a0%jcol(a0%nttbr),stat=ierr)
498 allocate(a0%val(a0%ndeg*a0%ndeg, a0%nttbr),stat=ierr)
499 if(ierr .ne. 0) then
500 call errtrp('stop due to allocation error.')
501 end if
502
503 kk = 0
504 ndof2 = ndof*ndof
505 do j= 1, numnp
506 !*Diagonal
507 kk = kk + 1
508 a0%irow(kk) = j
509 a0%jcol(kk) = j
510 call vlcpy(a0%val(:,kk),hecmat%D(ndof2*(j-1)+1:ndof2*j),ndof)
511 !*Lower
512 do k= hecmat%indexL(j-1)+1, hecmat%indexL(j)
513 i= hecmat%itemL(k)
514 kk = kk + 1
515 a0%irow(kk) = j
516 a0%jcol(kk) = i
517 call vlcpy(a0%val(:,kk),hecmat%AL(ndof2*(k-1)+1:ndof2*k),ndof)
518 enddo
519 enddo
520
521 return
522
523 contains
524
525 subroutine vlcpy(a,b,n)
526 implicit none
527 real(kind=kreal), intent(out) :: a(:)
528 real(kind=kreal), intent(in) :: b(:)
529 integer(kind=kint), intent(in) :: n
530
531 integer(kind=kint) :: i,j
532
533 do i = 1,n
534 do j = 1,n
535 a((j-1)*n+i) = b((i-1)*n+j) !transpose
536 end do
537 end do
538 return
539 end subroutine vlcpy
540
541 end subroutine geta0
542
543 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
544
545 subroutine sp_direct_child()
546 ! parallel direct solver child process
547
548 implicit none
549
550 include 'mpif.h'
551
552 type (child_matrix) :: cm
553
554 real(kind=kreal), allocatable :: b(:,:) ! (ndeg, neqns)
555 type (dsinfo) :: dsi
556 logical, save :: nusol_ready = .false.
557
558 ! for MPI
559 integer(kind=kint) :: istatus(mpi_status_size)
560 integer(kind=kint) :: imp, ierr
561 integer(kind=kint) :: i,j,k,l,m,n
562
563 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
564
565 call elapout('sp_direct_child: entered') !elap
566
567 imp=m_pds_procinfo%imp
568
569 if (.not. nusol_ready) then
570
571 call elapout('sp_direct_child: waiting matrix from parent via MPI') !elap
572 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
573 !
574 ! STEP01: get a,c from parent
575 ! C matrix is placed below A.
576 !
577 call mpi_recv(cm%a%ndeg, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
578 call mpi_recv(cm%a%neqns, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
579 call mpi_recv(cm%a%nttbr, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
580 allocate(cm%a%irow(cm%a%nttbr), stat=ierr)
581 if(ierr .ne. 0) then
582 call errtrp('stop due to allocation error.')
583 end if
584 allocate(cm%a%jcol(cm%a%nttbr), stat=ierr)
585 if(ierr .ne. 0) then
586 call errtrp('stop due to allocation error.')
587 end if
588 allocate(cm%a%val(cm%a%ndeg*cm%a%ndeg, cm%a%nttbr), stat=ierr)
589 if(ierr .ne. 0) then
590 call errtrp('stop due to allocation error.')
591 end if
592 call mpi_recv(cm%a%irow, cm%a%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
593 call mpi_recv(cm%a%jcol, cm%a%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
594 call mpi_recv(cm%a%val, cm%a%nttbr*cm%a%ndeg*cm%a%ndeg, mpi_real8,imp,1,mpi_comm_world,istatus,ierr)
595
596
597 call mpi_recv(cm%c%ndeg, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
598 call mpi_recv(cm%c%nttbr, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
599 call mpi_recv(cm%c%nrows, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
600 call mpi_recv(cm%c%ncols, 1,mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
601 allocate(cm%c%irow(cm%c%nttbr), stat=ierr)
602 if(ierr .ne. 0) then
603 call errtrp('stop due to allocation error.')
604 end if
605 allocate(cm%c%jcol(cm%c%nttbr), stat=ierr)
606 if(ierr .ne. 0) then
607 call errtrp('stop due to allocation error.')
608 end if
609 allocate(cm%c%val(cm%c%ndeg*cm%c%ndeg, cm%c%nttbr), stat=ierr)
610 if(ierr .ne. 0) then
611 call errtrp('stop due to allocation error.')
612 end if
613 call mpi_recv(cm%c%irow, cm%c%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
614 call mpi_recv(cm%c%jcol, cm%c%nttbr, mpi_integer,imp,1,mpi_comm_world,istatus,ierr)
615 call mpi_recv(cm%c%val, cm%c%nttbr*cm%c%ndeg*cm%c%ndeg, mpi_real8,imp,1,mpi_comm_world,istatus,ierr)
616
617 cm%ndeg = cm%a%ndeg
618 cm%ista_c = cm%a%neqns+1
619 cm%neqns_t = cm%a%neqns + cm%c%nrows
620 call elapout('sp_direct_child: end get matrix from parent via MPI') !elap
621 !call MPI_BARRIER(MPI_COMM_WORLD, ierr)
622
623
624
625 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
626 !
627 ! STEP1x: LDU decompose for given matrix
628 !
629
630 ! set up dsi for allocate matrix array for fill-in
631 call elapout('sp_direct_child: entering matini_para') !elap
632 call matini_para(cm, dsi, ierr)
633 call elapout('sp_direct_child: exit matini_para') !elap
634
635 call elapout('sp_direct_child: entering staij1') !elap
636 ! set real8 value
637 do i=1,cm%a%nttbr
638 call staij1(0, cm%a%irow(i), cm%a%jcol(i), cm%a%val(:,i), dsi, ierr)
639 end do
640 do i=1,cm%c%nttbr
641 ! call staij1(0, cm%c%irow(i)+cm%a%neqns, dsi%iperm(cm%c%jcol(i)), cm%c%val(:,i), dsi, ierr)
642 call staij1(0, cm%c%irow(i)+cm%a%neqns, cm%c%jcol(i), cm%c%val(:,i), dsi, ierr)
643 end do
644 call elapout('sp_direct_child: end staij1') !elap
645
646 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
647 !
648 ! following STEP12-15 will be done in nufct0_child()
649 !
650 ! STEP12: LDU decompose of A (1..nstop-1)
651 ! STEP13: LDU decompose of C (nstop..neqnsA+neqnsd)
652 ! STEP14: update D region.
653 ! STEP15: send D region to parent
654 call elapout('sp_direct_child: entering nufct0_child') !elap
655 call nufct0_child(dsi, ierr)
656 call elapout('sp_direct_child: exit nufct0_child') !elap
657
658 nusol_ready = .true.
659 end if
660
661 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
662 !
663 ! STEP2x: solve ax=b by using forward an backword substitution
664 !
665 ! STEP21: receive b from parent
666 !
667 allocate(b(cm%ndeg, cm%neqns_t), stat=ierr)
668 if(ierr .ne. 0) then
669 call errtrp('stop due to allocation error.')
670 end if
671
672 ! wait for right hand side vector
673 call mpi_recv(b, cm%ndeg*cm%a%neqns, mpi_real8, imp, 1,mpi_comm_world, istatus, ierr)
674
675 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
676 !
677 ! following STEP22-24 is done in nusol0-child()
678 !
679 ! STEP22: forward substitution for A
680 ! STEP23: forward substitution for C and send it (yi) to parent
681 ! STEP24: divide with diagonal matrix
682 ! STEP25: receive xd from parent and do backword substitution
683 call elapout('sp_direct_child: enter nusol0_child') !elap
684 call nusol0_child(b, dsi, ierr)
685 call elapout('sp_direct_child: exit nusol0_child') !elap
686
687 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
688 !
689 ! STEP26: send final result to parent
690 !
691 call elapout('sp_direct_child: begin send result to parent') !elap
692 call mpi_send(b, cm%ndeg*cm%a%neqns, mpi_real8, imp, 1,mpi_comm_world, ierr)
693 call elapout('sp_direct_child: end send result to parent') !elap
694
695 return
696
697 end subroutine sp_direct_child
698
699
700 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
701
702 subroutine initproc()
703 implicit none
704
705 include 'mpif.h'
706
707 integer(kind=kint) :: npe, myid
708 integer(kind=kint) :: ndiv
709 integer(kind=kint) :: ierr
710 integer(kind=kint) :: i,j,k,l
711
712 m_pds_procinfo%isparent=.false.
713 m_pds_procinfo%ischild=.false.
714
715 ! get process number and number of whole processes
716 call mpi_comm_size(mpi_comm_world, npe, ierr)
717 call mpi_comm_rank(mpi_comm_world, myid, ierr)
718
719 m_pds_procinfo%myid = myid
720
721 ndiv=0
722 do
723 if (2**(ndiv + 1) .gt. npe) then
724 exit
725 end if
726 ndiv = ndiv + 1
727 end do
728 m_pds_procinfo%ndiv = ndiv
729
730 if (npe .ne. 2**ndiv + 1) then
731 write(ilog,*) 'Error: please use 2**n+1 (3,5,9,17...) processes for parallel direct solver.'
732 write(6,*) 'Error: please use 2**n+1 (3,5,9,17...) processes for parallel direct solver.'
734 stop
735 end if
736
737 if (myid.eq.0) then
738 write(idbg,*)'parent process.'
739 m_pds_procinfo%isparent=.true.
740 m_pds_procinfo%nchildren=2**ndiv
741 allocate(m_pds_procinfo%ichildren(m_pds_procinfo%nchildren))
742 do i=1, 2**ndiv
743 m_pds_procinfo%ichildren(i)=i
744 end do
745 else
746 write(idbg,*)'child process.'
747 m_pds_procinfo%ischild=.true.
748 m_pds_procinfo%imp=0
749 end if
750
751 return
752 end subroutine initproc
753
754 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
755 subroutine errtrp(mes)
756 character(*) mes
757 write(6,*) 'Error in : process ', m_pds_procinfo%myid
758 write(ilog,*) mes
759
761 stop
762 end subroutine errtrp
763
764 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
765
766 subroutine matini_para(cm,dsi,ir)
767
768 !----------------------------------------------------------------------
769 !
770 ! matini initializes storage for sparse matrix solver.
771 ! this routine is used for both symmetric matrices
772 ! and must be called once at the beginning
773 !
774 ! (i)
775 ! neqns number of unknowns
776 ! nttbr number of non0s, pattern of non-zero elements are
777 ! given like following.
778 ! nonz(A)={(i,j);i=irow(l),j=jcol(l); 1<= l <= nttbr}
779 ! irow
780 ! jcol to define non-zero pattern
781 ! lenv length of the array v (iv)
782 !
783 ! (o)
784 ! dsi matrix informations
785 ! ir return code
786 ! =0 normal
787 ! =-1 non positive index
788 ! =1 too big index
789 ! =10 insufficient storage
790 !
791 !
792 ! stage 10 after initialization
793 ! 20 building up matrix
794 ! 30 after LU decomposition
795 ! 40 after solving
796 !
797 ! # coded by t.arakawa
798 ! # reviced by t.kitayama of Univ. Tokyo on 20071120
799 !
800 !----------------------------------------------------------------------
801
802 implicit none
803
804 type(child_matrix), intent(in) :: cm
805 type(dsinfo), intent(out) :: dsi
806 integer(kind=kint), intent(out) :: ir
807
808 integer(kind=kint), pointer :: irow_a(:), jcol_a(:)
809 integer(kind=kint), pointer :: irow_c(:), jcol_c(:)
810
811 integer(kind=kint), pointer :: ia(:) ! in stiaja() neqns+2
812 integer(kind=kint), pointer :: ja(:) ! in stiaja() 2*nttbr
813 integer(kind=kint), pointer :: jcpt(:) ! in stsmat() 2*nttbr
814 integer(kind=kint), pointer :: jcolno(:) ! in stsmat() 2*nttbr
815
816 integer(kind=kint), pointer :: iperm_a(:)
817 integer(kind=kint), pointer :: invp_a(:)
818
819 integer(kind=kint), pointer :: xlnzr_a(:)
820 integer(kind=kint), pointer :: colno_a(:)
821
822 integer(kind=kint), pointer :: xlnzr_c(:)
823 integer(kind=kint), pointer :: colno_c(:)
824
825
826 integer(kind=kint), pointer :: adjncy(:) ! in genqmd() 2*nttbr
827 integer(kind=kint), pointer :: qlink(:) ! in genqmd() neqne+2
828 integer(kind=kint), pointer :: qsize(:) ! in genqmd() neqne+2
829 integer(kind=kint), pointer :: nbrhd(:) ! in genqmd() neqne+2
830 integer(kind=kint), pointer :: rchset(:) ! in genqmd() neqne+2
831
832 integer(kind=kint), pointer :: cstr(:)
833
834 integer(kind=kint), pointer :: adjt(:) ! in rotate() neqne+2
835 integer(kind=kint), pointer :: anc(:) ! in rotate() neqne+2
836
837 integer(kind=kint), pointer :: lwk3arr(:)
838 integer(kind=kint), pointer :: lwk2arr(:)
839 integer(kind=kint), pointer :: lwk1arr(:)
840 integer(kind=kint), pointer :: lbtreearr(:,:) ! genbtq() (2,neqns+1)
841 integer(kind=kint), pointer :: lleafarr(:)
842 integer(kind=kint), pointer :: lxleafarr(:)
843 integer(kind=kint), pointer :: ladparr(:)
844 integer(kind=kint), pointer :: lpordrarr(:)
845
846 integer(kind=kint) :: neqns_a, nttbr_a, neqns_a1, nstop, neqns_t, neqns_d, nttbr_c, ndeg
847 integer(kind=kint) :: lncol_a, lncol_c
848 integer(kind=kint) :: neqnsz, nofsub, izz, izz0, lnleaf ! dummy variables
849 integer(kind=kint) :: ir1
850 integer(kind=kint) :: i, j, k , ipass, ks, ke, ierr
851
852 ndeg = cm%ndeg
853
854 neqns_t = cm%neqns_t
855 neqns_d = cm%c%nrows
856
857 neqns_a = cm%a%neqns
858 nttbr_a = cm%a%nttbr
859 irow_a => cm%a%irow
860 jcol_a => cm%a%jcol
861
862 nttbr_c = cm%c%nttbr
863 irow_c => cm%c%irow
864 jcol_c => cm%c%jcol
865
866 dsi%neqns=neqns_t ! because direct solver treat A + C as one matrix.
867 dsi%ndeg=ndeg
868
869 neqns_a1=neqns_a+2
870 ir=0
871 ierr=0
872 izz0 = 0.0d0
873 !
874 ! set z pivot
875 !
876 allocate(dsi%zpiv(neqns_a), stat=ierr)
877 if(ierr .ne. 0) then
878 call errtrp('stop due to allocation error.')
879 end if
880 call zpivot(neqns_a,neqnsz,nttbr_a,jcol_a,irow_a,dsi%zpiv,ir1)
881 if(ir1.ne.0) then
882 ir=ir1
883 goto 1000
884 endif
885 !
886 ! build jcpt,jcolno
887 !
888 allocate(jcpt(2*nttbr_a), jcolno(2*nttbr_a), stat=ierr)
889 if(ierr .ne. 0) then
890 call errtrp('stop due to allocation error.')
891 end if
892 call stsmat(neqns_a,nttbr_a,irow_a,jcol_a,jcpt,jcolno)
893 !
894 ! build ia,ja
895 !
896 allocate(ia(neqns_a1), ja(2*nttbr_a), stat=ierr)
897 if(ierr .ne. 0) then
898 call errtrp('stop due to allocation error.')
899 end if
900 call stiaja(neqns_a, neqns_a,ia,ja,jcpt,jcolno)
901 !
902 ! get permutation vector iperm,invp
903 !
904
905 ! setup identity permtation for C matrix
906 allocate(iperm_a(neqns_a), invp_a(neqns_a), stat=ierr)
907 if(ierr .ne. 0) then
908 call errtrp('stop due to allocation error.')
909 end if
910 call idntty(neqns_a,invp_a,iperm_a)
911
912 ! reorder A matrix
913 allocate(adjncy(2*nttbr_a),qlink(neqns_a1),qsize(neqns_a1),nbrhd(neqns_a1),rchset(neqns_a1), stat=ierr)
914 if(ierr .ne. 0) then
915 call errtrp('stop due to allocation error.')
916 end if
917 allocate(lwk2arr(neqns_a1),lwk1arr(neqns_a1), stat=ierr)
918 if(ierr .ne. 0) then
919 call errtrp('stop due to allocation error.')
920 end if
921 call genqmd(neqns_a,ia,ja,iperm_a,invp_a,lwk1arr,lwk2arr,rchset,nbrhd,qsize,qlink,nofsub,adjncy)
922 deallocate(adjncy, qlink, qsize, nbrhd, rchset)
923
924 ! set ia, ja
925 call stiaja(neqns_a, neqns_a, ia, ja, jcpt,jcolno)
926
927 ! build up the parent vector parent vector will be saved in
928 ! work2 for a while
929 allocate(cstr(neqns_a1),adjt(neqns_a1), stat=ierr)
930 if(ierr .ne. 0) then
931 call errtrp('stop due to allocation error.')
932 end if
933 10 continue
934 call genpaq(ia,ja,invp_a,iperm_a,lwk2arr,neqns_a,cstr)
935
936 ! build up the binary tree
937 allocate (lbtreearr(2,neqns_a1), stat=ierr)
938 if(ierr .ne. 0) then
939 call errtrp('stop due to allocation error.')
940 end if
941 call genbtq(ia, ja, invp_a, iperm_a,lwk2arr,lbtreearr,dsi%zpiv,izz,neqns_a)
942
943 ! rotate the binary tree to avoid a zero pivot
944 if(izz.eq.0) goto 20
945 if(izz0.eq.0) izz0=izz
946 if(izz0.ne.izz) goto 30
947 call rotate(ia, ja, invp_a, iperm_a, lwk2arr,lbtreearr,izz,neqns_a,anc,adjt,ir1)
948 goto 10
949 30 continue
950 call bringu(dsi%zpiv,iperm_a, invp_a, lwk2arr,izz,neqns_a,ir1)
951 goto 10
952
953 ! post ordering
954 20 continue
955 allocate(lwk3arr(0:neqns_a1),lpordrarr(neqns_a1),dsi%parent(neqns_a1), dsi%nch(neqns_a1), stat=ierr)
956 if(ierr .ne. 0) then
957 call errtrp('stop due to allocation error.')
958 end if
959 call posord(dsi%parent,lbtreearr,invp_a,iperm_a,lpordrarr,dsi%nch,neqns_a,lwk1arr,lwk2arr,lwk3arr)
960
961 ! generate skelton graph
962 allocate(lleafarr(nttbr_a),lxleafarr(neqns_a1),ladparr(neqns_a1), stat=ierr)
963 if(ierr .ne. 0) then
964 call errtrp('stop due to allocation error.')
965 end if
966 call gnleaf(ia, ja, invp_a, iperm_a, lpordrarr,dsi%nch,ladparr,lxleafarr,lleafarr,neqns_a,lnleaf)
967
968 ! build up xlnzr,colno (this is the symbolic fct.)
969 nstop = cm%ista_c
970 call countclno(dsi%parent, lxleafarr, lleafarr, neqns_a, nstop, lncol_a, ir1) ! only for A
971 allocate(colno_a(lncol_a),xlnzr_a(neqns_a1), stat=ierr)
972 if(ierr .ne. 0) then
973 call errtrp('stop due to allocation error.')
974 end if
975 call gnclno(dsi%parent,lpordrarr,lxleafarr,lleafarr,xlnzr_a, colno_a, neqns_a, nstop,lncol_a,ir1) ! only for A
976 if(ir1.ne.0) then
977 ir=10
978 goto 1000
979 endif
980
981 ! do symbolic LDU decomposition for C region.
982 call ldudecomposec(xlnzr_a,colno_a,invp_a,iperm_a, ndeg, nttbr_c, irow_c, &
983 jcol_c, cm%c%ncols, cm%c%nrows, xlnzr_c, colno_c, lncol_c)
984
985 ! set calculated informations to dsi.
986 allocate(dsi%xlnzr(neqns_t + 1), stat=ierr)
987 if(ierr .ne. 0) then
988 call errtrp('stop due to allocation error.')
989 end if
990 dsi%xlnzr(1:neqns_a)=xlnzr_a(:)
991 dsi%xlnzr(neqns_a+1:neqns_t+1)=xlnzr_c(:)+xlnzr_a(neqns_a+1)-1
992
993 dsi%lncol=lncol_a + lncol_c
994 allocate(dsi%colno(lncol_a + lncol_c), stat=ierr)
995 if(ierr .ne. 0) then
996 call errtrp('stop due to allocation error.')
997 end if
998 dsi%colno(1:lncol_a)=colno_a(:)
999 dsi%colno(lncol_a+1:lncol_a+lncol_c)=colno_c(:)
1000
1001 allocate(dsi%invp(neqns_t), dsi%iperm(neqns_t), stat=ierr)
1002 if(ierr .ne. 0) then
1003 call errtrp('stop due to allocation error.')
1004 end if
1005 dsi%invp(1:neqns_a)=invp_a(1:neqns_a)
1006 dsi%iperm(1:neqns_a)=iperm_a(1:neqns_a)
1007 do i=neqns_a+1,neqns_t
1008 dsi%invp(i)=i
1009 dsi%iperm(i)=i
1010 end do
1011
1012 deallocate(xlnzr_a, colno_a, xlnzr_c, colno_c, invp_a, iperm_a)
1013
1014 dsi%nstop=nstop
1015 dsi%stage=10
1016 1000 continue
1017
1018 return
1019 end subroutine matini_para
1020
1021 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1022
1023 subroutine nufct0_child(dsi,ir)
1024
1025 implicit none
1026 type(dsinfo), intent(inout) :: dsi
1027 integer(kind=kint), intent(out) :: ir
1028 !
1029 ! this performs Cholesky factorization
1030 !
1031
1032 if(dsi%stage.ne.20) then
1033 ir=40
1034 goto 1000
1035 else
1036 ir=0
1037 endif
1038 if(dsi%ndeg.eq.1) then
1039 call nufct1_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%neqns,dsi%parent,dsi%nch,dsi%nstop,ir)
1040 else if(dsi%ndeg.eq.2) then
1041 call nufct2_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%neqns,dsi%parent,dsi%nch,dsi%nstop,ir)
1042 else if(dsi%ndeg.eq.3) then
1043 call nufct3_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%neqns,dsi%parent,dsi%nch,dsi%nstop,ir)
1044 ! else if(dsi%ndeg.eq.6) then !TODO implement it
1045 ! call nufct6_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%neqns,dsi%parent,dsi%nch,dsi%nstop,ir)
1046 else
1047 call nufctx_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%neqns,dsi%parent,dsi%nch,dsi%nstop,dsi%ndeg,ir)
1048 end if
1049
1050 dsi%stage=30
1051 1000 continue
1052 return
1053 end subroutine nufct0_child
1054
1055 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1056
1057 subroutine nufct1_child(xlnzr,colno,zln,diag,neqns,parent,nch,nstop,ir)
1058
1059 implicit none
1060 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),parent(:)
1061 integer(kind=kint), intent(in) :: neqns, nstop, ir
1062 integer(kind=kint), intent(out) :: nch(:)
1063 real(kind=kreal), intent(out) :: zln(:,:),diag(:,:) !zln(1,:), diag(1,:)
1064
1065 integer(kind=kint) :: neqns_c
1066 integer(kind=kint) :: i,j,k,l, ic,ierr,imp
1067 integer(kind=kint) :: nspdsln
1068 integer(kind=kint), pointer :: spdslnidx(:)
1069 real(kind=kreal), pointer :: spdslnval(:,:)
1070 include 'mpif.h'
1071
1072 !----------------------------------------------------------------------
1073 !
1074 ! nufct1 performs cholesky factorization in row order for ndeg=1
1075 !
1076 ! (i) xlnzr,colno,zln,diag
1077 ! symbolicaly factorized
1078 !
1079 ! (o) zln,diag,dsln
1080 !
1081 ! #coded by t.arakawa
1082 !
1083 !----------------------------------------------------------------------
1084
1085 !
1086 ! phase I
1087 ! LDU decompose of A (1..nstop-1)
1088 !
1089 diag(1,1)=1.0d0/diag(1,1)
1090 l=parent(1)
1091 nch(l)=nch(l)-1
1092 nch(1)=-1
1093 do 100 ic=2,nstop-1
1094 call sum(ic,xlnzr,colno,zln(1,:),diag(1,:),nch,parent,neqns)
1095 100 continue
1096 !
1097 ! phase II
1098 ! LDU decompose of C (nstop..neqnsA+neqnsd)
1099 !
1100 do 200 ic=nstop,neqns
1101 call sum1(ic,xlnzr,colno,zln(1,:),diag(1,:),parent,neqns)
1102 200 continue
1103 !
1104 ! phase III
1105 ! Update D region.
1106 !
1107
1108 ! clear dummy diagonal value for D region
1109 do i=nstop,neqns
1110 diag(:,i)=0.0
1111 end do
1112
1113 neqns_c = neqns - nstop + 1
1114 call sum2_child(neqns,nstop,xlnzr,colno,zln(1,:),diag(1,:),spdslnidx,spdslnval,nspdsln)
1115 ! send D region to parent
1116 imp = m_pds_procinfo%imp
1117 call mpi_send(nspdsln, 1,mpi_integer,imp,1,mpi_comm_world,ierr)
1118 call mpi_send(spdslnidx, nspdsln,mpi_integer,imp,1,mpi_comm_world,ierr)
1119 call mpi_send(spdslnval, nspdsln,mpi_real8,imp,1,mpi_comm_world,ierr)
1120 call mpi_send(diag(1,nstop), neqns_c,mpi_real8,imp,1,mpi_comm_world,ierr)
1121 return
1122 end subroutine nufct1_child
1123
1124 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1125
1126 subroutine nufct2_child(xlnzr, colno, zln, diag, neqns, parent, nch, nstop, ir)
1127
1128 implicit none
1129
1130 integer(kind=kint), intent(in) :: xlnzr(:), colno(:), parent(:)
1131 integer(kind=kint), intent(out) :: nch(:)
1132 real(kind=kreal), intent(out) :: zln(:,:), diag(:,:) !zln(6,*), diag(3,*)
1133 integer(kind=kint), intent(in) :: neqns, nstop
1134 integer(kind=kint), intent(out) :: ir
1135
1136 integer(kind=kint) :: i,j,k,l, ic, imp, ierr, neqns_c
1137 integer(kind=kint) :: nspdsln
1138 integer(kind=kint), pointer :: spdslnidx(:)
1139 real(kind=kreal), pointer :: spdslnval(:,:)
1140 include 'mpif.h'
1141
1142 !----------------------------------------------------------------------
1143 !
1144 ! nufct2 performs cholesky factorization in row order for ndeg=2
1145 !
1146 ! (i) xlnzr,colno,zln,diag
1147 ! symbolicaly factorized
1148 !
1149 ! (o) zln,diag,dsln
1150 !
1151 ! #coded by t.arakawa
1152 !
1153 !----------------------------------------------------------------------
1154
1155
1156 ! For parallel calculation, factorization for A and C region
1157 ! will be done.
1158 !
1159 ! Creation for D region also done.
1160 !
1161 ! Factorization for D region is omitted.
1162
1163 !
1164 ! phase I
1165 ! LDU decompose of A (1..nstop-1)
1166 !
1167 call elapout('nufct2_child: begin phase I LDU decompose of A') !elap
1168 if(nstop.gt.1) call inv2(diag(:,1),ir)
1169 l=parent(1)
1170 nch(l)=nch(l)-1
1171 nch(1)=-1
1172 do ic=2,nstop-1
1173 call s2um(ic,xlnzr,colno,zln,diag,nch,parent,neqns)
1174 end do
1175 !
1176 ! phase II
1177 ! LDU decompose of C (nstop..neqnsA+neqnsd)
1178 !
1179 call elapout('nufct2_child: begin phase II LDU decompose of C') !elap
1180 do ic=nstop,neqns
1181 call s2um1(ic,xlnzr,colno,zln,diag,nch,parent,neqns)
1182 end do
1183
1184 !
1185 ! phase III
1186 ! Update D region.
1187 !
1188
1189 call elapout('nufct2_child: begin phase III update D region') !elap
1190
1191 ! clear dummy diagonal value for D region ! currently dummy value setting was not done
1192 do i=nstop,neqns
1193 diag(:,i)=0.0
1194 end do
1195
1196 neqns_c = neqns - nstop + 1
1197 call s2um2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln)
1198
1199 ! send D region to parent
1200 imp = m_pds_procinfo%imp
1201 call mpi_send(nspdsln, 1,mpi_integer,imp,1,mpi_comm_world,ierr)
1202 call mpi_send(spdslnidx, nspdsln,mpi_integer,imp,1,mpi_comm_world,ierr)
1203 call mpi_send(spdslnval, nspdsln*4,mpi_real8,imp,1,mpi_comm_world,ierr)
1204 call mpi_send(diag(1,nstop), neqns_c*3,mpi_real8,imp,1,mpi_comm_world,ierr)
1205
1206
1207 return
1208
1209 end subroutine nufct2_child
1210
1211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1212
1213 subroutine nufct3_child(xlnzr,colno,zln,diag,neqns,parent,nch,nstop,ir)
1214
1215 implicit none
1216
1217 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),parent(:)
1218 integer(kind=kint), intent(out) :: nch(:)
1219 real(kind=kreal), intent(out) :: zln(:,:),diag(:,:) !zln(9,*),diag(6,*)
1220 integer(kind=kint), intent(in) :: neqns, nstop, ir
1221
1222 integer(kind=kint) :: i,j,k,l, ic, imp, ierr, neqns_c
1223 integer(kind=kint) :: nspdsln
1224 integer(kind=kint), pointer :: spdslnidx(:)
1225 real(kind=kreal), pointer :: spdslnval(:,:)
1226 include 'mpif.h'
1227 !
1228 !----------------------------------------------------------------------
1229 !
1230 ! nufct3 performs cholesky factorization in row order for ndeg=3
1231 !
1232 ! (i) xlnzr,colno,zln,diag
1233 ! symbolicaly factorized
1234 !
1235 ! (o) zln,diag,dsln
1236 !
1237 ! #coded by t.arakawa
1238 !
1239 !----------------------------------------------------------------------
1240
1241 !
1242 ! For parallel calculation, factorization for A and C region
1243 ! will be done.
1244 !
1245 ! Creation for D region also done.
1246 !
1247 ! Factorization for D region is omitted.
1248 !
1249
1250 !
1251 ! phase I
1252 ! LDU decompose of A (1..nstop-1)
1253 !
1254 call elapout('nufct3_child: begin phase I LDU decompose of A') !elap
1255 if(nstop.gt.1) call inv3(diag(:,1),ir)
1256 l=parent(1)
1257 nch(l)=nch(l)-1
1258 nch(1)=-1
1259 do 100 ic=2,nstop-1
1260 call s3um(ic,xlnzr,colno,zln,diag,nch,parent,neqns)
1261 100 continue
1262 !
1263 ! phase II
1264 ! LDU decompose of C (nstop..neqnsA+neqnsd)
1265 !
1266 call elapout('nufct3_child: begin phase II LDU decompose of C') !elap
1267 do 200 ic=nstop,neqns
1268 call s3um1(ic,xlnzr,colno,zln,diag,nch,parent,neqns)
1269 200 continue
1270
1271 !
1272 ! phase III
1273 ! Update D region.
1274 !
1275
1276 call elapout('nufct3_child: begin phase III update D region') !elap
1277
1278 ! clear dummy diagonal value for D region ! currently dummy value setting was not done
1279 do i=nstop,neqns
1280 diag(:,i)=0.0
1281 end do
1282
1283 neqns_c = neqns - nstop + 1
1284 call s3um2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln)
1285
1286 ! send D region to parent
1287 imp = m_pds_procinfo%imp
1288 call mpi_send(nspdsln, 1,mpi_integer,imp,1,mpi_comm_world,ierr)
1289 call mpi_send(spdslnidx, nspdsln,mpi_integer,imp,1,mpi_comm_world,ierr)
1290 call mpi_send(spdslnval, nspdsln*9,mpi_real8,imp,1,mpi_comm_world,ierr)
1291 call mpi_send(diag(1,nstop), neqns_c*6,mpi_real8,imp,1,mpi_comm_world,ierr)
1292
1293 return
1294 end subroutine nufct3_child
1295
1296 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1297
1298 !subroutine nufct6_child(xlnzr, colno, zln, diag, neqns, parent, nch, nstop, ir)
1299 !TODO caution! currently this routine is not implemented because of s6um1() is not implemented.
1300 !
1301 !implicit none
1302 !
1303 !integer :: xlnzr(:), colno(:), parent(:), nch(:)
1304 !real(8) :: zln(:,:), diag(:,:) !zln(36,*), diag(21,*)
1305 !
1306 !integer :: neqns, nstop, ir, imp, ierr
1307 !
1308 !integer :: neqns_c
1309 !integer :: l, ic
1310 !integer :: i
1311 !
1312 !real(8), allocatable :: dsln(:,:)
1313 !include 'mpif.h'
1314 !
1315 !!----------------------------------------------------------------------
1316 !!
1317 !! nufct6 performs cholesky factorization in row order for ndeg=3
1318 !!
1319 !! (i) xlnzr,colno,zln,diag
1320 !! symbolicaly factorized
1321 !!
1322 !! (o) zln,diag,dsln
1323 !!
1324 !! #coded by t.arakawa
1325 !!
1326 !!----------------------------------------------------------------------
1327 !
1328 !
1329 !! For parallel calculation, factorization for A and C region
1330 !! will be done.
1331 !!
1332 !! Creation for D region also done.
1333 !!
1334 !! Factorization for D region is omitted.
1335 !!
1336 !
1337 !!
1338 !! phase I
1339 !! LDU decompose of A (1..nstop-1)
1340 !!
1341 !call elapout('nufct6_child: begin phase I LDU decompose of A') !elap
1342 !if(nstop.gt.1) call inv6(diag(:,1),ir)
1343 !l=parent(1)
1344 !nch(l)=nch(l)-1
1345 !nch(1)=-1
1346 !do ic=2,nstop-1
1347 ! call s6um(ic,xlnzr,colno,zln,diag,nch,parent,neqns)
1348 !end do
1349 !
1350 !!
1351 !! phase II
1352 !! LDU decompose of C (nstop..neqnsA+neqnsd)
1353 !!
1354 !call elapout('nufct6_child: begin phase II LDU decompose of C') !elap
1355 !do ic=nstop,neqns
1356 ! call s6um1(ic,xlnzr,colno,zln,diag,nch,parent,neqns)
1357 !end do
1358 !end subroutine nufct6_child
1359
1360 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1361
1362 subroutine nufctx_child(xlnzr,colno,zln,diag,neqns,parent,nch,nstop,ndeg,ir)
1363
1364 implicit none
1365
1366 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),parent(:)
1367 integer(kind=kint), intent(out) :: nch(:)
1368 real(kind=kreal), intent(out) :: zln(:,:),diag(:,:) !zln(9,*),diag(6,*)
1369 integer(kind=kint), intent(in) :: neqns, nstop, ndeg
1370 integer(kind=kint), intent(out) :: ir
1371
1372 integer(kind=kint) :: i,j,k,l, ic, neqns_c, ndeg2, ndegl, imp, ierr
1373 integer(kind=kint) :: nspdsln
1374 integer(kind=kint), pointer :: spdslnidx(:)
1375 real(kind=kreal), pointer :: spdslnval(:,:)
1376 include 'mpif.h'
1377 !----------------------------------------------------------------------
1378 !
1379 ! nufctx performs cholesky factorization in row order for every ndeg
1380 !
1381 ! (i) xlnzr,colno,zln,diag
1382 ! symbolicaly factorized
1383 !
1384 ! (o) zln,diag,dsln
1385 !
1386 ! #coded by t.arakawa
1387 !
1388 !----------------------------------------------------------------------
1389
1390 !
1391 ! For parallel calculation, factorization for A and C region
1392 ! will be done.
1393 !
1394 ! Creation for D region also done.
1395 !
1396 ! Factorization for D region is omitted.
1397 !
1398
1399 ndeg2=ndeg*ndeg
1400 ndegl=(ndeg+1)*ndeg/2
1401 !
1402 ! phase I
1403 ! LDU decompose of A (1..nstop-1)
1404 if(nstop.gt.1) call invx(diag,ndeg,ir)
1405 l=parent(1)
1406 nch(l)=nch(l)-1
1407 nch(1)=-1
1408 do ic=2,nstop-1
1409 call sxum(ic,xlnzr,colno,zln,diag,nch,parent,neqns,ndeg,ndegl)
1410 end do
1411
1412 !
1413 ! phase II
1414 ! LDU decompose of C (nstop..neqnsA+neqnsd)
1415 !
1416 do ic=nstop,neqns
1417 call sxum1(ic,xlnzr,colno,zln,diag,nch,parent,neqns,ndeg,ndegl)
1418 end do
1419
1420 !
1421 ! phase III
1422 ! Update D region.
1423 !
1424 neqns_c = neqns - nstop + 1
1425 call sxum2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln,ndeg,ndegl)
1426
1427 ! send D region to parent
1428 imp = m_pds_procinfo%imp
1429 call mpi_send(nspdsln, 1,mpi_integer,imp,1,mpi_comm_world,ierr)
1430 call mpi_send(spdslnidx, nspdsln,mpi_integer,imp,1,mpi_comm_world,ierr)
1431 call mpi_send(spdslnval, nspdsln*ndeg2,mpi_real8,imp,1,mpi_comm_world,ierr)
1432 call mpi_send(diag(1,nstop), neqns_c*ndegl,mpi_real8,imp,1,mpi_comm_world,ierr)
1433 return
1434 end subroutine nufctx_child
1435
1436 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1437
1438 subroutine nusol0_child(b,dsi,ir)
1439
1440 !----------------------------------------------------------------------
1441 !
1442 ! this performs forward elimination and backward substitution
1443 !
1444 ! (i/o)
1445 ! b on entry right hand side vector
1446 ! on exit solution vector
1447 !
1448 ! #coded by t.arakawa
1449 !
1450 !----------------------------------------------------------------------
1451
1452 implicit none
1453
1454 real(kind=kreal), intent(inout) :: b(:,:)
1455 type(dsinfo), intent(inout) :: dsi
1456 integer(kind=kint), intent(out) :: ir
1457
1458 integer(kind=kint) :: neqns, nstop, ndeg
1459
1460 if(dsi%stage.ne.30 .and. dsi%stage.ne.40) then
1461 ir=50
1462 goto 1000
1463 else
1464 ir=0
1465 end if
1466 if(dsi%ndeg.eq.1) then
1467 call nusol1_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%iperm,b,dsi%neqns,dsi%nstop)
1468 else if(dsi%ndeg .eq. 2) then
1469 call nusol2_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%iperm,b,dsi%neqns,dsi%nstop)
1470 else if(dsi%ndeg .eq. 3) then
1471 call nusol3_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%iperm,b,dsi%neqns,dsi%nstop)!org
1472 ! else if(dsi%ndeg .eq. 6) then !TODO implement it
1473 ! call nusol6_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%iperm,b,dsi%neqns,dsi%nstop)
1474 else
1475 call nusolx_child(dsi%xlnzr,dsi%colno,dsi%zln,dsi%diag,dsi%iperm,b,dsi%neqns,dsi%nstop,dsi%ndeg)
1476 endif
1477 dsi%stage=40
1478 1000 continue
1479 return
1480 end subroutine nusol0_child
1481
1482 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1483
1484 subroutine nusol1_child(xlnzr,colno,zln,diag,iperm,b,neqns,nstop)
1485
1486 implicit none
1487
1488 integer(kind=kint), intent(in) :: xlnzr(:), colno(:), iperm(:)
1489 real(kind=kreal), intent(in) :: zln(:,:), diag(:,:)
1490 real(kind=kreal), intent(inout) :: b(:,:)
1491 integer(kind=kint), intent(in) :: neqns, nstop
1492
1493 integer(kind=kint) :: neqns_a, neqns_c
1494 integer(kind=kint) :: k, ks, ke, i, j, imp, ierr
1495 real(kind=kreal), allocatable :: wk(:), wk_d(:)
1496
1497 include 'mpif.h'
1498
1499 ! forward
1500
1501 ! now nstop is begining point of C
1502 neqns_a = nstop - 1
1503 neqns_c = neqns - nstop + 1
1504
1505 allocate(wk(neqns), stat=ierr)
1506 if(ierr .ne. 0) then
1507 call errtrp('stop due to allocation error.')
1508 end if
1509 wk = 0
1510
1511 do 10 i=1,neqns_a
1512 wk(i)=b(1,iperm(i))
1513 10 continue
1514
1515 ! STEP22: forward substitution for A
1516 do 100 i=1,neqns_a
1517 ks=xlnzr(i)
1518 ke=xlnzr(i+1)-1
1519 if(ke.lt.ks) goto 110
1520 wk(i)=wk(i)-spdot2(wk,zln(1,:),colno,ks,ke)
1521 110 continue
1522 100 continue
1523
1524 ! STEP23: forward substitution for C and send it (yi) to parent
1525 allocate(wk_d(nstop:neqns), stat=ierr)
1526 if(ierr .ne. 0) then
1527 call errtrp('stop due to allocation error.')
1528 end if
1529 wk_d=0
1530
1531 do 101 i=nstop,neqns
1532 ks=xlnzr(i)
1533 ke=xlnzr(i+1)-1
1534 if(ke.lt.ks) goto 111
1535 wk_d(i)=wk_d(i)-spdot2(wk,zln(1,:),colno,ks,ke)
1536 111 continue
1537 101 continue
1538 imp = m_pds_procinfo%imp
1539 call mpi_send(wk_d, neqns_c, mpi_real8, imp, 1,mpi_comm_world, ierr)
1540
1541 ! STEP24: divide with diagonal matrix
1542 do 120 i=1,neqns
1543 wk(i)=wk(i)*diag(1,i)
1544 120 continue
1545
1546 ! STEP25: receive xd from parent and do backword substitution
1547 call mpi_bcast(wk_d, neqns_c, mpi_real8, imp, mpi_comm_world, ierr)
1548
1549 wk(nstop:neqns)=wk_d(nstop:neqns)
1550
1551 do 200 i=neqns,1,-1
1552 ks=xlnzr(i)
1553 ke=xlnzr(i+1)-1
1554 if(ke.lt.ks) goto 200
1555
1556 do 210 k=ks,ke
1557 j=colno(k)
1558 wk(j)=wk(j)-wk(i)*zln(1,k)
1559 210 continue
1560 200 continue
1561
1562 ! permutaion
1563 do 300 i=1,neqns
1564 b(1,iperm(i))=wk(i)
1565 300 continue
1566 return
1567
1568 end subroutine nusol1_child
1569
1570 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1571
1572 subroutine nusol2_child(xlnzr, colno, zln, diag, iperm, b, neqns, nstop)
1573 ! perform forward substitution for sparse matrix,
1574 ! send bd to parent and receive xd from parent,
1575 ! backword substitution using xd ad send final result x to parent.
1576
1577 use m_elap
1578
1579 implicit none
1580
1581 integer(kind=kint), intent(in) :: xlnzr(:), colno(:), iperm(:)
1582 real(kind=kreal), intent(in) :: zln(:,:), diag(:,:) !zln(4,*), diag(3,*), b(2,*)
1583 real(kind=kreal), intent(inout) :: b(:,:) ! b(2,*)
1584
1585 real(kind=kreal), allocatable :: wk(:,:), wk_d(:,:)
1586 integer(kind=kint) :: neqns_c, neqns_a, nstop, neqns, ks, ke
1587 integer(kind=kint) :: i, j, k, l, imp, ierr
1588
1589 include 'mpif.h'
1590
1591 !now nstop is begining point of C
1592 neqns_a = nstop - 1
1593 neqns_c = neqns - nstop + 1
1594
1595 allocate(wk(2,neqns), stat=ierr)
1596 if(ierr .ne. 0) then
1597 call errtrp('stop due to allocation error.')
1598 end if
1599 wk = 0
1600 do i=1,neqns_a
1601 wk(1,i) = b(1,iperm(i))
1602 wk(2,i) = b(2,iperm(i))
1603 end do
1604
1605 ! STEP22: forward substitution for A
1606 call elapout('nusol2_child: begin forward substitution for A') !elap
1607 do i=1, neqns_a
1608 ks=xlnzr(i)
1609 ke=xlnzr(i+1)-1
1610 if(ke.ge.ks) then
1611 call s2pdot(wk(:,i),wk,zln,colno,ks,ke)
1612 end if
1613 end do
1614
1615 ! STEP23: forward substitution for C and send it (yi) to parent
1616 call elapout('nusol2_child: begin forward substitution for C') !elap
1617 allocate(wk_d(2,nstop:neqns), stat=ierr)
1618 if(ierr .ne. 0) then
1619 call errtrp('stop due to allocation error.')
1620 end if
1621 wk_d=0
1622
1623 do i=nstop,neqns
1624 ks=xlnzr(i)
1625 ke=xlnzr(i+1)-1
1626 if(ke.ge.ks) then
1627 call s2pdot(wk_d(:,i),wk,zln,colno,ks,ke)
1628 end if
1629 end do
1630
1631 call elapout('nusol2_child: wait to send wk_d') !elap
1632 imp = m_pds_procinfo%imp
1633 call mpi_send(wk_d, 2*neqns_c, mpi_real8, imp, 1,mpi_comm_world, ierr)
1634
1635 ! STEP24: divide with diagonal matrix
1636 do i=1,neqns_a
1637 wk(2,i)=wk(2,i)-wk(1,i)*diag(2,i)
1638 wk(1,i)=wk(1,i)*diag(1,i)
1639 wk(2,i)=wk(2,i)*diag(3,i)
1640 wk(1,i)=wk(1,i)-wk(2,i)*diag(2,i)
1641 end do
1642
1643 ! STEP25: receive xd from parent and do backword substitution
1644 call elapout('nusol2_child: wait until receive wk_d') !elap
1645 call mpi_bcast(wk_d, 2*neqns_c, mpi_real8, imp, mpi_comm_world, ierr)
1646 call elapout('nusol2_child: end receive wk_d') !elap
1647 call elapout('nusol2_child: begin backword substitution') !elap
1648
1649 wk(:,nstop:neqns)=wk_d(:,nstop:neqns)
1650 do i=neqns,1,-1
1651 ks=xlnzr(i)
1652 ke=xlnzr(i+1)-1
1653 if(ke.ge.ks) then
1654 do k=ks,ke
1655 j=colno(k)
1656 wk(1,j)=wk(1,j)-wk(1,i)*zln(1,k)-wk(2,i)*zln(2,k)
1657 wk(2,j)=wk(2,j)-wk(1,i)*zln(3,k)-wk(2,i)*zln(4,k)
1658 end do
1659 end if
1660 end do
1661 call elapout('nusol2_child: end backword substitution') !elap
1662
1663 ! permutaion
1664 do i=1,neqns_a
1665 b(1,iperm(i))=wk(1,i)
1666 b(2,iperm(i))=wk(2,i)
1667 end do
1668
1669 call elapout('nusol2_child: end') !elap
1670 return
1671
1672 end subroutine nusol2_child
1673
1674 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1675
1676 subroutine nusol3_child(xlnzr,colno,zln,diag,iperm,b,neqns, nstop)
1677
1678 ! perform forward substitution for sparse matrix,
1679 ! send bd to parent and receive xd from parent,
1680 ! backword substitution using xd ad send final result x to parent.
1681
1682 use m_elap
1683
1684 implicit none
1685
1686 integer(kind=kint), intent(in) :: xlnzr(:), colno(:), iperm(:)
1687 real(kind=kreal), intent(in) :: zln(:,:), diag(:,:) !zln(9,*),diag(6,*)
1688 real(kind=kreal), intent(inout) :: b(:,:) ! b(3,*)
1689
1690 real(kind=kreal), allocatable :: wk(:,:), wk_d(:,:)
1691 integer(kind=kint) :: neqns_c, neqns_a, nstop, neqns, ks, ke
1692 integer(kind=kint) :: i, j, k, l, imp, ierr
1693
1694 include 'mpif.h'
1695
1696 ! now nstop is begining point of C
1697
1698 neqns_a = nstop - 1
1699 neqns_c = neqns - nstop + 1
1700
1701 call elapout('nusol3_child: entered') !elap
1702
1703 allocate(wk(3,neqns), stat=ierr)
1704 if(ierr .ne. 0) then
1705 call errtrp('stop due to allocation error.')
1706 end if
1707 wk = 0
1708 do 10 i=1,neqns_a
1709 wk(1,i)=b(1,iperm(i))
1710 wk(2,i)=b(2,iperm(i))
1711 wk(3,i)=b(3,iperm(i))
1712 10 continue
1713
1714
1715 ! STEP22: forward substitution for A
1716 call elapout('nusol3_child: begin forward substitution for A') !elap
1717
1718 do 100 i=1, neqns_a
1719 ks=xlnzr(i)
1720 ke=xlnzr(i+1)-1
1721 if(ke.lt.ks) goto 110
1722 call s3pdot(wk(:,i),wk,zln,colno,ks,ke)
1723 110 continue
1724 100 continue
1725
1726
1727 ! STEP23: forward substitution for C and send it (yi) to parent
1728 call elapout('nusol3_child: begin forward substitution for C') !elap
1729 allocate(wk_d(3,nstop:neqns), stat=ierr)
1730 if(ierr .ne. 0) then
1731 call errtrp('stop due to allocation error.')
1732 end if
1733 wk_d=0
1734
1735 do 101 i=nstop,neqns
1736 ks=xlnzr(i)
1737 ke=xlnzr(i+1)-1
1738 if(ke.lt.ks) goto 111
1739 call s3pdot(wk_d(:,i),wk,zln,colno,ks,ke)
1740 111 continue
1741 101 continue
1742
1743 call elapout('nusol3_child: wait to send wk_d') !elap
1744 imp = m_pds_procinfo%imp
1745 call mpi_send(wk_d, 3*neqns_c, mpi_real8, imp, 1,mpi_comm_world, ierr)
1746
1747 ! STEP24: divide with diagonal matrix
1748 call elapout('nusol3_child: divide with diagonal matrix') !elap
1749 do 120 i=1,neqns_a
1750 wk(2,i)=wk(2,i)-wk(1,i)*diag(2,i)
1751 wk(3,i)=wk(3,i)-wk(1,i)*diag(4,i)-wk(2,i)*diag(5,i)
1752 wk(1,i)=wk(1,i)*diag(1,i)
1753 wk(2,i)=wk(2,i)*diag(3,i)
1754 wk(3,i)=wk(3,i)*diag(6,i)
1755 wk(2,i)=wk(2,i)-wk(3,i)*diag(5,i)
1756 wk(1,i)=wk(1,i)-wk(2,i)*diag(2,i)-wk(3,i)*diag(4,i)
1757 120 continue
1758
1759 ! STEP25: receive xd from parent and do backword substitution
1760 call elapout('nusol3_child: wait until receive wk_d') !elap
1761 call mpi_bcast(wk_d, 3*neqns_c, mpi_real8, imp, mpi_comm_world, ierr)
1762 call elapout('nusol3_child: end receive wk_d') !elap
1763 call elapout('nusol3_child: begin backword substitution') !elap
1764
1765 wk(:,nstop:neqns)=wk_d(:,nstop:neqns)
1766
1767 do 200 i=neqns,1,-1
1768 ks=xlnzr(i)
1769 ke=xlnzr(i+1)-1
1770 if(ke.lt.ks) goto 200
1771
1772 do 210 k=ks,ke
1773 j=colno(k)
1774
1775 wk(1,j)=wk(1,j)-wk(1,i)*zln(1,k)-wk(2,i)*zln(2,k)-wk(3,i)*zln(3,k)
1776 wk(2,j)=wk(2,j)-wk(1,i)*zln(4,k)-wk(2,i)*zln(5,k)-wk(3,i)*zln(6,k)
1777 wk(3,j)=wk(3,j)-wk(1,i)*zln(7,k)-wk(2,i)*zln(8,k)-wk(3,i)*zln(9,k)
1778 210 continue
1779 200 continue
1780 call elapout('nusol3_child: end backword substitution') !elap
1781
1782
1783 ! permutaion
1784 do 300 i=1,neqns_a
1785 b(1,iperm(i))=wk(1,i)
1786 b(2,iperm(i))=wk(2,i)
1787 b(3,iperm(i))=wk(3,i)
1788 300 continue
1789
1790 call elapout('nusol3_child: end') !elap
1791 return
1792 end subroutine nusol3_child
1793
1794 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1795 subroutine nusolx_child(xlnzr, colno, zln, diag, iperm, b, neqns, nstop, ndeg)
1796
1797 implicit none
1798
1799 integer(kind=kint), intent(in) :: xlnzr(:), colno(:), iperm(:)
1800 real(kind=kreal), intent(in) :: zln(:,:), diag(:,:)
1801 real(kind=kreal), intent(inout) :: b(:,:)
1802 integer(kind=kint), intent(in) :: neqns, nstop, ndeg
1803
1804 real(kind=kreal), allocatable :: wk(:,:), wk_d(:,:)
1805 integer(kind=kint) :: neqns_c, neqns_a, ks, ke, locd, loc1
1806 integer(kind=kint) :: i, j, k, l, m, n, imp, ierr
1807
1808 include 'mpif.h'
1809
1810 !now nstop is begining point of C
1811 neqns_a = nstop - 1
1812 neqns_c = neqns - nstop + 1
1813
1814 allocate(wk(ndeg,neqns), stat=ierr)
1815 if(ierr .ne. 0) then
1816 call errtrp('stop due to allocation error.')
1817 end if
1818 wk = 0
1819 do i=1,neqns_a
1820 wk(1,i)=b(1,iperm(i))
1821 wk(2,i)=b(2,iperm(i))
1822 wk(3,i)=b(3,iperm(i))
1823 end do
1824
1825 ! STEP22: forward substitution for A
1826 do i=1,neqns_a
1827 ks=xlnzr(i)
1828 ke=xlnzr(i+1)-1
1829 if(ke.ge.ks) then
1830 call sxpdot(ndeg,wk(1,i),wk,zln,colno,ks,ke)
1831 end if
1832 end do
1833
1834 ! STEP23: forward substitution for C and send it (yi) to parent
1835 allocate(wk_d(ndeg,nstop:neqns), stat=ierr)
1836 if(ierr .ne. 0) then
1837 call errtrp('stop due to allocation error.')
1838 end if
1839 wk_d=0
1840 do i=nstop,neqns
1841 ks=xlnzr(i)
1842 ke=xlnzr(i+1)-1
1843 if(ke.ge.ks) then
1844 call sxpdot(ndeg,wk_d(:,i),wk,zln,colno,ks,ke)
1845 end if
1846 end do
1847
1848 imp = m_pds_procinfo%imp
1849 call mpi_send(wk_d, ndeg*neqns_c, mpi_real8, imp, 1,mpi_comm_world, ierr)
1850
1851
1852 ! STEP24: divide with diagonal matrix
1853 do i=1,neqns_a
1854 locd=0
1855 do m=1,ndeg-1
1856 locd=locd+m
1857 loc1=locd+m
1858 do n=m+1,ndeg
1859 wk(n,i)=wk(n,i)-wk(m,i)*diag(loc1,i)
1860 loc1=loc1+n
1861 end do
1862 end do
1863
1864 locd=0
1865 do m=1,ndeg
1866 locd=locd+m
1867 wk(m,i)=wk(m,i)*diag(locd,i)
1868 end do
1869
1870 do n=ndeg,2,-1
1871 locd=locd-1
1872 do m=n-1,1,-1
1873 wk(m,i)=wk(m,i)-wk(n,i)*diag(locd,i)
1874 locd=locd-1
1875 end do
1876 end do
1877 end do
1878
1879 call mpi_bcast(wk_d, ndeg*neqns_c, mpi_real8, imp, mpi_comm_world, ierr)
1880 wk(:,nstop:neqns)=wk_d(:,nstop:neqns)
1881
1882 ! back ward
1883 do i=neqns,1,-1
1884 ks=xlnzr(i)
1885 ke=xlnzr(i+1)-1
1886 if(ke.ge.ks) then
1887 do k=ks,ke
1888 j=colno(k)
1889 do m=1,ndeg
1890 do n=1,ndeg
1891 wk(m,j)=wk(m,j)-wk(n,i)*zln(n+(m-1)*ndeg,k)
1892 end do
1893 end do
1894 end do
1895 end if
1896 end do
1897
1898 ! permutaion
1899 do l=1,ndeg
1900 do i=1,neqns_a
1901 b(l,iperm(i))=wk(l,i)
1902 end do
1903 end do
1904
1905 return
1906 end subroutine nusolx_child
1907 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1908
1909 subroutine nufct0_parent(dsln, diag, neqns, ndeg)
1910 ! select LDU decomposer for dens matrix according to ndeg
1911
1912 implicit none
1913
1914 real(kind=kreal), intent(inout) :: dsln(:,:)
1915 real(kind=kreal), intent(inout) :: diag(:,:)
1916 integer(kind=kint), intent(in) :: neqns, ndeg
1917
1918 integer(kind=kint) :: ndegl
1919
1920 if (ndeg .eq. 1) then
1921 call sum3(neqns, dsln(1,:), diag(1,:))
1922 else if (ndeg .eq. 3) then
1923 call s3um3(neqns, dsln, diag)
1924 else
1925 ndegl = (ndeg+1)*ndeg/2
1926 call sxum3(neqns, dsln, diag, ndeg, ndegl)
1927 end if
1928
1929 return
1930 end subroutine nufct0_parent
1931
1932 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1933
1934 subroutine nusol0_parent(dsln, diag, b, neqns, ndeg)
1935 ! select solvers according to ndeg
1936
1937 implicit none
1938
1939 real(kind=kreal), intent(in) :: dsln(:,:)
1940 real(kind=kreal), intent(in) :: diag(:,:)
1941 real(kind=kreal), intent(inout) :: b(:,:)
1942
1943 integer(kind=kint), intent(in) :: neqns, ndeg
1944
1945 if (ndeg .eq. 1) then
1946 call nusol1_parent(dsln(1,:), diag(1,:), b(1,:), neqns)
1947 else if (ndeg .eq. 3) then
1948 call nusol3_parent(dsln, diag, b, neqns)
1949 else
1950 call nusolx_parent(dsln, diag, b, neqns, ndeg)
1951 end if
1952
1953 return
1954 end subroutine nusol0_parent
1955
1956 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1957
1958 subroutine nusol1_parent(dsln, diag, b, neqns)
1959 ! solve Ax=b for dens matrix with ndeg=1
1960 ! require dsln, diag is already LDU decomposed.
1961 ! currently not tested 20071121
1962
1963 implicit none
1964
1965 real(kind=kreal), intent(in) :: dsln(:) !((neqns+1)*neqns/2)
1966 real(kind=kreal), intent(in) :: diag(:) !(neqns)
1967 real(kind=kreal), intent(inout) :: b(:) !(3,neqns)
1968 integer(kind=kint), intent(in) :: neqns
1969
1970 integer(kind=kint) :: i,j,k,l,loc
1971
1972 ! forward substitution
1973 do i=2,neqns
1974 k=(i-1)*(i-2)/2 + 1 ! first element of i'th row.
1975 b(i)=b(i)-dot_product(b(1:i-1),dsln(k:k+i-2))
1976 end do
1977
1978 ! divide by D (because of diag is already inverted (1/Dii))
1979 b(:)=b(:)*diag(:)
1980
1981 ! Backword substitution.
1982 ! Substitute Zi into D and get Xd results.
1983 loc=(neqns-1)*neqns/2
1984 do i=neqns,1,-1
1985 do j=i-1,1,-1
1986 b(j)=b(j)-b(i)*dsln(loc)
1987 loc=loc-1
1988 end do
1989 end do
1990
1991 return
1992 end subroutine nusol1_parent
1993
1994 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1995
1996 subroutine nusol2_parent(dsln, diag, b, neqns)
1997 ! solve Ax=b for dens matrix with ndeg=3
1998 ! require dsln, diag is already LDU decomposed.
1999
2000 implicit none
2001
2002 real(kind=kreal), intent(in) :: dsln(:,:) !(4, (neqns+1)*newns/2)
2003 real(kind=kreal), intent(in) :: diag(:,:) !(3,neqns)
2004 real(kind=kreal), intent(inout) :: b(:,:) !(2,neqns)
2005 integer(kind=kint), intent(in) :: neqns
2006
2007 integer(kind=kint) :: i,j,k,l,loc
2008
2009 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2010 !
2011 ! STEP22 forward substitution
2012 !
2013 do i=2,neqns
2014 k=(i-1)*(i-2)/2 + 1 ! first element of i'th row.
2015 call d2sdot(b(:,i),b,dsln(:, k:k+i-2),i-1)
2016 end do
2017
2018 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2019 !
2020 ! STEP23 divide Yd by diagonal elment of D and get Zi=Yi/Di
2021 !
2022 do i=1,neqns
2023 b(2,i)=b(2,i)-b(1,i)*diag(2,i)
2024 b(1,i)=b(1,i)*diag(1,i)
2025 b(2,i)=b(2,i)*diag(3,i)
2026 b(1,i)=b(1,i)-b(2,i)*diag(2,i)
2027 end do
2028
2029 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2030 !
2031 ! STEP24 Backword substitution.
2032 ! Substitute Zi into D and get Xd results.
2033 !
2034 loc=(neqns-1)*neqns/2
2035 do i=neqns,1,-1
2036 do j=i-1,1,-1
2037 b(1,j)=b(1,j)-b(1,i)*dsln(1,loc)-b(2,i)*dsln(2,loc)
2038 b(2,j)=b(2,j)-b(1,i)*dsln(3,loc)-b(2,i)*dsln(4,loc)
2039 loc=loc-1
2040 end do
2041 end do
2042
2043 return
2044
2045 end subroutine nusol2_parent
2046
2047 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2048
2049 subroutine nusol3_parent(dsln, diag, b, neqns)
2050 ! solve Ax=b for dens matrix with ndeg=3
2051 ! require dsln, diag is already LDU decomposed.
2052
2053 implicit none
2054
2055 real(kind=kreal), intent(in) :: dsln(:,:) !(9,(neqns+1)*neqns/2)
2056 real(kind=kreal), intent(in) :: diag(:,:) !(6,neqns)
2057 real(kind=kreal), intent(inout) :: b(:,:) !(3,neqns)
2058 integer(kind=kint), intent(in) :: neqns
2059
2060 integer(kind=kint) :: i,j,k,l,loc
2061
2062 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2063 !
2064 ! STEP22 forward substitution
2065 !
2066 do i=2,neqns
2067 k=(i-1)*(i-2)/2 + 1 ! first element of i'th row.
2068 call d3sdot(b(:,i),b,dsln(:, k:k+i-2),i-1)
2069 end do
2070
2071 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2072 !
2073 ! STEP23 divide Yd by diagonal elment of D and get Zi=Yi/Di
2074 !
2075 do i=1,neqns
2076 b(2,i)=b(2,i)-b(1,i)*diag(2,i)
2077 b(3,i)=b(3,i)-b(1,i)*diag(4,i)-b(2,i)*diag(5,i)
2078 b(1,i)=b(1,i)*diag(1,i)
2079 b(2,i)=b(2,i)*diag(3,i)
2080 b(3,i)=b(3,i)*diag(6,i)
2081 b(2,i)=b(2,i)-b(3,i)*diag(5,i)
2082 b(1,i)=b(1,i)-b(2,i)*diag(2,i)-b(3,i)*diag(4,i)
2083 end do
2084
2085 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2086 !
2087 ! STEP24 Backword substitution.
2088 ! Substitute Zi into D and get Xd results.
2089 !
2090 loc=(neqns-1)*neqns/2
2091 do i=neqns,1,-1
2092 do j=i-1,1,-1
2093 b(1,j)=b(1,j)-b(1,i)*dsln(1,loc)-b(2,i)*dsln(2,loc)-b(3,i)*dsln(3,loc)
2094 b(2,j)=b(2,j)-b(1,i)*dsln(4,loc)-b(2,i)*dsln(5,loc)-b(3,i)*dsln(6,loc)
2095 b(3,j)=b(3,j)-b(1,i)*dsln(7,loc)-b(2,i)*dsln(8,loc)-b(3,i)*dsln(9,loc)
2096 loc=loc-1
2097 end do
2098 end do
2099
2100 return
2101 end subroutine nusol3_parent
2102
2103 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2104
2105 subroutine nusolx_parent (dsln,diag,b,neqns,ndeg)
2106
2107 implicit none
2108
2109 real(kind=kreal), intent(in) :: diag(:,:), dsln(:,:)
2110 real(kind=kreal), intent(inout) :: b(:,:)
2111 integer(kind=kint), intent(in) :: neqns, ndeg
2112
2113 integer(kind=kint) :: i,j,k,l,m,n,loc, locd, loc1
2114
2115 ! forward
2116 do i=2,neqns
2117 k=(i-1)*(i-2)/2 + 1 ! first element of i'th row.
2118 call dxsdot(ndeg,b(:,i),b,dsln(:,k:k+i-2),i-1)
2119 end do
2120
2121 ! divide
2122 do i=1,neqns
2123 locd=0
2124 do m=1,ndeg-1
2125 locd=locd+m
2126 loc1=locd+m
2127 do n=m+1,ndeg
2128 b(n,i)=b(n,i)-b(m,i)*diag(loc1,i)
2129 loc1=loc1+n
2130 end do
2131 end do
2132
2133 locd=0
2134 do m=1,ndeg
2135 locd=locd+m
2136 b(m,i)=b(m,i)*diag(locd,i)
2137 end do
2138
2139 do n=ndeg,2,-1
2140 locd=locd-1
2141 do m=n-1,1,-1
2142 b(m,i)=b(m,i)-b(n,i)*diag(locd,i)
2143 locd=locd-1
2144 end do
2145 end do
2146 end do
2147 ! back ward
2148 loc=(neqns-1)*neqns/2
2149 do i=neqns,1,-1
2150 do j=i-1,1,-1
2151 do m=1,ndeg
2152 do n=1,ndeg
2153 b(m,j)=b(m,j)-b(n,i)*dsln((m-1)*ndeg+n,loc)
2154 end do
2155 end do
2156 loc=loc-1
2157 end do
2158 end do
2159 return
2160 end subroutine nusolx_parent
2161
2162 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2163
2164 subroutine s3um2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln)
2165
2166 implicit none
2167
2168 integer(kind=kint), intent(in) :: neqns, nstop
2169 integer(kind=kint), intent(in) :: xlnzr(:),colno(:)
2170 real(kind=kreal), intent(inout) :: zln(:,:),diag(:,:) !zln(9,*),diag(6,*),dsln(9,*) !
2171 integer(kind=kint), pointer :: spdslnidx(:)
2172 real(kind=kreal), pointer :: spdslnval(:,:)
2173 integer(kind=kint), intent(out) :: nspdsln
2174
2175 real(kind=kreal), allocatable :: temp(:,:)
2176 integer(kind=kint), allocatable :: indx(:)
2177 logical :: ftflag
2178 integer(kind=kint) :: i,j,k,l,m,n, ic,ks,ke,ii,jj,jc,j1,j2,loc,ispdsln, ierr
2179
2180 allocate(temp(neqns,9),indx(neqns), stat=ierr)
2181 if(ierr .ne. 0) then
2182 call errtrp('stop due to allocation error.')
2183 end if
2184
2185 nspdsln=0
2186 do ic=nstop,neqns
2187 ks=xlnzr(ic)
2188 ke=xlnzr(ic+1)-1
2189 do k=ks,ke
2190 jj=colno(k)
2191 indx(jj)=ic
2192 end do
2193 do jc=nstop,ic-1
2194 j1=xlnzr(jc)
2195 j2=xlnzr(jc+1)
2196 do jj=xlnzr(jc),xlnzr(jc+1)-1
2197 j=colno(jj)
2198 if(indx(j).eq.ic) then
2199 nspdsln=nspdsln+1
2200 exit
2201 endif
2202 end do
2203 end do
2204 end do
2205 allocate(spdslnidx(nspdsln), stat=ierr)
2206 if(ierr .ne. 0) then
2207 call errtrp('stop due to allocation error.')
2208 end if
2209 allocate(spdslnval(9,nspdsln), stat=ierr)
2210 if(ierr .ne. 0) then
2211 call errtrp('stop due to allocation error.')
2212 end if
2213
2214 loc=0
2215 ispdsln=0
2216 spdslnval=0
2217 ftflag = .true.
2218 do 100 ic=nstop,neqns
2219 ! do 105 m=1,9
2220 ! do 105 jj=1,nstop
2221 ! temp(jj,m)=0.0d0
2222 ! 105 continue
2223 ks=xlnzr(ic)
2224 ke=xlnzr(ic+1)-1
2225 do 110 k=ks,ke
2226 jj=colno(k)
2227 temp(jj,1)=zln(1,k)
2228 temp(jj,2)=zln(2,k)
2229 temp(jj,3)=zln(3,k)
2230 temp(jj,4)=zln(4,k)
2231 temp(jj,5)=zln(5,k)
2232 temp(jj,6)=zln(6,k)
2233 temp(jj,7)=zln(7,k)
2234 temp(jj,8)=zln(8,k)
2235 temp(jj,9)=zln(9,k)
2236 indx(jj)=ic
2237 110 continue
2238 do 111 k=ks,ke
2239 jj=colno(k)
2240 zln(4,k)=temp(jj,4)-temp(jj,1)*diag(2,jj)
2241 zln(7,k)=temp(jj,7)-temp(jj,1)*diag(4,jj)-zln(4,k)*diag(5,jj)
2242 zln(1,k)=temp(jj,1)*diag(1,jj)
2243 zln(4,k)=zln(4,k)*diag(3,jj)
2244 zln(7,k)=zln(7,k)*diag(6,jj)
2245 zln(4,k)=zln(4,k)-zln(7,k)*diag(5,jj)
2246 zln(1,k)=zln(1,k)-zln(4,k)*diag(2,jj)-zln(7,k)*diag(4,jj)
2247 !
2248 zln(5,k)=temp(jj,5)-temp(jj,2)*diag(2,jj)
2249 zln(8,k)=temp(jj,8)-temp(jj,2)*diag(4,jj)-zln(5,k)*diag(5,jj)
2250 zln(2,k)=temp(jj,2)*diag(1,jj)
2251 zln(5,k)=zln(5,k)*diag(3,jj)
2252 zln(8,k)=zln(8,k)*diag(6,jj)
2253 zln(5,k)=zln(5,k)-zln(8,k)*diag(5,jj)
2254 zln(2,k)=zln(2,k)-zln(5,k)*diag(2,jj)-zln(8,k)*diag(4,jj)
2255 !
2256 zln(6,k)=temp(jj,6)-temp(jj,3)*diag(2,jj)
2257 zln(9,k)=temp(jj,9)-temp(jj,3)*diag(4,jj)-zln(6,k)*diag(5,jj)
2258 zln(3,k)=temp(jj,3)*diag(1,jj)
2259 zln(6,k)=zln(6,k)*diag(3,jj)
2260 zln(9,k)=zln(9,k)*diag(6,jj)
2261 zln(6,k)=zln(6,k)-zln(9,k)*diag(5,jj)
2262 zln(3,k)=zln(3,k)-zln(6,k)*diag(2,jj)-zln(9,k)*diag(4,jj)
2263 ! write(60,6000) k,(zln(llll,k),llll=1,9)
2264 !6000 format(i6,3d20.10/6x,3d20.10/6x,3d20.10)
2265 111 continue
2266 !
2267 do 112 k=ks,ke
2268 jj=colno(k)
2269 diag(1,ic)=diag(1,ic)-temp(jj,1)*zln(1,k)-temp(jj,4)*zln(4,k)-temp(jj,7)*zln(7,k)
2270 diag(2,ic)=diag(2,ic)-temp(jj,1)*zln(2,k)-temp(jj,4)*zln(5,k)-temp(jj,7)*zln(8,k)
2271 diag(3,ic)=diag(3,ic)-temp(jj,2)*zln(2,k)-temp(jj,5)*zln(5,k)-temp(jj,8)*zln(8,k)
2272 diag(4,ic)=diag(4,ic)-temp(jj,1)*zln(3,k)-temp(jj,4)*zln(6,k)-temp(jj,7)*zln(9,k)
2273 diag(5,ic)=diag(5,ic)-temp(jj,2)*zln(3,k)-temp(jj,5)*zln(6,k)-temp(jj,8)*zln(9,k)
2274 diag(6,ic)=diag(6,ic)-temp(jj,3)*zln(3,k)-temp(jj,6)*zln(6,k)-temp(jj,9)*zln(9,k)
2275 112 continue
2276 do 120 jc=nstop,ic-1
2277 loc=loc+1
2278 j1=xlnzr(jc)
2279 j2=xlnzr(jc+1)
2280 do 220 jj=xlnzr(jc),xlnzr(jc+1)-1
2281 j=colno(jj)
2282 if(indx(j).eq.ic) then
2283 if (ftflag) then
2284 ispdsln=ispdsln+1
2285 ftflag=.false.
2286 end if
2287 spdslnidx(ispdsln)=loc
2288 spdslnval(1,ispdsln)=spdslnval(1,ispdsln)-temp(j,1)*zln(1,jj)-temp(j,4)*zln(4,jj)-temp(j,7)*zln(7,jj)
2289 spdslnval(2,ispdsln)=spdslnval(2,ispdsln)-temp(j,2)*zln(1,jj)-temp(j,5)*zln(4,jj)-temp(j,8)*zln(7,jj)
2290 spdslnval(3,ispdsln)=spdslnval(3,ispdsln)-temp(j,3)*zln(1,jj)-temp(j,6)*zln(4,jj)-temp(j,9)*zln(7,jj)
2291 spdslnval(4,ispdsln)=spdslnval(4,ispdsln)-temp(j,1)*zln(2,jj)-temp(j,4)*zln(5,jj)-temp(j,7)*zln(8,jj)
2292 spdslnval(5,ispdsln)=spdslnval(5,ispdsln)-temp(j,2)*zln(2,jj)-temp(j,5)*zln(5,jj)-temp(j,8)*zln(8,jj)
2293 spdslnval(6,ispdsln)=spdslnval(6,ispdsln)-temp(j,3)*zln(2,jj)-temp(j,6)*zln(5,jj)-temp(j,9)*zln(8,jj)
2294 spdslnval(7,ispdsln)=spdslnval(7,ispdsln)-temp(j,1)*zln(3,jj)-temp(j,4)*zln(6,jj)-temp(j,7)*zln(9,jj)
2295 spdslnval(8,ispdsln)=spdslnval(8,ispdsln)-temp(j,2)*zln(3,jj)-temp(j,5)*zln(6,jj)-temp(j,8)*zln(9,jj)
2296 spdslnval(9,ispdsln)=spdslnval(9,ispdsln)-temp(j,3)*zln(3,jj)-temp(j,6)*zln(6,jj)-temp(j,9)*zln(9,jj)
2297 endif
2298 220 continue
2299 ftflag = .true.
2300 120 continue
2301 100 continue
2302 return
2303 end subroutine s3um2_child
2304
2305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2306
2307 subroutine zpivot(neqns,neqnsz,nttbr,jcol,irow,zpiv,ir)
2308
2309 implicit none
2310
2311 integer(kind=kint), intent(in) :: jcol(:),irow(:)
2312 integer(kind=kint), intent(out) :: zpiv(:)
2313 integer(kind=kint), intent(in) :: neqns,nttbr
2314 integer(kind=kint), intent(out) :: neqnsz,ir
2315
2316 integer(kind=kint) :: i,j,k,l
2317
2318 ir=0
2319 do 100 l=1,neqns
2320 zpiv(l)=1
2321 100 continue
2322
2323 do 200 l=1,nttbr
2324 i=irow(l)
2325 j=jcol(l)
2326 if(i.le.0.or.j.le.0) then
2327 ir=-1
2328 goto 1000
2329 elseif(i.gt.neqns.or.j.gt.neqns) then
2330 ir=1
2331 goto 1000
2332 endif
2333 if(i.eq.j) zpiv(i)=0
2334 200 continue
2335
2336 do 310 i=neqns,1,-1
2337 if(zpiv(i).eq.0) then
2338 neqnsz=i
2339 goto 320
2340 endif
2341 310 continue
2342 320 continue
2343 1000 continue
2344 if(ldbg) write(idbg,*) '# zpivot ########################'
2345 if(ldbg) write(idbg,60) (zpiv(i),i=1,neqns)
2346 60 format(20i3)
2347 return
2348 end subroutine zpivot
2349
2350 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2351
2352 subroutine stsmat(neqns,nttbr,irow,jcol,jcpt,jcolno)
2353
2354 implicit none
2355
2356 integer(kind=kint), intent(in) :: irow(:), jcol(:)
2357 integer(kind=kint), intent(out) :: jcpt(:), jcolno(:)
2358 integer(kind=kint), intent(in) :: neqns, nttbr
2359
2360 integer(kind=kint) :: i,j,k,l,loc,locr
2361
2362 do 10 i=1,2*nttbr
2363 jcpt(i)=0
2364 jcolno(i)=0
2365 10 continue
2366 do 20 i=1,neqns
2367 jcpt(i)=i+neqns
2368 jcolno(i+neqns)=i
2369 20 continue
2370
2371 k=2*neqns
2372 do 100 l=1,nttbr
2373 i=irow(l)
2374 j=jcol(l)
2375 if(i.eq.j) goto 100
2376 loc=jcpt(i)
2377 locr=i
2378 110 continue
2379 if(loc.eq.0) goto 120
2380 if(jcolno(loc).eq.j) then
2381 goto 100
2382 elseif(jcolno(loc).gt.j) then
2383 goto 130
2384 endif
2385 locr=loc
2386 loc=jcpt(loc)
2387 goto 110
2388 120 continue
2389 k=k+1
2390 jcpt(locr)=k
2391 jcolno(k)=j
2392 goto 150
2393 130 continue
2394 k=k+1
2395 jcpt(locr)=k
2396 jcpt(k)=loc
2397 jcolno(k)=j
2398 150 continue
2399 loc=jcpt(j)
2400 locr=j
2401 160 continue
2402 if(loc.eq.0) goto 170
2403 if(jcolno(loc).eq.i) then
2404 goto 100
2405 elseif(jcolno(loc).gt.i) then
2406 goto 180
2407 endif
2408 locr=loc
2409 loc=jcpt(loc)
2410 goto 160
2411 170 continue
2412 k=k+1
2413 jcpt(locr)=k
2414 jcolno(k)=i
2415 goto 100
2416 180 continue
2417 k=k+1
2418 jcpt(locr)=k
2419 jcpt(k)=loc
2420 jcolno(k)=i
2421 100 continue
2422 if(ldbg) then
2423 write(idbg,*) 'jcolno'
2424 write(idbg,60) (jcolno(i),i=1,k)
2425 write(idbg,*) 'jcpt'
2426 write(idbg,60) (jcpt(i),i=1,k)
2427 60 format(10i7)
2428 endif
2429 return
2430 end subroutine stsmat
2431
2432 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2433 subroutine stiaja(neqns,neqnsz,ia,ja,jcpt,jcolno)
2434
2435 implicit none
2436 !
2437 ! coded by t.arakawa
2438 !
2439 integer(kind=kint), intent(in) :: jcpt(:),jcolno(:)
2440 integer(kind=kint), intent(out) :: ia(:),ja(:)
2441 integer(kind=kint), intent(in) :: neqns, neqnsz
2442
2443 integer(kind=kint) :: i,j,k,l,ii,loc
2444 !
2445
2446 ia(1)=1
2447 l=0
2448 do 100 k=1,neqns
2449 loc=jcpt(k)
2450 110 continue
2451 if(loc.eq.0) goto 120
2452 ii=jcolno(loc)
2453 if(ii.eq.k.or.ii.gt.neqnsz) goto 130
2454 l=l+1
2455 ja(l)=ii
2456 130 continue
2457 loc=jcpt(loc)
2458 goto 110
2459 120 ia(k+1)=l+1
2460 100 continue
2461 if(ldbg) then
2462 write(idbg,*) 'stiaja(): ia '
2463 write(idbg,60) (ia(i),i=1,neqns+1)
2464 write(idbg,*) 'stiaja(): ja '
2465 write(idbg,60) (ja(i),i=1,ia(neqns+1))
2466 endif
2467 60 format(10i7)
2468 return
2469 end subroutine stiaja
2470
2471 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2472
2473 subroutine idntty(neqns,invp,iperm)
2474
2475 implicit none
2476
2477 integer(kind=kint), intent(out) :: invp(:),iperm(:)
2478 integer(kind=kint), intent(in) :: neqns
2479
2480 integer(kind=kint) :: i
2481
2482 do 100 i=1,neqns
2483 invp(i)=i
2484 iperm(i)=i
2485 100 continue
2486 return
2487 end subroutine idntty
2488
2489 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2490
2491 subroutine genqmd(neqns,xadj,adj0,perm,invp,deg,marker,rchset,nbrhd,qsize,qlink,nofsub,adjncy)
2492
2493 implicit none
2494
2495 integer(kind=kint), intent(in) :: adj0(:),xadj(:)
2496 integer(kind=kint), intent(out) :: rchset(:),nbrhd(:),adjncy(:),perm(:),invp(:),deg(:),marker(:),qsize(:),qlink(:)
2497 integer(kind=kint), intent(in) :: neqns
2498 integer(kind=kint), intent(out) :: nofsub
2499
2500 integer(kind=kint) :: inode,ip,irch,mindeg,nhdsze,node,np,num,nump1,nxnode,rchsze,search,thresh,ndeg
2501 integer(kind=kint) :: i,j,k,l
2502
2503 mindeg=neqns
2504 nofsub=0
2505 do 10 i=1,xadj(neqns+1)-1
2506 adjncy(i)=adj0(i)
2507 10 continue
2508 do 100 node=1,neqns
2509 perm(node)=node
2510 invp(node)=node
2511 marker(node)=0
2512 qsize(node)=1
2513 qlink(node)=0
2514 ndeg=xadj(node+1)-xadj(node)
2515 deg(node)=ndeg
2516 if(ndeg.lt.mindeg) mindeg=ndeg
2517 100 continue
2518
2519 num=0
2520 200 search=1
2521 thresh=mindeg
2522 mindeg=neqns
2523 300 nump1=num+1
2524 if(nump1.gt.search) search=nump1
2525 do 400 j=search,neqns
2526 node=perm(j)
2527 if(marker(node).lt.0) goto 400
2528 ndeg=deg(node)
2529 if(ndeg.le.thresh) goto 500
2530 if(ndeg.lt.mindeg) mindeg=ndeg
2531 400 continue
2532 goto 200
2533
2534 500 search=j
2535 nofsub=nofsub+deg(node)
2536 marker(node)=1
2537 call qmdrch(node,xadj,adjncy,deg,marker,rchsze,rchset,nhdsze,nbrhd)
2538 nxnode=node
2539 600 num=num+1
2540 np=invp(nxnode)
2541 ip=perm(num)
2542 perm(np)=ip
2543 invp(ip)=np
2544 perm(num)=nxnode
2545 invp(nxnode)=num
2546 deg(nxnode)=-1
2547 nxnode=qlink(nxnode)
2548 if(nxnode.gt.0) goto 600
2549 if(rchsze.le.0) goto 800
2550 !
2551 call qmdupd(xadj,adjncy,rchsze,rchset,deg,qsize,qlink,marker,rchset(rchsze+1:),nbrhd(nhdsze+1:))
2552 marker(node)=0
2553 do 700 irch=1,rchsze
2554 inode=rchset(irch)
2555 if(marker(inode).lt.0) goto 700
2556 marker(inode)=0
2557 ndeg=deg(inode)
2558 if(ndeg.lt.mindeg) mindeg=ndeg
2559 if(ndeg.gt.thresh) goto 700
2560 mindeg=thresh
2561 thresh=ndeg
2562 search=invp(inode)
2563 700 continue
2564 if(nhdsze.gt.0) call qmdot(node,xadj,adjncy,marker,rchsze,rchset,nbrhd)
2565 800 if(num.lt.neqns) goto 300
2566 return
2567 end subroutine genqmd
2568
2569 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2570
2571 subroutine genpaq(xadj,adjncy,invp,iperm,parent,neqns,ancstr)
2572
2573 implicit none
2574
2575 integer(kind=kint), intent(in) :: xadj(:),adjncy(:),invp(:),iperm(:)
2576 integer(kind=kint), intent(out) :: parent(:),ancstr(:)
2577 integer(kind=kint), intent(in) :: neqns
2578
2579 integer(kind=kint) :: i,j,k,l,ip,it
2580
2581 do 100 i=1,neqns
2582 parent(i)=0
2583 ancstr(i)=0
2584 ip=iperm(i)
2585 do 110 k=xadj(ip),xadj(ip+1)-1
2586 l=invp(adjncy(k))
2587 if(l.ge.i) goto 110
2588 112 continue
2589 if(ancstr(l).eq.0) goto 111
2590 if(ancstr(l).eq.i) goto 110
2591 it=ancstr(l)
2592 ancstr(l)=i
2593 l=it
2594 goto 112
2595 111 continue
2596 ancstr(l)=i
2597 parent(l)=i
2598 110 continue
2599 100 continue
2600 do 200 i=1,neqns
2601 if(parent(i).eq.0) parent(i)=neqns+1
2602 200 continue
2603 parent(neqns+1)=0
2604 if(ldbg) write(idbg,6010)
2605 if(ldbg) write(idbg,6000) (i,parent(i),i=1,neqns)
2606 6000 format(2i6)
2607 6010 format(' parent')
2608 return
2609 end subroutine genpaq
2610
2611 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2612
2613 subroutine genbtq(xadj,adjncy,invp,iperm,parent,btree,zpiv,izz,neqns)
2614
2615 implicit none
2616
2617 integer(kind=kint), intent(in) :: xadj(:),adjncy(:),parent(:),invp(:),iperm(:),zpiv(:)
2618 integer(kind=kint), intent(out) :: btree(:,:) ! btree is (2,:)
2619 integer(kind=kint), intent(in) :: neqns
2620 integer(kind=kint), intent(out) :: izz
2621
2622 integer(kind=kint) :: i,j,k,l,ip,ib,inext
2623
2624 do 10 i=1,neqns+1
2625 btree(1,i)=0
2626 btree(2,i)=0
2627 10 continue
2628 do 100 i=1,neqns+1
2629 ip=parent(i)
2630 if(ip.le.0) goto 100
2631 ib=btree(1,ip)
2632 if(ib.eq.0) then
2633 btree(1,ip)=i
2634 else
2635 101 continue
2636 inext=btree(2,ib)
2637 if(inext.eq.0) then
2638 btree(2,ib)=i
2639 else
2640 ib=inext
2641 goto 101
2642 endif
2643 endif
2644 100 continue
2645 !
2646 ! find zeropivot
2647 !
2648 do 200 i=1,neqns
2649 if(zpiv(i).ne.0) then
2650 if(btree(1,invp(i)).eq.0) then
2651 izz=i
2652 goto 210
2653 endif
2654 endif
2655 200 continue
2656 izz=0
2657 210 continue
2658 if(ldbg) write(idbg,6010)
2659 if(ldbg) write(idbg,6000) (i,btree(1,i),btree(2,i),i=1,neqns)
2660 if(ldbg) write(idbg,6020) izz
2661 ! if(idbg1.ge.2) write(10,6100) neqns
2662 ! if(idbg1.ge.2) write(10,6100) (btree(1,i),btree(2,i),i=1,neqns)
2663 6000 format(i6,'(',2i6,')')
2664 6010 format(' binary tree')
2665 6020 format(' the first zero pivot is ',i4)
2666 6100 format(2i8)
2667 return
2668 end subroutine genbtq
2669
2670 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2671
2672 subroutine rotate(xadj,adjncy,invp,iperm,parent,btree,izz,neqns,anc,adjt,irr)
2673
2674 implicit none
2675
2676 integer(kind=kint), intent(in) :: xadj(:),adjncy(:),parent(:),btree(:,:)
2677 integer(kind=kint), intent(out) :: anc(:),adjt(:),invp(:),iperm(:)
2678 integer(kind=kint), intent(in) :: neqns,izz
2679 integer(kind=kint), intent(out) :: irr
2680
2681 integer(kind=kint) :: i,j,k,l,izzz,nanc,loc,locc,ll,kk,iy
2682
2683 !----------------------------------------------------------------------
2684 ! irr return code irr=0 node izz is not a bottom node
2685 ! irr=1 is a bottom node then rotation is
2686 ! performed
2687 !
2688 !----------------------------------------------------------------------
2689 if(izz.eq.0) then
2690 irr=0
2691 return
2692 endif
2693 izzz=invp(izz)
2694 if(btree(1,izzz).ne.0) then
2695 irr=0
2696 ! return
2697 endif
2698 irr=1
2699 !
2700 ! ancestors of izzz
2701 !
2702 nanc=0
2703 loc=izzz
2704 100 continue
2705 nanc=nanc+1
2706 anc(nanc)=loc
2707 loc=parent(loc)
2708 if(loc.ne.0) goto 100
2709 !
2710 ! to find the eligible node from ancestors of izz
2711 !
2712 ! adjt = Adj(Tree(y))
2713 l=1
2714 200 continue
2715 do 210 i=1,neqns
2716 adjt(i)=0
2717 210 continue
2718 locc=anc(l)
2719 220 continue
2720 loc=locc
2721 locc=btree(1,loc)
2722 if(locc.ne.0) goto 220
2723 230 continue
2724 do 240 k=xadj(iperm(loc)),xadj(iperm(loc)+1)-1
2725 adjt(invp(adjncy(k)))=1
2726 240 continue
2727 if(loc.ge.anc(l)) goto 250
2728 locc=btree(2,loc)
2729 if(locc.ne.0) goto 220
2730 loc=parent(loc)
2731 goto 230
2732 250 continue
2733 do 260 ll=l+1,nanc
2734 if(adjt(anc(ll)).eq.0) then
2735 l=l+1
2736 goto 200
2737 endif
2738 260 continue
2739 if(l.eq.1) goto 500
2740
2741 !
2742 ! anc(l-1) is the eligible node
2743 !
2744 ! (1) number the node not in Ancestor(iy)
2745 iy=anc(l-1)
2746 do 300 i=1,neqns
2747 adjt(i)=0
2748 300 continue
2749 do 310 ll=l,nanc
2750 adjt(anc(ll))=1
2751 310 continue
2752 k=0
2753 do 320 ll=1,neqns
2754 if(adjt(ll).eq.0) then
2755 k=k+1
2756 invp(iperm(ll))=k
2757 endif
2758 320 continue
2759 ! (2) followed by nodes in Ancestor(iy)-Adj(T(iy))
2760 330 continue
2761 do 340 i=1,neqns
2762 adjt(i)=0
2763 340 continue
2764 locc=iy
2765 350 continue
2766 loc=locc
2767 locc=btree(1,loc)
2768 if(locc.ne.0) goto 350
2769 360 continue
2770 do 370 kk=xadj(iperm(loc)),xadj(iperm(loc)+1)-1
2771 adjt(invp(adjncy(kk)))=1
2772 370 continue
2773 if(loc.ge.iy) goto 380
2774 locc=btree(2,loc)
2775 if(locc.ne.0) goto 350
2776 loc=parent(loc)
2777 goto 360
2778 380 continue
2779 do 390 ll=l,nanc
2780 if(adjt(anc(ll)).eq.0) then
2781 k=k+1
2782 invp(iperm(anc(ll)))=k
2783 endif
2784 390 continue
2785 ! (3) and finally number the node in Adj(t(iy))
2786 do 400 ll=l,nanc
2787 if(adjt(anc(ll)).ne.0) then
2788 k=k+1
2789 invp(iperm(anc(ll)))=k
2790 endif
2791 400 continue
2792 goto 600
2793 !
2794 ! izz can be numbered last
2795 !
2796 500 continue
2797 k=0
2798 do 510 i=1,neqns
2799 if(i.eq.izzz) goto 510
2800 k=k+1
2801 invp(iperm(i))=k
2802 510 continue
2803 invp(iperm(izzz))=neqns
2804 !
2805 ! set iperm
2806 !
2807 600 continue
2808 do 610 i=1,neqns
2809 iperm(invp(i))=i
2810 610 continue
2811 if(ldbg) write(idbg,6000) (invp(i),i=1,neqns)
2812 6000 format(10i6)
2813 return
2814 end subroutine rotate
2815
2816 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2817
2818 subroutine bringu(zpiv,iperm,invp,parent,izz,neqns,irr)
2819
2820 implicit none
2821
2822 integer(kind=kint), intent(in) :: zpiv(:),parent(:)
2823 integer(kind=kint), intent(out) :: iperm(:),invp(:)
2824 integer(kind=kint), intent(in) :: neqns,izz
2825 integer(kind=kint), intent(out) :: irr
2826
2827 integer(kind=kint) :: i,j,k,l,ib0,ib,ibp,izzp
2828
2829 !----------------------------------------------------------------------
2830 !
2831 ! bringu brings up zero pivots from bottom of the elimination tree
2832 ! to higher nodes
2833 !
2834 ! irr = 0 complete
2835 ! = 1 impossible
2836 !
2837 ! #coded by t.arakawa
2838 !
2839 !----------------------------------------------------------------------
2840
2841 irr=0
2842 ib0=invp(izz)
2843 ib=ib0
2844 100 continue
2845 if(ib.le.0) goto 1000
2846 ibp=parent(ib)
2847 izzp=iperm(ibp)
2848 if(zpiv(izzp).eq.0) goto 110
2849 ib=ibp
2850 goto 100
2851 110 continue
2852 invp(izz)=ibp
2853 invp(izzp)=ib0
2854 iperm(ibp)=izz
2855 iperm(ib0)=izzp
2856 if(ldbg) then
2857 do 200 i=1,neqns
2858 if(invp(iperm(i)).ne.i) goto 210
2859 if(iperm(invp(i)).ne.i) goto 210
2860 200 continue
2861 goto 220
2862 210 continue
2863 write(20,*) 'permutation error'
2864 stop
2865 endif
2866 220 continue
2867 return
2868 1000 continue
2869 irr=1
2870 end subroutine bringu
2871
2872 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2873
2874 subroutine posord(parent,btree,invp,iperm,pordr,nch,neqns,iw,qarent,mch)
2875
2876 implicit none
2877
2878 integer(kind=kint), intent(in) :: btree(:,:),qarent(:)
2879 integer(kind=kint), intent(out) :: pordr(:),invp(:),iperm(:),nch(:),iw(:),parent(:),mch(0:neqns+1)
2880 integer(kind=kint), intent(in) :: neqns
2881
2882 integer(kind=kint) :: i,j,k,l,locc,loc,locp,invpos,ipinv,ii
2883
2884 do 5 i=1,neqns
2885 mch(i)=0
2886 pordr(i)=0
2887 5 continue
2888 l=1
2889 locc=neqns+1
2890 10 continue
2891 loc=locc
2892 locc=btree(1,loc)
2893 if(locc.ne.0) goto 10
2894 locp=qarent(loc)
2895 mch(locp)=mch(locp)+1
2896 20 continue
2897 pordr(loc)=l
2898 if(l.ge.neqns) goto 1000
2899 l=l+1
2900 locc=btree(2,loc)
2901 if(locc.ne.0) goto 10
2902 loc=qarent(loc)
2903 locp=qarent(loc)
2904 mch(locp)=mch(locp)+mch(loc)+1
2905 goto 20
2906 1000 continue
2907 do 100 i=1,neqns
2908 ipinv=pordr(invp(i))
2909 invp(i)=ipinv
2910 iperm(ipinv)=i
2911 iw(pordr(i))=i
2912 100 continue
2913 do 110 i=1,neqns
2914 invpos=iw(i)
2915 nch(i)=mch(invpos)
2916 ii=qarent(invpos)
2917 if(ii.gt.0.and.ii.le.neqns) then
2918 parent(i)=pordr(ii)
2919 else
2920 parent(i)=qarent(invpos)
2921 endif
2922 110 continue
2923 if(ldbg) write(idbg,6020)
2924 if(ldbg) write(idbg,6000) (pordr(i),i=1,neqns)
2925 if(ldbg) write(idbg,6030)
2926 if(ldbg) write(idbg,6050)
2927 if(ldbg) write(idbg,6000) (parent(i),i=1,neqns)
2928 if(ldbg) write(idbg,6000) (invp(i),i=1,neqns)
2929 if(ldbg) write(idbg,6040)
2930 if(ldbg) write(idbg,6000) (iperm(i),i=1,neqns)
2931 if(ldbg) write(idbg,6010)
2932 if(ldbg) write(idbg,6000) (nch(i),i=1,neqns)
2933 6000 format(10i6)
2934 6010 format(' nch')
2935 6020 format(' post order')
2936 6030 format(/' invp ')
2937 6040 format(/' iperm ')
2938 6050 format(/' parent')
2939 return
2940 end subroutine posord
2941
2942 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2943
2944 subroutine gnleaf(xadj,adjncy,invp,iperm,pordr,nch,adjncp,xleaf,leaf,neqns,lnleaf)
2945
2946 implicit none
2947
2948 integer(kind=kint), intent(in) :: xadj(:),adjncy(:),pordr(:),nch(:),invp(:),iperm(:)
2949 integer(kind=kint), intent(out) :: xleaf(:),leaf(:),adjncp(:)
2950 integer(kind=kint), intent(in) :: neqns
2951
2952 integer(kind=kint) i,j,k,l,m,n,ik,istart,ip,iq,lnleaf,lc1,lc
2953
2954 l=1
2955 ik=0
2956 istart=0
2957 do 100 i=1,neqns
2958 xleaf(i)=l
2959 ip=iperm(i)
2960 do 105 k=xadj(ip),xadj(ip+1)-1
2961 iq=invp(adjncy(k))
2962 if(iq.lt.i) then
2963 ik=ik+1
2964 adjncp(ik)=iq
2965 endif
2966 105 continue
2967 m=ik-istart
2968 if(m.eq.0) goto 131
2969 call qqsort(adjncp(istart+1:),m)
2970 lc1=adjncp(istart+1)
2971 if(lc1.ge.i) goto 100
2972 leaf(l)=lc1
2973 l=l+1
2974 do 130 k=istart+2,ik
2975 lc=adjncp(k)
2976 ! if(lc.ge.i) goto 125
2977 if(lc1.lt.lc-nch(lc)) then
2978 leaf(l)=lc
2979 l=l+1
2980 endif
2981 125 continue
2982 lc1=lc
2983 130 continue
2984 ik=1
2985 istart=ik
2986 131 continue
2987 100 continue
2988 xleaf(neqns+1)=l
2989 lnleaf=l-1
2990 if(ldbg) write(idbg,6020)
2991 if(ldbg) write(idbg,6000) (xleaf(i),i=1,neqns+1)
2992 if(ldbg) write(idbg,6010) lnleaf
2993 if(ldbg) write(idbg,6000) (leaf(i),i=1,lnleaf)
2994 return
2995 6000 format(10i6)
2996 6010 format(' leaf (len = ',i6,')')
2997 6020 format(' xleaf')
2998 end subroutine gnleaf
2999
3000 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3001
3002 subroutine countclno(parent,xleaf,leaf,neqns,nstop,lncol,ir)
3003
3004 implicit none
3005
3006 ! Count total number of non-zero elements
3007 ! which include fill-in.
3008 ! A and C region of given sparse matrix will consider.
3009 ! D region will not consider because of D is treat as
3010 ! dens matrix.
3011 !
3012 integer(kind=kint), intent(in) :: parent(:),xleaf(:),leaf(:)
3013 integer(kind=kint), intent(in) :: neqns, nstop
3014 integer(kind=kint), intent(out) :: lncol, ir
3015
3016 integer(kind=kint) :: i,j,k,l,nc,ks,ke,nxleaf
3017
3018 nc=0
3019 ir=0
3020 l=1
3021 do 100 i=1,neqns
3022 ks=xleaf(i)
3023 ke=xleaf(i+1)-1
3024 if(ke.lt.ks) goto 100
3025 nxleaf=leaf(ks)
3026 do 110 k=ks,ke-1
3027 j=nxleaf
3028 nxleaf=leaf(k+1)
3029 105 continue
3030 if(j.ge.nxleaf) goto 110
3031 if(j.ge.nstop) then
3032 goto 100
3033 endif
3034 l=l+1
3035 j=parent(j)
3036 goto 105
3037 110 continue
3038 j=leaf(ke)
3039 115 continue
3040 if(j.ge.nstop) goto 100
3041 if(j.ge.i.or.j.eq.0) goto 100
3042 l=l+1
3043 j=parent(j)
3044 goto 115
3045 100 continue
3046 lncol=l-1
3047 return
3048 end subroutine countclno
3049
3050 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3051
3052 subroutine gnclno(parent,pordr,xleaf,leaf,xlnzr,colno,neqns,nstop,lncol,ir)
3053
3054 implicit none
3055
3056 integer(kind=kint), intent(in) :: parent(:),pordr(:),xleaf(:),leaf(:)
3057 integer(kind=kint), intent(out) :: colno(:),xlnzr(:)
3058 integer(kind=kint), intent(in) :: neqns, nstop
3059 integer(kind=kint), intent(out) :: lncol,ir
3060
3061 integer(kind=kint) :: i,j,k,l,nc,ks,ke,nxleaf
3062
3063 nc=0
3064 ir=0
3065 l=1
3066 do 100 i=1,neqns
3067 xlnzr(i)=l
3068 ks=xleaf(i)
3069 ke=xleaf(i+1)-1
3070 if(ke.lt.ks) goto 100
3071 nxleaf=leaf(ks)
3072 do 110 k=ks,ke-1
3073 j=nxleaf
3074 nxleaf=leaf(k+1)
3075 105 continue
3076 if(j.ge.nxleaf) goto 110
3077 if(j.ge.nstop) then
3078 goto 100
3079 endif
3080 colno(l)=j
3081 l=l+1
3082 j=parent(j)
3083 goto 105
3084 110 continue
3085 j=leaf(ke)
3086 115 continue
3087 if(j.ge.nstop) goto 100
3088 if(j.ge.i.or.j.eq.0) goto 100
3089 colno(l)=j
3090 l=l+1
3091 j=parent(j)
3092 goto 115
3093 100 continue
3094 xlnzr(neqns+1)=l
3095 lncol=l-1
3096 if(ldbg) write(idbg,6010)
3097 ! if(idbg1.ne.0) write(6,6000) (xlnzr(i),i=1,neqns+1)
3098 if(ldbg) write(idbg,6020) lncol
3099 if(ldbg) then
3100 do 200 k=1,neqns
3101 write(idbg,6100) k
3102 write(idbg,6000) (colno(i),i=xlnzr(k),xlnzr(k+1)-1)
3103 200 continue
3104 endif
3105 6000 format(10i4)
3106 6010 format(' xlnzr')
3107 6020 format(' colno (lncol =',i10,')')
3108 6100 format(/' row = ',i6)
3109 return
3110 end subroutine gnclno
3111
3112 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3113
3114 subroutine qmdrch(root,xadj,adjncy,deg,marker,rchsze,rchset,nhdsze,nbrhd)
3115
3116 implicit none
3117
3118 integer(kind=kint), intent(in) :: deg(:),xadj(:),adjncy(:)
3119 integer(kind=kint), intent(out) :: rchset(:),marker(:),nbrhd(:)
3120 integer(kind=kint), intent(in) :: root
3121 integer(kind=kint), intent(out) :: nhdsze,rchsze
3122
3123 integer(kind=kint) :: i,j,k,l, istrt, istop, jstrt, jstop, nabor, node
3124
3125 nhdsze=0
3126 rchsze=0
3127 istrt=xadj(root)
3128 istop=xadj(root+1)-1
3129 if(istop.lt.istrt) return
3130 do 600 i=istrt,istop
3131 nabor=adjncy(i)
3132 if(nabor.eq.0) return
3133 if(marker(nabor).ne.0) goto 600
3134 if(deg(nabor).lt.0) goto 200
3135 rchsze=rchsze+1
3136 rchset(rchsze)=nabor
3137 marker(nabor)=1
3138 goto 600
3139 200 marker(nabor)=-1
3140 nhdsze=nhdsze+1
3141 nbrhd(nhdsze)=nabor
3142 300 jstrt=xadj(nabor)
3143 jstop=xadj(nabor+1)-1
3144 do 500 j=jstrt,jstop
3145 node=adjncy(j)
3146 nabor=-node
3147 if(node) 300,600,400
3148 400 if(marker(node).ne.0) goto 500
3149 rchsze=rchsze+1
3150 rchset(rchsze)=node
3151 marker(node)=1
3152 500 continue
3153 600 continue
3154 return
3155 end subroutine qmdrch
3156
3157 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3158
3159 subroutine qmdupd(xadj,adjncy,nlist,list,deg,qsize,qlink,marker,rchset,nbrhd)
3160
3161 implicit none
3162
3163 integer(kind=kint), intent(in) :: adjncy(:),list(:),xadj(:)
3164 integer(kind=kint), intent(out) :: marker(:),nbrhd(:),rchset(:),deg(:),qsize(:),qlink(:)
3165 integer(kind=kint), intent(in) :: nlist
3166
3167 integer(kind=kint) :: i,j,k,l, deg0,deg1,il,inhd,inode,irch,jstrt,jstop,mark,nabor,nhdsze,node,rchsze
3168
3169 if(nlist.le.0) return
3170 deg0=0
3171 nhdsze=0
3172 do 200 il=1,nlist
3173 node=list(il)
3174 deg0=deg0+qsize(node)
3175 jstrt=xadj(node)
3176 jstop=xadj(node+1)-1
3177
3178 do 100 j=jstrt,jstop
3179 nabor=adjncy(j)
3180 if(marker(nabor).ne.0.or.deg(nabor).ge.0) goto 100
3181 marker(nabor)=-1
3182 nhdsze=nhdsze+1
3183 nbrhd(nhdsze)=nabor
3184 100 continue
3185 200 continue
3186
3187 if(nhdsze.gt.0) call qmdmrg(xadj,adjncy,deg,qsize,qlink,marker,deg0,nhdsze,nbrhd,rchset,nbrhd(nhdsze+1:))
3188 do 600 il=1,nlist
3189 node=list(il)
3190 mark=marker(node)
3191 if(mark.gt.1.or.mark.lt.0) goto 600
3192 call qmdrch(node,xadj,adjncy,deg,marker,rchsze,rchset,nhdsze,nbrhd)
3193 deg1=deg0
3194 if(rchsze.le.0) goto 400
3195 do 300 irch=1,rchsze
3196 inode=rchset(irch)
3197 deg1=deg1+qsize(inode)
3198 marker(inode)=0
3199 300 continue
3200 400 deg(node)=deg1-1
3201 if(nhdsze.le.0) goto 600
3202 do 500 inhd=1,nhdsze
3203 inode=nbrhd(inhd)
3204 marker(inode)=0
3205 500 continue
3206 600 continue
3207 return
3208 end subroutine qmdupd
3209
3210 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3211
3212 subroutine qmdot(root,xadj,adjncy,marker,rchsze,rchset,nbrhd)
3213
3214 implicit none
3215
3216 integer(kind=kint), intent(in) :: marker(:),rchset(:),nbrhd(:),xadj(:)
3217 integer(kind=kint), intent(out) :: adjncy(:)
3218 integer(kind=kint), intent(in) :: rchsze,root
3219
3220 integer(kind=kint) :: i,j,k,l,irch,inhd,node,jstrt,jstop,link,nabor
3221
3222 irch=0
3223 inhd=0
3224 node=root
3225 100 jstrt=xadj(node)
3226 jstop=xadj(node+1)-2
3227 if(jstop.lt.jstrt) goto 300
3228 do 200 j=jstrt,jstop
3229 irch=irch+1
3230 adjncy(j)=rchset(irch)
3231 if(irch.ge.rchsze) goto 400
3232 200 continue
3233 300 link=adjncy(jstop+1)
3234 node=-link
3235 if(link.lt.0) goto 100
3236 inhd=inhd+1
3237 node=nbrhd(inhd)
3238 adjncy(jstop+1)=-node
3239 goto 100
3240 400 adjncy(j+1)=0
3241 do 600 irch=1,rchsze
3242 node=rchset(irch)
3243 if(marker(node).lt.0) goto 600
3244 jstrt=xadj(node)
3245 jstop=xadj(node+1)-1
3246 do 500 j=jstrt,jstop
3247 nabor=adjncy(j)
3248 if(marker(nabor).ge.0) goto 500
3249 adjncy(j)=root
3250 goto 600
3251 500 continue
3252 600 continue
3253 return
3254 end subroutine qmdot
3255
3256 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3257
3258 subroutine qmdmrg(xadj,adjncy,deg,qsize,qlink,marker,deg0,nhdsze,nbrhd,rchset,ovrlp)
3259
3260 implicit none
3261
3262 integer(kind=kint), intent(in) :: adjncy(:),nbrhd(:),xadj(:)
3263 integer(kind=kint), intent(out) :: deg(:),marker(:),rchset(:),ovrlp(:),qsize(:),qlink(:)
3264 integer(kind=kint), intent(in) :: nhdsze
3265
3266 integer(kind=kint) :: i,j,k,l, deg0,deg1,head,inhd,iov,irch,jstrt,jstop,link,lnode,mark,mrgsze,nabor,node,novrlp,rchsze,root
3267
3268
3269 if(nhdsze.le.0) return
3270 do 100 inhd=1,nhdsze
3271 root=nbrhd(inhd)
3272 marker(root)=0
3273 100 continue
3274 do 1400 inhd=1,nhdsze
3275 root=nbrhd(inhd)
3276 marker(root)=-1
3277 rchsze=0
3278 novrlp=0
3279 deg1=0
3280 200 jstrt=xadj(root)
3281 jstop=xadj(root+1)-1
3282 do 600 j=jstrt,jstop
3283 nabor=adjncy(j)
3284 root=-nabor
3285 if(nabor) 200,700,300
3286 300 mark=marker(nabor)
3287 if(mark)600,400,500
3288 400 rchsze=rchsze+1
3289 rchset(rchsze)=nabor
3290 deg1=deg1+qsize(nabor)
3291 marker(nabor)=1
3292 goto 600
3293 500 if(mark.gt.1) goto 600
3294 novrlp=novrlp+1
3295 ovrlp(novrlp)=nabor
3296 marker(nabor)=2
3297 600 continue
3298 700 head=0
3299 mrgsze=0
3300 do 1100 iov=1,novrlp
3301 node=ovrlp(iov)
3302 jstrt=xadj(node)
3303 jstop=xadj(node+1)-1
3304 do 800 j=jstrt,jstop
3305 nabor=adjncy(j)
3306 if(marker(nabor).ne.0) goto 800
3307 marker(node)=1
3308 goto 1100
3309 800 continue
3310 mrgsze=mrgsze+qsize(node)
3311 marker(node)=-1
3312 lnode=node
3313 900 link=qlink(lnode)
3314 if(link.le.0) goto 1000
3315 lnode=link
3316 goto 900
3317 1000 qlink(lnode)=head
3318 head=node
3319 1100 continue
3320 if(head.le.0) goto 1200
3321 qsize(head)=mrgsze
3322 deg(head)=deg0+deg1-1
3323 marker(head)=2
3324 1200 root=nbrhd(inhd)
3325 marker(root)=0
3326 if(rchsze.le.0) goto 1400
3327 do 1300 irch=1,rchsze
3328 node=rchset(irch)
3329 marker(node)=0
3330 1300 continue
3331 1400 continue
3332 return
3333 end subroutine qmdmrg
3334
3335 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3336
3337 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)
3338 ! find fill-in position in C which placed under A and set it in xlnzr_c, colno_c
3339
3340 use m_crs_matrix
3341 implicit none
3342 !input
3343 integer(kind=kint), intent(in) :: xlnzr_a(:)
3344 integer(kind=kint), intent(in) :: colno_a(:)
3345 integer(kind=kint), intent(in) :: iperm_a(:)
3346 integer(kind=kint), intent(in) :: invp_a(:)
3347 integer(kind=kint), intent(in) :: ndeg
3348 integer(kind=kint), intent(in) :: nttbr_c
3349 integer(kind=kint), intent(in) :: irow_c(:)
3350 integer(kind=kint), intent(inout) :: jcol_c(:)
3351 integer(kind=kint), intent(in) :: ncol
3352 integer(kind=kint), intent(in) :: nrow
3353
3354 !output
3355 integer(kind=kint), pointer :: xlnzr_c(:)
3356 integer(kind=kint), pointer :: colno_c(:)
3357 integer(kind=kint), intent(out) :: lncol_c
3358
3359 ! internal
3360 integer(kind=kint) :: i,j,k,l,m,n
3361 integer(kind=kint) :: ks, ke, ipass, ierr
3362 logical, allocatable :: cnz(:)
3363 type(crs_matrix) :: crs_c
3364
3365 !permtate column in C for crs_c
3366 do i=1,nttbr_c
3367 jcol_c(i)=invp_a(jcol_c(i))
3368 end do
3369
3370 ! make Compact Column Storoge using symbolic information.
3371 call symbolicirjctocrs(ndeg, nttbr_c, irow_c, jcol_c, ncol, nrow, crs_c)
3372
3373 ! symbolic LDU factorization for C matrix
3374 allocate(cnz(ncol), stat=ierr)
3375 if(ierr .ne. 0) then
3376 call errtrp('stop due to allocation error.')
3377 end if
3378 do ipass = 1,2
3379 lncol_c = 0
3380 do k=1,nrow
3381 ! set cnz as non-zero pattern of C
3382 cnz = .false.
3383 ks = crs_c%ia(k)
3384 ke = crs_c%ia(k+1)-1
3385 if (ke .lt. ks) then
3386 if (ipass .eq. 2) then
3387 xlnzr_c(k+1)=lncol_c+1
3388 end if
3389 cycle ! in case of zero vector, no need to check dot product.
3390 end if
3391
3392 do i=ks,ke
3393 cnz(crs_c%ja(i)) = .true.
3394 end do
3395
3396 ! check for non-zero dot product and update cnz for each point of cnz
3397 do i=2,ncol
3398 ks = xlnzr_a(i)
3399 ke = xlnzr_a(i+1)-1
3400 if (ke .lt. ks) then ! in case of column of A is zero vector.
3401 cycle
3402 end if
3403 do j=ks,ke
3404 if (cnz(colno_a(j))) then
3405 cnz(i) = .true.
3406 exit
3407 end if
3408 end do
3409 end do
3410
3411 do i=1,ncol
3412 if (cnz(i)) then
3413 lncol_c = lncol_c + 1
3414 if (ipass .eq. 2) then
3415 colno_c(lncol_c) = i
3416 end if
3417 end if
3418 end do
3419 if (ipass .eq. 2) then
3420 xlnzr_c(k+1)=lncol_c + 1
3421 end if
3422 end do
3423
3424 if (ipass .eq. 1) then
3425 allocate(xlnzr_c(nrow+1),colno_c(lncol_c), stat=ierr)
3426 if(ierr .ne. 0) then
3427 call errtrp('stop due to allocation error.')
3428 end if
3429 xlnzr_c(1)=1
3430 end if
3431 end do
3432
3433 ! restore order of C column.
3434 do i=1,nttbr_c
3435 jcol_c(i)=iperm_a(jcol_c(i))
3436 end do
3437
3438 return
3439
3440 end subroutine ldudecomposec
3441
3442 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3443
3444 subroutine qqsort(iw,ik)
3445
3446 implicit none
3447
3448 integer(kind=kint), intent(out) :: iw(:)
3449 integer(kind=kint), intent(in) :: ik
3450
3451 integer(kind=kint) :: l,m,itemp
3452
3453 !----------------------------------------------------------------------
3454 ! sort in increasing order up to i
3455 !
3456 ! iw array
3457 ! ik number of input/output
3458 ! i deal with numbers less than this numberi
3459 !
3460 !----------------------------------------------------------------------
3461
3462 if(ik.le.1) return
3463 do 100 l=1,ik-1
3464 do 110 m=l+1,ik
3465 if(iw(l).lt.iw(m)) goto 110
3466 itemp=iw(l)
3467 iw(l)=iw(m)
3468 iw(m)=itemp
3469 110 continue
3470 100 continue
3471 200 continue
3472 return
3473 end subroutine qqsort
3474
3475
3476 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3477
3478 subroutine staij1(isw,i,j,aij,dsi,ir)
3479
3480 implicit none
3481
3482 !----------------------------------------------------------------------
3483 !
3484 ! this routine sets an non-zero entry of the matrix.
3485 ! (symmetric version)
3486 !
3487 ! (i)
3488 ! isw =0 set the value
3489 ! =1 add the value
3490 ! i row entry
3491 ! j column entry
3492 ! aij value
3493 !
3494 ! (o)
3495 ! iv communication array
3496 !
3497 ! #coded by t.arakawa
3498 !
3499 !----------------------------------------------------------------------
3500 !
3501 type(dsinfo) :: dsi
3502 real(kind=kreal), intent(out) :: aij(:) ! ndeg*ndeg
3503 integer(kind=kint), intent(in) :: isw, i, j
3504 integer(kind=kint), intent(out) :: ir
3505
3506 integer(kind=kint) :: ndeg, neqns, nstop, ndeg2, ndeg2l, ierr
3507 ndeg=dsi%ndeg
3508 neqns=dsi%neqns
3509 nstop=dsi%nstop
3510 ndeg2=ndeg*ndeg
3511 ndeg2l=ndeg*(ndeg+1)/2
3512
3513 ir=0
3514 ierr=0
3515
3516
3517 ! array allocation
3518 if(dsi%stage.ne.20) then
3519 if(dsi%stage.eq.30) write(ilog,*) 'Warning a matrix was build up but never solved.'
3520 !
3521 ! for diagonal
3522 !
3523 allocate(dsi%diag(ndeg2l,neqns), stat=ierr)
3524 if(ierr .ne. 0) then
3525 call errtrp('stop due to allocation error.')
3526 end if
3527 dsi%diag=0
3528 !
3529 ! for lower triangle
3530 !
3531 allocate(dsi%zln(ndeg2,dsi%lncol), stat=ierr)
3532 if(ierr .ne. 0) then
3533 call errtrp('stop due to allocation error.')
3534 end if
3535 dsi%zln=0
3536 !
3537 ! for dense window !TODO delete this and corresponding line in addr3()
3538 !
3539 ! allocate(dsi%dsln(ndeg2,dsi%lndsln))! because of there is no dense window
3540 ! dsi%dsln=0
3541
3542 dsi%stage=20
3543 endif
3544
3545 ! set value
3546 if(ndeg.le.2) then
3547 call addr0(isw,i,j,aij,dsi%invp,dsi%xlnzr,dsi%colno,dsi%diag,dsi%zln,dsi%dsln,nstop,dsi%ndeg,ir)
3548 elseif(ndeg.eq.3) then
3549 call addr3(isw,i,j,aij,dsi%invp,dsi%xlnzr,dsi%colno,dsi%diag,dsi%zln,dsi%dsln,nstop,ir)
3550 else
3551 call addrx(isw,i,j,aij,dsi%invp,dsi%xlnzr,dsi%colno,dsi%diag,dsi%zln,dsi%dsln,nstop,ndeg,ndeg2,ndeg2l,ir)
3552 endif
3553 1000 continue
3554 return
3555 end subroutine staij1
3556
3557 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3558 !
3559 ! After here, routines specilized for ndeg = 1
3560 !
3561 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3562
3563 ! LDU decompose of A (1..nstop-1) region
3564 subroutine sum(ic,xlnzr,colno,zln,diag,nch,par,neqns)
3565
3566 implicit none
3567
3568 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
3569 integer(kind=kint), intent(in) :: ic, neqns
3570 real(kind=kreal), intent(inout) :: zln(:),diag(:)
3571 integer(kind=kint), intent(out) :: nch(:)
3572
3573 real(kind=kreal) :: s, t, zz, piv
3574 integer(kind=kint) :: ks, ke, kk, k, jc, jj, j, ierr
3575 integer(kind=kint) :: isem
3576 real(kind=kreal),allocatable :: temp(:)
3577 integer(kind=kint),allocatable :: indx(:)
3578 allocate(temp(neqns),indx(neqns), stat=ierr)
3579 if(ierr .ne. 0) then
3580 call errtrp('stop due to allocation error.')
3581 end if
3582
3583 isem = 0
3584
3585 2 continue
3586 ks=xlnzr(ic)
3587 ke=xlnzr(ic+1)
3588 t=0.0d0
3589 ! do 100 i=1,ic
3590 ! temp(i)=0.0d0
3591 ! 100 continue
3592 do 200 k=ks,ke-1
3593 jc=colno(k)
3594 indx(jc)=ic
3595 s=0.0d0
3596 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
3597 j=colno(jj)
3598 if(indx(j).eq.ic) then
3599 s=s+temp(j)*zln(jj)
3600 endif
3601 310 continue
3602 ! j1=xlnzr(jc)
3603 ! jj=xlnzr(jc+1)-j1
3604 ! ss=ddoti(jj,zln(j1),colno(j1),temp)
3605 ! zz=zln(k)-ddoti(jj,zln(j1),colno(j1),temp)
3606 zz=zln(k)-s
3607 zln(k)=zz*diag(jc)
3608 temp(jc)=zz
3609 t=t+zz*zln(k)
3610 200 continue
3611 piv=diag(ic)-t
3612 if(dabs(piv).gt.rmin) then
3613 diag(ic)=1.0d0/piv
3614 endif
3615 1 continue
3616 if(isem.eq.1) then
3617 isem=0
3618 nch(ic)=-1
3619 kk=par(ic)
3620 nch(kk)=nch(kk)-1
3621 isem=1
3622 else
3623 goto 1
3624 endif
3625 return
3626 end subroutine sum
3627
3628 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3629
3630 ! LDU decompose of C (nstop..neqnsA+neqnsd) region
3631 subroutine sum1(ic,xlnzr,colno,zln,diag,par,neqns)
3632
3633 implicit none
3634
3635 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
3636 integer(kind=kint), intent(in) :: ic, neqns
3637 real(kind=kreal), intent(inout) :: zln(:),diag(:)
3638
3639 real(kind=kreal) :: s, t, zz
3640 integer(kind=kint) :: ks, ke, k, jc, j, jj, ierr
3641 real(kind=kreal),allocatable :: temp(:)
3642 integer(kind=kint),allocatable :: indx(:)
3643
3644 ierr=0
3645
3646 allocate(temp(neqns),indx(neqns), stat=ierr)
3647 if(ierr .ne. 0) then
3648 call errtrp('stop due to allocation error.')
3649 end if
3650
3651 ks=xlnzr(ic)
3652 ke=xlnzr(ic+1)
3653 t=0.0d0
3654 ! do 100 i=1,ic
3655 ! temp(i)=0.0d0
3656 ! 100 continue
3657 do 200 k=ks,ke-1
3658 jc=colno(k)
3659 indx(jc)=ic
3660 s=0.0d0
3661 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
3662 j=colno(jj)
3663 if(indx(j).eq.ic) then
3664 s=s+temp(j)*zln(jj)
3665 endif
3666 310 continue
3667 zz=zln(k)-s
3668 ! j1=xlnzr(jc)
3669 ! jj=xlnzr(jc+1)-j1
3670 ! zz=zln(k)-ddoti(jj,zln(j1),colno(j1),temp)
3671 zln(k)=zz
3672 temp(jc)=zz
3673 ! t=t+zz*zz*diag(jc)
3674 200 continue
3675 ! diag(ic)=diag(ic)-t
3676 return
3677 end subroutine sum1
3678
3679 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3680
3681 ! LDU decompose and Update D region.
3682 subroutine sum2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln)
3683
3684 implicit none
3685
3686 integer(kind=kint), intent(in) :: neqns, nstop
3687 integer(kind=kint), intent(in) :: xlnzr(:),colno(:)
3688 real(kind=kreal), intent(inout) :: zln(:),diag(:)
3689 integer(kind=kint), pointer :: spdslnidx(:)
3690 real(kind=kreal), pointer :: spdslnval(:,:)
3691 integer(kind=kint), intent(out) :: nspdsln
3692
3693 real(kind=kreal) :: s, t
3694 integer(kind=kint) :: ks, ke, kk, k, jc, jj, j, j1,j2
3695 integer(kind=kint) :: ic, i, loc, ierr
3696 integer(kind=kint) :: ispdsln
3697 logical :: ftflag
3698 real(kind=kreal),allocatable :: temp(:)
3699 integer(kind=kint),allocatable :: indx(:)
3700 ierr=0
3701 allocate(temp(neqns),indx(neqns), stat=ierr)
3702 if(ierr .ne. 0) then
3703 call errtrp('stop due to allocation error.')
3704 end if
3705
3706 nspdsln=0
3707 do ic=nstop,neqns
3708 ks=xlnzr(ic)
3709 ke=xlnzr(ic+1)-1
3710 do k=ks,ke
3711 jj=colno(k)
3712 indx(jj)=ic
3713 end do
3714 do jc=nstop,ic-1
3715 j1=xlnzr(jc)
3716 j2=xlnzr(jc+1)
3717 do jj=xlnzr(jc),xlnzr(jc+1)-1
3718 j=colno(jj)
3719 if(indx(j).eq.ic) then
3720 nspdsln=nspdsln+1
3721 exit
3722 endif
3723 end do
3724 end do
3725 end do
3726 allocate(spdslnidx(nspdsln),spdslnval(1,nspdsln), stat=ierr)
3727 if(ierr .ne. 0) then
3728 call errtrp('stop due to allocation error.')
3729 end if
3730
3731 loc=0
3732 ispdsln=0
3733 spdslnval=0
3734 ftflag = .true.
3735 do 100 ic=nstop,neqns
3736 do 105 i=1,nstop
3737 temp(i)=0.0d0
3738 105 continue
3739 ks=xlnzr(ic)
3740 ke=xlnzr(ic+1)-1
3741 do 110 k=ks,ke
3742 jj=colno(k)
3743 temp(jj)=zln(k)
3744 zln(k)=temp(jj)*diag(jj)
3745 indx(jj)=ic
3746 diag(ic)=diag(ic)-temp(jj)*zln(k)
3747 110 continue
3748 do 120 jc=nstop,ic-1
3749 s=0.0d0
3750 loc=loc+1
3751 do 220 jj=xlnzr(jc),xlnzr(jc+1)-1
3752 j=colno(jj)
3753 if(indx(j).eq.ic) then
3754 if (ftflag) then
3755 ispdsln=ispdsln+1
3756 ftflag=.false.
3757 end if
3758 s=s+temp(j)*zln(jj)
3759 endif
3760 220 continue
3761 spdslnidx(ispdsln)=loc
3762 spdslnval(1,ispdsln)=spdslnval(1,ispdsln)-s
3763 ftflag = .true.
3764 ! j1=xlnzr(jc)
3765 ! jj=xlnzr(jc+1)-j1
3766 ! dsln(loc)=dsln(loc)-ddoti(jj,zln(j1),colno(j1),temp)
3767 120 continue
3768 100 continue
3769 return
3770 end subroutine sum2_child
3771
3772 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3773
3774 subroutine sum3(n,dsln,diag)
3775
3776 implicit none
3777
3778 real(kind=kreal), intent(inout) :: dsln(:),diag(:)
3779 integer(kind=kint), intent(in) :: n
3780
3781 integer(kind=kint) :: i, j, loc, ierr
3782 real(kind=kreal),allocatable :: temp(:)
3783 integer(kind=kint),allocatable :: indx(:)
3784 allocate(temp(n),indx(n), stat=ierr)
3785 if(ierr .ne. 0) then
3786 call errtrp('stop due to allocation error.')
3787 end if
3788
3789 if(n.le.0) goto 1000
3790 indx(1)=0
3791 loc=1
3792 diag(1)=1.0d0/diag(1)
3793 do 100 i=2,n
3794 indx(i)=loc
3795 do 110 j=1,i-1
3796 dsln(loc)=dsln(loc)-dot_product(dsln(indx(i):indx(i)+j-2),dsln(indx(j):indx(j)+j-2))
3797 loc=loc+1
3798 110 continue
3799 temp(1:i-1)=dsln(indx(i):indx(i)+i-2)*diag(1:i-1)
3800 diag(i)=diag(i)-dot_product(temp(1:i-1),dsln(indx(i):indx(i)+i-2))
3801 dsln(indx(i):indx(i)+i-2)=temp(1:i-1)
3802 diag(i)=1.0d0/diag(i)
3803 100 continue
3804 1000 continue
3805 return
3806 end subroutine sum3
3807
3808 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3809
3810 real(kind=kreal) function spdot2(b,zln,colno,ks,ke)
3811
3812 implicit none
3813
3814 integer(kind=kint), intent(in) :: colno(:)
3815 integer(kind=kint), intent(in) :: ks,ke
3816 real(kind=kreal), intent(in) :: zln(:),b(:)
3817
3818 integer(kind=kint) :: j,jj
3819 real(kind=kreal) :: s
3820
3821 !----------------------------------------------------------------------
3822 !
3823 ! spdot1 performs inner product of sparse vectors
3824 !
3825 !
3826 ! #coded by t.arakawa
3827 !
3828 !----------------------------------------------------------------------
3829 !
3830 s=0.0d0
3831 do 100 jj=ks,ke
3832 j=colno(jj)
3833 s=s+zln(jj)*b(j)
3834 100 continue
3835 spdot2=s
3836 end function spdot2
3837
3838 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3839
3840 real(kind=kreal) function ddot(a,b,n)
3841
3842 implicit none
3843
3844 real(kind=kreal), intent(in) :: a(n),b(n)
3845 integer(kind=kint), intent(in) :: n
3846
3847 real(kind=kreal) :: s
3848 integer(kind=kint) :: i
3849
3850 s=0.0d0
3851 do 100 i=1,n
3852 s=s+a(i)*b(i)
3853 100 continue
3854 ddot=s
3855 return
3856 end function ddot
3857
3858 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3859
3860 subroutine addr0(isw,i,j,aij,invp,xlnzr,colno,diag,zln,dsln,nstop,ndeg,ir)
3861
3862 implicit none
3863
3864 integer(kind=kint), intent(in) :: isw ! 0: renew diag, dsln, zln other: add to diag, dsln, zln
3865 integer(kind=kint), intent(in) :: i,j,nstop, ndeg, invp(:),xlnzr(:),colno(:)
3866 real(kind=kreal), intent(inout) :: zln(:,:),diag(:,:),dsln(:,:),aij(:)
3867 integer(kind=kint), intent(out) :: ir
3868
3869 integer(kind=kint) :: ndeg2, ii, jj, itrans, k, i0, j0, l, ks, ke
3870 integer(kind=kint), parameter :: idbg=0
3871 ndeg2=ndeg*ndeg
3872
3873 ir=0
3874 ii=invp(i)
3875 jj=invp(j)
3876 if(idbg.ne.0) write(idbg,*) 'addr0',ii,jj,aij
3877 if(ii.eq.jj) then
3878 if(ndeg2.eq.1) then
3879 if(isw.eq.0) then
3880 diag(1,ii)=aij(1)
3881 else
3882 diag(1,ii)=diag(1,ii)+aij(1)
3883 endif
3884 elseif(ndeg2.eq.4) then
3885 if(isw.eq.0) then
3886 diag(1,ii)=aij(1)
3887 diag(2,ii)=aij(2)
3888 diag(3,ii)=aij(4)
3889 else
3890 diag(1,ii)=diag(1,ii)+aij(1)
3891 diag(2,ii)=diag(2,ii)+aij(2)
3892 diag(3,ii)=diag(3,ii)+aij(4)
3893 endif
3894 endif
3895 goto 1000
3896 endif
3897 itrans=0
3898 if(jj.gt.ii) then
3899 k=jj
3900 jj=ii
3901 ii=k
3902 itrans=1
3903 endif
3904 if(jj.ge.nstop) then
3905 i0=ii-nstop
3906 j0=jj-nstop+1
3907 k=i0*(i0-1)/2+j0
3908 if(ndeg2.eq.1) then
3909 dsln(1,k)=aij(1)
3910 goto 1000
3911 elseif(ndeg2.eq.4) then
3912 if(itrans.eq.0) then
3913 do 3 l=1,ndeg2
3914 dsln(l,k)=aij(l)
3915 3 continue
3916 goto 1000
3917 else
3918 dsln(1,k)=aij(1)
3919 dsln(2,k)=aij(3)
3920 dsln(3,k)=aij(2)
3921 dsln(4,k)=aij(4)
3922 goto 1000
3923 endif
3924 endif
3925 endif
3926 ks=xlnzr(ii)
3927 ke=xlnzr(ii+1)-1
3928 do 100 k=ks,ke
3929 if(colno(k).eq.jj) then
3930 if(isw.eq.0) then
3931 if(ndeg2.eq.1) then
3932 zln(1,k)=aij(1)
3933 elseif(ndeg2.eq.4) then
3934 if(itrans.eq.0) then
3935 do 4 l=1,ndeg2
3936 zln(l,k)=aij(l)
3937 4 continue
3938 else
3939 zln(1,k)=aij(1)
3940 zln(2,k)=aij(3)
3941 zln(3,k)=aij(2)
3942 zln(4,k)=aij(4)
3943 endif
3944 endif
3945 else
3946 if(ndeg2.eq.1) then
3947 zln(1,k)=zln(1,k)+aij(1)
3948 elseif(ndeg2.eq.4) then
3949 if(itrans.eq.0) then
3950 do 5 l=1,ndeg2
3951 zln(l,k)=zln(l,k)+aij(l)
3952 5 continue
3953 else
3954 zln(1,k)=zln(1,k)+aij(1)
3955 zln(2,k)=zln(2,k)+aij(3)
3956 zln(3,k)=zln(3,k)+aij(2)
3957 zln(4,k)=zln(4,k)+aij(4)
3958 endif
3959 endif
3960 endif
3961 goto 1000
3962 endif
3963 100 continue
3964 ir=20
3965 1000 continue
3966 return
3967 end subroutine addr0
3968
3969 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3970 !
3971 ! After here, routines specilized for ndeg = 3
3972 !
3973 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3974
3975 subroutine s3um(ic,xlnzr,colno,zln,diag,nch,par,neqns)
3976
3977 implicit none
3978
3979 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
3980 integer(kind=kint), intent(out) :: nch(:)
3981 real(kind=kreal), intent(out) :: zln(:,:), diag(:,:) ! zln(9,*),diag(6,*),temp(9,*),indx(*)
3982 integer(kind=kint), intent(in) :: ic,neqns
3983
3984 real(kind=kreal), allocatable :: temp(:,:)
3985 integer(kind=kint), allocatable :: indx(:)
3986 real(kind=kreal) :: zz(9),t(6)
3987 integer(kind=kint) :: i,j,k,l,ks,ke,kk,jc,jj,ir, ierr
3988
3989 allocate(temp(9,neqns),indx(neqns), stat=ierr)
3990 if(ierr .ne. 0) then
3991 call errtrp('stop due to allocation error.')
3992 end if
3993
3994 ks=xlnzr(ic)
3995 ke=xlnzr(ic+1)
3996 !$dir max_trips(6)
3997 do 100 l=1,6
3998 t(l)=0.0d0
3999 100 continue
4000 do 200 k=ks,ke-1
4001 jc=colno(k)
4002 indx(jc)=ic
4003 !$dir max_trips(9)
4004 do 210 l=1,9
4005 zz(l)=zln(l,k)
4006 210 continue
4007 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
4008 j=colno(jj)
4009 if(indx(j).eq.ic) then
4010 zz(1)=zz(1)-temp(1,j)*zln(1,jj)-temp(4,j)*zln(4,jj)-temp(7,j)*zln(7,jj)
4011 zz(2)=zz(2)-temp(2,j)*zln(1,jj)-temp(5,j)*zln(4,jj)-temp(8,j)*zln(7,jj)
4012 zz(3)=zz(3)-temp(3,j)*zln(1,jj)-temp(6,j)*zln(4,jj)-temp(9,j)*zln(7,jj)
4013 zz(4)=zz(4)-temp(1,j)*zln(2,jj)-temp(4,j)*zln(5,jj)-temp(7,j)*zln(8,jj)
4014 zz(5)=zz(5)-temp(2,j)*zln(2,jj)-temp(5,j)*zln(5,jj)-temp(8,j)*zln(8,jj)
4015 zz(6)=zz(6)-temp(3,j)*zln(2,jj)-temp(6,j)*zln(5,jj)-temp(9,j)*zln(8,jj)
4016 zz(7)=zz(7)-temp(1,j)*zln(3,jj)-temp(4,j)*zln(6,jj)-temp(7,j)*zln(9,jj)
4017 zz(8)=zz(8)-temp(2,j)*zln(3,jj)-temp(5,j)*zln(6,jj)-temp(8,j)*zln(9,jj)
4018 zz(9)=zz(9)-temp(3,j)*zln(3,jj)-temp(6,j)*zln(6,jj)-temp(9,j)*zln(9,jj)
4019 endif
4020 310 continue
4021
4022 call inv33(zln(:,k),zz,diag(:,jc))
4023
4024 !$dir max_trips(9)
4025 do 220 l=1,9
4026 temp(l,jc)=zz(l)
4027 220 continue
4028 t(1)=t(1)+zz(1)*zln(1,k)+zz(4)*zln(4,k)+zz(7)*zln(7,k)
4029 t(2)=t(2)+zz(1)*zln(2,k)+zz(4)*zln(5,k)+zz(7)*zln(8,k)
4030 t(3)=t(3)+zz(2)*zln(2,k)+zz(5)*zln(5,k)+zz(8)*zln(8,k)
4031 t(4)=t(4)+zz(1)*zln(3,k)+zz(4)*zln(6,k)+zz(7)*zln(9,k)
4032 t(5)=t(5)+zz(2)*zln(3,k)+zz(5)*zln(6,k)+zz(8)*zln(9,k)
4033 t(6)=t(6)+zz(3)*zln(3,k)+zz(6)*zln(6,k)+zz(9)*zln(9,k)
4034 200 continue
4035 !$dir max_trips(6)
4036 do 320 l=1,6
4037 diag(l,ic)=diag(l,ic)-t(l)
4038 320 continue
4039
4040 call inv3(diag(:,ic),ir)
4041 nch(ic)=-1
4042 kk=par(ic)
4043 nch(kk)=nch(kk)-1
4044
4045 return
4046 end subroutine s3um
4047
4048 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4049
4050 subroutine s3um1(ic,xlnzr,colno,zln,diag,nch,par,neqns)
4051
4052 implicit none
4053
4054 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),nch(:),par(:)
4055 real(kind=kreal), intent(in) :: diag(:,:)! zln(9,*),diag(6,*)
4056 real(kind=kreal), intent(out) :: zln(:,:)
4057 integer(kind=kint), intent(in) :: ic,neqns
4058
4059 integer(kind=kint) :: i,j,k,l,ks,ke,jc,jj,ierr
4060 real(kind=kreal) :: s(9),zz(9)
4061 real(kind=kreal),allocatable :: temp(:,:)
4062 integer(kind=kint),allocatable :: indx(:)
4063
4064 ierr=0
4065 allocate(temp(9,neqns),indx(neqns), stat=ierr)
4066 temp = 0.0d0
4067 indx = 0
4068
4069 ks=xlnzr(ic)
4070 ke=xlnzr(ic+1)
4071 !$dir max_trip(9)
4072 do 100 l=1,9
4073 s(l)=0.0d0
4074 100 continue
4075 do 200 k=ks,ke-1
4076 jc=colno(k)
4077 indx(jc)=ic
4078
4079 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
4080 j=colno(jj)
4081 if(indx(j).eq.ic) then
4082 s(1)=s(1)+temp(1,j)*zln(1,jj)+temp(4,j)*zln(4,jj)+temp(7,j)*zln(7,jj)
4083 s(2)=s(2)+temp(2,j)*zln(1,jj)+temp(5,j)*zln(4,jj)+temp(8,j)*zln(7,jj)
4084 s(3)=s(3)+temp(3,j)*zln(1,jj)+temp(6,j)*zln(4,jj)+temp(9,j)*zln(7,jj)
4085 s(4)=s(4)+temp(1,j)*zln(2,jj)+temp(4,j)*zln(5,jj)+temp(7,j)*zln(8,jj)
4086 s(5)=s(5)+temp(2,j)*zln(2,jj)+temp(5,j)*zln(5,jj)+temp(8,j)*zln(8,jj)
4087 s(6)=s(6)+temp(3,j)*zln(2,jj)+temp(6,j)*zln(5,jj)+temp(9,j)*zln(8,jj)
4088 s(7)=s(7)+temp(1,j)*zln(3,jj)+temp(4,j)*zln(6,jj)+temp(7,j)*zln(9,jj)
4089 s(8)=s(8)+temp(2,j)*zln(3,jj)+temp(5,j)*zln(6,jj)+temp(8,j)*zln(9,jj)
4090 s(9)=s(9)+temp(3,j)*zln(3,jj)+temp(6,j)*zln(6,jj)+temp(9,j)*zln(9,jj)
4091 endif
4092 310 continue
4093 !$dir max_trip(9)
4094 do 320 l=1,9
4095 temp(l,jc)=zln(l,k)-s(l)
4096 zln(l,k)=temp(l,jc)
4097 s(l)=0.0d0
4098 320 continue
4099
4100 200 continue
4101
4102 deallocate(temp,indx)
4103 end subroutine s3um1
4104
4105
4106 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4107
4108 subroutine s3um3(n,dsln,diag)
4109
4110 implicit none
4111
4112 real(kind=kreal), intent(out) :: dsln(:,:),diag(:,:)!dsln(9,*),diag(6,*)
4113 integer(kind=kint), intent(in) :: n
4114
4115 real(kind=kreal) :: t(9)
4116 integer(kind=kint) :: i,j,k,l,loc,ir, ierr
4117 real(kind=kreal), allocatable :: temp(:,:)
4118 integer(kind=kint), allocatable :: indx(:)
4119
4120 allocate(temp(9,n),indx(n), stat=ierr)
4121 if(ierr .ne. 0) then
4122 call errtrp('stop due to allocation error.')
4123 end if
4124
4125 if(n.le.0) goto 1000
4126 indx(1)=1
4127 loc=1
4128 call inv3(diag(:,1),ir)
4129 do 100 i=2,n
4130 indx(i)=loc
4131 do 110 j=1,i-1
4132 call d3dot(t,dsln(:,indx(i):indx(i)+j-2), dsln(:,indx(j):indx(j)+j-2),j-1)
4133 !$dir max_trips(9)
4134 ! do 111 l=1,9
4135 ! dsln(l,loc)=dsln(l,loc)-t(l)
4136 ! 111 continue
4137 dsln(:,loc)=dsln(:,loc)-t(:)
4138 loc=loc+1
4139 110 continue
4140 call v3prod(dsln(:,indx(i):indx(i)+i-2), diag,temp,i-1)
4141 call d3dotl(t,temp,dsln(:,indx(i):indx(i)+i-2),i-1)
4142 !$dir max_trips(6)
4143 ! do 112 l=1,6
4144 ! diag(l,i)=diag(l,i)-t(l)
4145 ! 112 continue
4146 diag(:,i)=diag(:,i)-t(1:6)
4147 dsln(:,indx(i):indx(i)+i-2)=temp(:,1:i-1)
4148 call inv3(diag(:,i),ir)
4149 100 continue
4150 1000 continue
4151 return
4152 end subroutine s3um3
4153
4154 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4155
4156 subroutine d3sdot(wi,a,b,n)
4157
4158 implicit none
4159
4160 real(kind=kreal), intent(in) :: a(:,:),b(:,:) !wi(3),a(3,*),b(9,*) !
4161 real(kind=kreal), intent(out) :: wi(:)
4162 integer(kind=kint), intent(in) :: n
4163
4164 integer(kind=kint) :: jj
4165 !
4166 !----------------------------------------------------------------------
4167 !
4168 ! spdot1 performs inner product of sparse vectors
4169 !
4170 !
4171 ! #coded by t.arakawa
4172 !
4173 !----------------------------------------------------------------------
4174 !
4175 do 100 jj=1,n
4176 wi(1)=wi(1)-a(1,jj)*b(1,jj)-a(2,jj)*b(4,jj)-a(3,jj)*b(7,jj)
4177 wi(2)=wi(2)-a(1,jj)*b(2,jj)-a(2,jj)*b(5,jj)-a(3,jj)*b(8,jj)
4178 wi(3)=wi(3)-a(1,jj)*b(3,jj)-a(2,jj)*b(6,jj)-a(3,jj)*b(9,jj)
4179 100 continue
4180 return
4181 end subroutine d3sdot
4182
4183 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4184
4185 subroutine s3pdot(bi,b,zln,colno,ks,ke)
4186
4187 implicit none
4188
4189 integer(kind=kint), intent(in) :: colno(:)
4190 real(kind=kreal), intent(in) :: zln(:,:),b(:,:) !zln(9,*),b(3,*),bi(3) !
4191 real(kind=kreal), intent(out) :: bi(:)
4192 integer(kind=kint), intent(in) :: ks,ke
4193
4194 integer(kind=kint) :: j,jj
4195 !
4196 !----------------------------------------------------------------------
4197 !
4198 ! spdot1 performs inner product of sparse vectors
4199 !
4200 !
4201 ! #coded by t.arakawa
4202 !
4203 !----------------------------------------------------------------------
4204 !
4205 do 100 jj=ks,ke
4206 j=colno(jj)
4207 bi(1)=bi(1)-zln(1,jj)*b(1,j)-zln(4,jj)*b(2,j)-zln(7,jj)*b(3,j)
4208 bi(2)=bi(2)-zln(2,jj)*b(1,j)-zln(5,jj)*b(2,j)-zln(8,jj)*b(3,j)
4209 bi(3)=bi(3)-zln(3,jj)*b(1,j)-zln(6,jj)*b(2,j)-zln(9,jj)*b(3,j)
4210 100 continue
4211 return
4212 end subroutine s3pdot
4213
4214 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4215
4216 subroutine inv33(zln,zz,diag)
4217
4218 implicit none
4219
4220 real(kind=kreal), intent(in) :: zz(9),diag(6)
4221 real(kind=kreal), intent(out) :: zln(9)
4222
4223 zln(4)=zz(4)-zz(1)*diag(2)
4224 zln(7)=zz(7)-zz(1)*diag(4)-zln(4)*diag(5)
4225 zln(1)=zz(1)*diag(1)
4226 zln(4)=zln(4)*diag(3)
4227 zln(7)=zln(7)*diag(6)
4228 zln(4)=zln(4)-zln(7)*diag(5)
4229 zln(1)=zln(1)-zln(4)*diag(2)-zln(7)*diag(4)
4230 !
4231 zln(5)=zz(5)-zz(2)*diag(2)
4232 zln(8)=zz(8)-zz(2)*diag(4)-zln(5)*diag(5)
4233 zln(2)=zz(2)*diag(1)
4234 zln(5)=zln(5)*diag(3)
4235 zln(8)=zln(8)*diag(6)
4236 zln(5)=zln(5)-zln(8)*diag(5)
4237 zln(2)=zln(2)-zln(5)*diag(2)-zln(8)*diag(4)
4238 !
4239 zln(6)=zz(6)-zz(3)*diag(2)
4240 zln(9)=zz(9)-zz(3)*diag(4)-zln(6)*diag(5)
4241 zln(3)=zz(3)*diag(1)
4242 zln(6)=zln(6)*diag(3)
4243 zln(9)=zln(9)*diag(6)
4244 zln(6)=zln(6)-zln(9)*diag(5)
4245 zln(3)=zln(3)-zln(6)*diag(2)-zln(9)*diag(4)
4246 return
4247 end subroutine inv33
4248
4249 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4250
4251 subroutine inv3(dsln,ir)
4252
4253 implicit none
4254
4255 real(kind=kreal) :: dsln(6),t(2)
4256 integer(kind=kint) :: ir
4257
4258 ir=0
4259 if(dabs(dsln(1)).lt.rmin) then
4260 goto 999
4261 endif
4262 dsln(1)=1.0d0/dsln(1)
4263 t(1)=dsln(2)*dsln(1)
4264 dsln(3)=dsln(3)-t(1)*dsln(2)
4265 dsln(2)=t(1)
4266 if(dabs(dsln(3)).lt.rmin) then
4267 goto 999
4268 endif
4269 dsln(3)=1.0d0/dsln(3)
4270 t(1)=dsln(4)*dsln(1)
4271 dsln(5)=dsln(5)-dsln(2)*dsln(4)
4272 t(2)=dsln(5)*dsln(3)
4273 dsln(6)=dsln(6)-t(1)*dsln(4)-t(2)*dsln(5)
4274 dsln(4)=t(1)
4275 dsln(5)=t(2)
4276 if(dabs(dsln(6)).lt.rmin) then
4277 goto 999
4278 endif
4279 dsln(6)=1.0d0/dsln(6)
4280 ! write(6,*) "dsln",dsln(1),dsln(3),dsln(6)
4281 return
4282 999 continue
4283 write(ilog,*) "singular"
4284 dsln(1)=1.0d0
4285 dsln(2)=0.0d0
4286 dsln(3)=1.0d0
4287 dsln(4)=0.0d0
4288 dsln(5)=0.0d0
4289 dsln(6)=1.0d0
4290 return
4291 end subroutine inv3
4292
4293 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4294
4295 subroutine d3dot(t,a,b,n)
4296 implicit none
4297
4298 real(kind=kreal), intent(in) :: a(:,:),b(:,:)
4299 real(kind=kreal), intent(out) :: t(:)
4300 integer(kind=kint),intent(in) :: n
4301
4302 integer(kind=kint) :: l,jj
4303 ! double precision a(9,n),b(9,n)
4304 ! double precision t(9)
4305 !
4306 !----------------------------------------------------------------------
4307 !
4308 ! spdot1 performs inner product of sparse vectors
4309 !
4310 ! it might be 'DENS' kitayama
4311 !
4312 !
4313 ! #coded by t.arakawa
4314 !
4315 !----------------------------------------------------------------------
4316 !
4317 !$dir max_trips(9)
4318 do 10 l=1,9
4319 t(l)=0.0d0
4320 10 continue
4321 do 100 jj=1,n
4322 t(1)=t(1)+a(1,jj)*b(1,jj)+a(4,jj)*b(4,jj)+a(7,jj)*b(7,jj)
4323 t(2)=t(2)+a(2,jj)*b(1,jj)+a(5,jj)*b(4,jj)+a(8,jj)*b(7,jj)
4324 t(3)=t(3)+a(3,jj)*b(1,jj)+a(6,jj)*b(4,jj)+a(9,jj)*b(7,jj)
4325 t(4)=t(4)+a(1,jj)*b(2,jj)+a(4,jj)*b(5,jj)+a(7,jj)*b(8,jj)
4326 t(5)=t(5)+a(2,jj)*b(2,jj)+a(5,jj)*b(5,jj)+a(8,jj)*b(8,jj)
4327 t(6)=t(6)+a(3,jj)*b(2,jj)+a(6,jj)*b(5,jj)+a(9,jj)*b(8,jj)
4328 t(7)=t(7)+a(1,jj)*b(3,jj)+a(4,jj)*b(6,jj)+a(7,jj)*b(9,jj)
4329 t(8)=t(8)+a(2,jj)*b(3,jj)+a(5,jj)*b(6,jj)+a(8,jj)*b(9,jj)
4330 t(9)=t(9)+a(3,jj)*b(3,jj)+a(6,jj)*b(6,jj)+a(9,jj)*b(9,jj)
4331 100 continue
4332 return
4333 end subroutine d3dot
4334
4335 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4336
4337 subroutine v3prod(zln,diag,zz,n)
4338
4339 implicit none
4340
4341 real(kind=kreal), intent(in) :: zln(:,:),diag(:,:)
4342 real(kind=kreal), intent(out) :: zz(:,:)
4343 integer(kind=kint), intent(in) :: n
4344
4345 integer(kind=kint) :: i
4346
4347 do 100 i=1,n
4348 zz(4,i)=zln(4,i)-zln(1,i)*diag(2,i)
4349 zz(7,i)=zln(7,i)-zln(1,i)*diag(4,i)-zz(4,i)*diag(5,i)
4350 zz(1,i)=zln(1,i)*diag(1,i)
4351 zz(4,i)=zz(4,i)*diag(3,i)
4352 zz(7,i)=zz(7,i)*diag(6,i)
4353 zz(4,i)=zz(4,i)-zz(7,i)*diag(5,i)
4354 zz(1,i)=zz(1,i)-zz(4,i)*diag(2,i)-zz(7,i)*diag(4,i)
4355 !
4356 zz(5,i)=zln(5,i)-zln(2,i)*diag(2,i)
4357 zz(8,i)=zln(8,i)-zln(2,i)*diag(4,i)-zz(5,i)*diag(5,i)
4358 zz(2,i)=zln(2,i)*diag(1,i)
4359 zz(5,i)=zz(5,i)*diag(3,i)
4360 zz(8,i)=zz(8,i)*diag(6,i)
4361 zz(5,i)=zz(5,i)-zz(8,i)*diag(5,i)
4362 zz(2,i)=zz(2,i)-zz(5,i)*diag(2,i)-zz(8,i)*diag(4,i)
4363 !
4364 zz(6,i)=zln(6,i)-zln(3,i)*diag(2,i)
4365 zz(9,i)=zln(9,i)-zln(3,i)*diag(4,i)-zz(6,i)*diag(5,i)
4366 zz(3,i)=zln(3,i)*diag(1,i)
4367 zz(6,i)=zz(6,i)*diag(3,i)
4368 zz(9,i)=zz(9,i)*diag(6,i)
4369 zz(6,i)=zz(6,i)-zz(9,i)*diag(5,i)
4370 zz(3,i)=zz(3,i)-zz(6,i)*diag(2,i)-zz(9,i)*diag(4,i)
4371 100 continue
4372 return
4373 end subroutine v3prod
4374 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4375
4376 subroutine d3dotl(t,a,b,n)
4377 implicit none
4378
4379 real(kind=kreal), intent(in) :: a(:,:),b(:,:)
4380 real(kind=kreal), intent(out) :: t(:)
4381 integer(kind=kint), intent(in) :: n
4382
4383 integer(kind=kint) :: l,jj
4384 ! double precision t(6),a(9,n),b(9,n)
4385 !
4386 !----------------------------------------------------------------------
4387 !
4388 ! spdot1 performs inner product of sparse vectors
4389 !
4390 !
4391 ! #coded by t.arakawa
4392 !
4393 !----------------------------------------------------------------------
4394 !
4395 !$dir max_trips(6)
4396 do 10 l=1,6
4397 t(l)=0.0d0
4398 10 continue
4399 do 100 jj=1,n
4400 t(1)=t(1)+a(1,jj)*b(1,jj)+a(4,jj)*b(4,jj)+a(7,jj)*b(7,jj)
4401 t(2)=t(2)+a(2,jj)*b(1,jj)+a(5,jj)*b(4,jj)+a(8,jj)*b(7,jj)
4402 t(3)=t(3)+a(2,jj)*b(2,jj)+a(5,jj)*b(5,jj)+a(8,jj)*b(8,jj)
4403 t(4)=t(4)+a(3,jj)*b(1,jj)+a(6,jj)*b(4,jj)+a(9,jj)*b(7,jj)
4404 t(5)=t(5)+a(3,jj)*b(2,jj)+a(6,jj)*b(5,jj)+a(9,jj)*b(8,jj)
4405 t(6)=t(6)+a(3,jj)*b(3,jj)+a(6,jj)*b(6,jj)+a(9,jj)*b(9,jj)
4406 100 continue
4407 return
4408 end subroutine d3dotl
4409
4410 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4411
4412 subroutine addr3(isw,i,j,aij,invp,xlnzr,colno,diag,zln,dsln,nstop,ir)
4413
4414 implicit none
4415
4416 integer(kind=kint), intent(in) :: invp(:),xlnzr(:),colno(:)
4417 real(kind=kreal), intent(in) :: aij(:) !zln(9,*),diag(6,*),dsln(9,*),aij(9)
4418 real(kind=kreal), intent(out) :: zln(:,:),diag(:,:),dsln(:,:) !zln(9,*),diag(6,*),dsln(9,*),aij(9)
4419 integer(kind=kint), intent(in) :: isw,i,j,nstop
4420 integer(kind=kint), intent(out) :: ir
4421
4422 integer(kind=kint), parameter :: ndeg2=9
4423 integer(kind=kint), parameter :: ndeg2l=6
4424 integer(kind=kint) :: k,l,ii,jj,itrans,i0,j0,ks,ke
4425
4426 ir=0
4427 ii=invp(i)
4428 jj=invp(j)
4429 if(ldbg) write(idbg,*) 'addr3',ii,jj,aij
4430
4431 ! diagonal
4432 if(ii.eq.jj) then
4433 diag(1,ii)=aij(1)
4434 diag(2,ii)=aij(2)
4435 diag(3,ii)=aij(5)
4436 diag(4,ii)=aij(3)
4437 diag(5,ii)=aij(6)
4438 diag(6,ii)=aij(9)
4439 goto 1000
4440 endif
4441 itrans=0
4442 if(jj.gt.ii) then
4443 k=jj
4444 jj=ii
4445 ii=k
4446 itrans=1
4447 endif
4448
4449 ! D region
4450 if(jj.ge.nstop) then
4451 i0=ii-nstop
4452 j0=jj-nstop+1
4453 k=i0*(i0-1)/2+j0
4454 if(itrans.eq.0) then
4455 do 110 l=1,ndeg2
4456 dsln(l,k)=aij(l)
4457 110 continue
4458 goto 1000
4459 else
4460 dsln(1,k)=aij(1)
4461 dsln(2,k)=aij(4)
4462 dsln(3,k)=aij(7)
4463 dsln(4,k)=aij(2)
4464 dsln(5,k)=aij(5)
4465 dsln(6,k)=aij(8)
4466 dsln(7,k)=aij(3)
4467 dsln(8,k)=aij(6)
4468 dsln(9,k)=aij(9)
4469 goto 1000
4470 endif
4471 endif
4472
4473 ! A and C region
4474 ks=xlnzr(ii)
4475 ke=xlnzr(ii+1)-1
4476 do 100 k=ks,ke
4477 if(colno(k).eq.jj) then
4478 if(itrans.eq.0) then
4479 do 120 l=1,ndeg2
4480 zln(l,k)=aij(l)
4481 120 continue
4482 else
4483 zln(1,k)=aij(1)
4484 zln(2,k)=aij(4)
4485 zln(3,k)=aij(7)
4486 zln(4,k)=aij(2)
4487 zln(5,k)=aij(5)
4488 zln(6,k)=aij(8)
4489 zln(7,k)=aij(3)
4490 zln(8,k)=aij(6)
4491 zln(9,k)=aij(9)
4492 endif
4493 goto 1000
4494 endif
4495 100 continue
4496 ir=20
4497 1000 continue
4498 return
4499 end subroutine addr3
4500
4501 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4502
4503 subroutine s2um(ic,xlnzr,colno,zln,diag,nch,par,neqns)
4504
4505 implicit none
4506
4507 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),par(:)
4508 integer(kind=kint), intent(out) :: nch(:)
4509 real(kind=kreal), intent(out) :: zln(:,:),diag(:,:) ! zln(6,*), diag(3,*)
4510 integer(kind=kint), intent(in) :: ic,neqns
4511
4512 integer(kind=kint) :: i,j,k,l,ks,ke,jj,jc,ir,kk, ierr
4513 real(kind=kreal), allocatable :: temp(:,:)
4514 integer(kind=kint), allocatable :: indx(:)
4515 real(kind=kreal) :: s(4),zz(4),t(3)
4516
4517 allocate(temp(4,neqns),indx(neqns), stat=ierr)
4518 if(ierr .ne. 0) then
4519 call errtrp('stop due to allocation error.')
4520 end if
4521
4522 ks=xlnzr(ic)
4523 ke=xlnzr(ic+1)
4524 t(1)=0.0d0
4525 t(2)=0.0d0
4526 t(3)=0.0d0
4527 do 200 k=ks,ke-1
4528 jc=colno(k)
4529 indx(jc)=ic
4530 do 210 l=1,4
4531 s(l)=0.0d0
4532 zz(l)=zln(l,k)
4533 210 continue
4534 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
4535 j=colno(jj)
4536 if(indx(j).eq.ic) then
4537 zz(1)=zz(1)-temp(1,j)*zln(1,jj)-temp(3,j)*zln(3,jj)
4538 zz(2)=zz(2)-temp(2,j)*zln(1,jj)-temp(4,j)*zln(3,jj)
4539 zz(3)=zz(3)-temp(1,j)*zln(2,jj)-temp(3,j)*zln(4,jj)
4540 zz(4)=zz(4)-temp(2,j)*zln(2,jj)-temp(4,j)*zln(4,jj)
4541 endif
4542 310 continue
4543 call inv22(zln(:,k),zz,diag(:,jc))
4544 do 220 l=1,4
4545 temp(l,jc)=zz(l)
4546 220 continue
4547 t(1)=t(1)+zz(1)*zln(1,k)+zz(3)*zln(3,k)
4548 t(2)=t(2)+zz(2)*zln(1,k)+zz(4)*zln(3,k)
4549 t(3)=t(3)+zz(2)*zln(2,k)+zz(4)*zln(4,k)
4550 200 continue
4551 diag(1,ic)=diag(1,ic)-t(1)
4552 diag(2,ic)=diag(2,ic)-t(2)
4553 diag(3,ic)=diag(3,ic)-t(3)
4554 call inv2(diag(:,ic),ir)
4555 nch(ic)=-1
4556 kk=par(ic)
4557 nch(kk)=nch(kk)-1
4558 return
4559 end subroutine s2um
4560
4561 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4562
4563 subroutine s2um1(ic,xlnzr,colno,zln,diag,nch,par,neqns)
4564
4565 implicit none
4566
4567 integer(kind=kint), intent(in) :: xlnzr(:),colno(:),nch(:),par(:)
4568 real(kind=kreal), intent(in) :: diag(:,:) !zln(4,*),diag(3,*)
4569 real(kind=kreal), intent(out) :: zln(:,:)
4570 integer(kind=kint), intent(in) :: ic,neqns
4571
4572 integer(kind=kint) :: i,j,k,l,ks,ke,jc,jj, ierr
4573 real(kind=kreal) :: s(4),zz(4)
4574 real(kind=kreal), allocatable :: temp(:,:)
4575 integer(kind=kint), allocatable :: indx(:)
4576
4577 allocate(temp(4,neqns),indx(neqns), stat=ierr)
4578 if(ierr .ne. 0) then
4579 call errtrp('stop due to allocation error.')
4580 end if
4581
4582 ks=xlnzr(ic)
4583 ke=xlnzr(ic+1)
4584 do 100 l=1,4
4585 s(l)=0.0d0
4586 100 continue
4587 do 200 k=ks,ke-1
4588 jc=colno(k)
4589 indx(jc)=ic
4590 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
4591 j=colno(jj)
4592 if(indx(j).eq.ic) then
4593 s(1)=s(1)+temp(1,j)*zln(1,jj)+temp(3,j)*zln(3,jj)
4594 s(2)=s(2)+temp(2,j)*zln(1,jj)+temp(4,j)*zln(3,jj)
4595 s(3)=s(3)+temp(1,j)*zln(2,jj)+temp(3,j)*zln(4,jj)
4596 s(4)=s(4)+temp(2,j)*zln(2,jj)+temp(4,j)*zln(4,jj)
4597 endif
4598 310 continue
4599 do 320 l=1,4
4600 temp(l,jc)=zln(l,k)-s(l)
4601 zln(l,k)=temp(l,jc)
4602 s(l)=0.0d0
4603 320 continue
4604 200 continue
4605 return
4606 end subroutine s2um1
4607
4608 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4609
4610 subroutine s2um2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln)
4611
4612 implicit none
4613
4614 integer(kind=kint), intent(in) :: neqns, nstop
4615 integer(kind=kint), intent(in) :: xlnzr(:),colno(:)
4616 real(kind=kreal), intent(inout) :: zln(:,:),diag(:,:)!zln(4,*),diag(3,*)
4617 integer(kind=kint), pointer :: spdslnidx(:)
4618 real(kind=kreal), pointer :: spdslnval(:,:)
4619 integer(kind=kint), intent(out) :: nspdsln
4620
4621 integer(kind=kint) :: i,j,k,l,m,n, ic,ks,ke,ii,jj,jc,j1,j2,loc,ispdsln, ierr
4622 real(kind=kreal), allocatable :: temp(:,:)
4623 integer(kind=kint), allocatable :: indx(:)
4624 logical :: ftflag
4625
4626 allocate(temp(4,neqns),indx(neqns), stat=ierr)
4627 if(ierr .ne. 0) then
4628 call errtrp('stop due to allocation error.')
4629 end if
4630
4631 nspdsln=0
4632 do ic=nstop,neqns
4633 ks=xlnzr(ic)
4634 ke=xlnzr(ic+1)-1
4635 do k=ks,ke
4636 jj=colno(k)
4637 indx(jj)=ic
4638 end do
4639 do jc=nstop,ic-1
4640 j1=xlnzr(jc)
4641 j2=xlnzr(jc+1)
4642 do jj=xlnzr(jc),xlnzr(jc+1)-1
4643 j=colno(jj)
4644 if(indx(j).eq.ic) then
4645 nspdsln=nspdsln+1
4646 exit
4647 endif
4648 end do
4649 end do
4650 end do
4651 allocate(spdslnidx(nspdsln),spdslnval(4,nspdsln), stat=ierr)
4652 if(ierr .ne. 0) then
4653 call errtrp('stop due to allocation error.')
4654 end if
4655
4656 loc=0
4657 ispdsln=0
4658 spdslnval=0
4659 ftflag = .true.
4660 do 100 ic=nstop,neqns
4661 ks=xlnzr(ic)
4662 ke=xlnzr(ic+1)-1
4663 do 110 k=ks,ke
4664 jj=colno(k)
4665 temp(1,jj)=zln(1,k)
4666 temp(2,jj)=zln(2,k)
4667 temp(3,jj)=zln(3,k)
4668 temp(4,jj)=zln(4,k)
4669 ! call inv22(zln(1,k),temp(1,jj),diag(1,jj))
4670 zln(3,k)=temp(3,jj)-temp(1,jj)*diag(2,jj)
4671 zln(1,k)=temp(1,jj)*diag(1,jj)
4672 zln(3,k)=zln(3,k)*diag(3,jj)
4673 zln(1,k)=zln(1,k)-zln(3,k)*diag(2,jj)
4674 !
4675 zln(4,k)=temp(4,jj)-temp(2,jj)*diag(2,jj)
4676 zln(2,k)=temp(2,jj)*diag(1,jj)
4677 zln(4,k)=zln(4,k)*diag(3,jj)
4678 zln(2,k)=zln(2,k)-zln(4,k)*diag(2,jj)
4679 !
4680 diag(1,ic)=diag(1,ic)-(temp(1,jj)*zln(1,k)+temp(3,jj)*zln(3,k))
4681 diag(2,ic)=diag(2,ic)-(temp(1,jj)*zln(2,k)+temp(3,jj)*zln(4,k))
4682 diag(3,ic)=diag(3,ic)-(temp(2,jj)*zln(2,k)+temp(4,jj)*zln(4,k))
4683 indx(jj)=ic
4684 110 continue
4685 do 120 jc=nstop,ic-1
4686 loc=loc+1
4687 do 220 jj=xlnzr(jc),xlnzr(jc+1)-1
4688 j=colno(jj)
4689 if(indx(j).eq.ic) then
4690 if (ftflag) then
4691 ispdsln=ispdsln+1
4692 ftflag=.false.
4693 end if
4694 spdslnidx(ispdsln)=loc
4695 spdslnval(1,ispdsln)=spdslnval(1,ispdsln)-(temp(1,j)*zln(1,jj)+temp(3,j)*zln(3,jj))
4696 spdslnval(2,ispdsln)=spdslnval(2,ispdsln)-(temp(2,j)*zln(1,jj)+temp(4,j)*zln(3,jj))
4697 spdslnval(3,ispdsln)=spdslnval(3,ispdsln)-(temp(1,j)*zln(2,jj)+temp(3,j)*zln(4,jj))
4698 spdslnval(4,ispdsln)=spdslnval(4,ispdsln)-(temp(2,j)*zln(2,jj)+temp(4,j)*zln(4,jj))
4699 endif
4700 220 continue
4701 ftflag = .true.
4702 120 continue
4703 100 continue
4704 return
4705 end subroutine s2um2_child
4706
4707 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4708
4709 subroutine inv22(zln,zz,diag)
4710
4711 implicit none
4712
4713 real(kind=kreal), intent(in) :: zz(4),diag(3)
4714 real(kind=kreal), intent(out) :: zln(4)
4715
4716 zln(3)=zz(3)-zz(1)*diag(2)
4717 zln(1)=zz(1)*diag(1)
4718 zln(3)=zln(3)*diag(3)
4719 zln(1)=zln(1)-zln(3)*diag(2)
4720
4721 zln(4)=zz(4)-zz(2)*diag(2)
4722 zln(2)=zz(2)*diag(1)
4723 zln(4)=zln(4)*diag(3)
4724 zln(2)=zln(2)-zln(4)*diag(2)
4725
4726 return
4727 end subroutine inv22
4728
4729 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4730
4731 subroutine inv2(dsln,ir)
4732
4733 implicit none
4734
4735 real(kind=kreal), intent(out) :: dsln(3)
4736 integer(kind=kint), intent(out) :: ir
4737
4738 real(kind=kreal) :: t
4739
4740 ir=0
4741 if(dabs(dsln(1)).lt.rmin) then
4742 ir=10
4743 return
4744 endif
4745 dsln(1)=1.0d0/dsln(1)
4746 t=dsln(2)*dsln(1)
4747 dsln(3)=dsln(3)-t*dsln(2)
4748 dsln(2)=t
4749 if(dabs(dsln(3)).lt.rmin) then
4750 ir=10
4751 return
4752 endif
4753 dsln(3)=1.0d0/dsln(3)
4754 return
4755 end subroutine inv2
4756
4757 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4758
4759 subroutine d2sdot(wi,a,b,n)
4760
4761 implicit none
4762
4763 real(kind=kreal), intent(in) :: a(:,:),b(:,:)!wi(2),a(2,*),b(4,*)
4764 real(kind=kreal), intent(inout) :: wi(:)
4765 integer(kind=kint), intent(in) :: n
4766
4767 integer(kind=kint) :: jj
4768 !
4769 !----------------------------------------------------------------------
4770 !
4771 ! d2sdot performs inner product of dens vectors
4772 !
4773 !
4774 ! #coded by t.arakawa
4775 !
4776 !----------------------------------------------------------------------
4777 !
4778 do 100 jj=1,n
4779 wi(1)=wi(1)-a(1,jj)*b(1,jj)-a(2,jj)*b(3,jj)
4780 wi(2)=wi(2)-a(1,jj)*b(2,jj)-a(2,jj)*b(4,jj)
4781 100 continue
4782 return
4783 end subroutine d2sdot
4784
4785 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4786
4787 subroutine s2pdot(bi,b,zln,colno,ks,ke)
4788
4789 implicit none
4790
4791 integer(kind=kint), intent(in) :: colno(:)
4792 integer(kind=kint), intent(in) :: ks,ke
4793 real(kind=kreal), intent(in) :: zln(:,:),b(:,:) !zln(4,*),b(2,*),bi(2)
4794 real(kind=kreal), intent(out) :: bi(:) !zln(4,*),b(2,*),bi(2)
4795
4796 integer(kind=kint) :: jj,j
4797
4798 !----------------------------------------------------------------------
4799 !
4800 ! s2pdot performs inner product of sparse vectors
4801 !
4802 !
4803 ! #coded by t.arakawa
4804 !
4805 !----------------------------------------------------------------------
4806
4807 do 100 jj=ks,ke
4808 j=colno(jj)
4809 bi(1)=bi(1)-zln(1,jj)*b(1,j)-zln(3,jj)*b(2,j)
4810 bi(2)=bi(2)-zln(2,jj)*b(1,j)-zln(4,jj)*b(2,j)
4811 100 continue
4812 return
4813 end subroutine s2pdot
4814
4815 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4816
4817 subroutine addrx(isw,i,j,aij,invp,xlnzr,colno,diag,zln,dsln,nstop,ndeg,ndeg2,ndeg2l,ir)
4818
4819 implicit none
4820
4821 integer(kind=kint), intent(in) :: invp(*),xlnzr(*),colno(*)
4822 real(kind=kreal), intent(in) :: aij(ndeg,ndeg)
4823 real(kind=kreal), intent(out) :: zln(ndeg,ndeg,*),diag(ndeg2l,*),dsln(ndeg,ndeg,*)
4824 integer(kind=kint), intent(in) :: isw,i,j,nstop,ndeg,ndeg2,ndeg2l
4825 integer(kind=kint), intent(out) :: ir
4826
4827 integer(kind=kint) :: ii,jj,k,l,m,n,ks,ke,itrans,i0,j0
4828
4829 ir=0
4830 ii=invp(i)
4831 jj=invp(j)
4832 if(ldbg) write(idbg,*) 'addrx',ii,jj,aij
4833 if(ii.eq.jj) then
4834 l=0
4835 do 100 n=1,ndeg
4836 do 110 m=1,n
4837 l=l+1
4838 diag(l,ii)=aij(n,m)
4839 110 continue
4840 100 continue
4841 goto 1000
4842 endif
4843 itrans=0
4844 if(jj.gt.ii) then
4845 k=jj
4846 jj=ii
4847 ii=k
4848 itrans=1
4849 endif
4850 if(jj.ge.nstop) then
4851 i0=ii-nstop
4852 j0=jj-nstop+1
4853 k=i0*(i0-1)/2+j0
4854 if(itrans.eq.0) then
4855 do 120 m=1,ndeg
4856 do 130 n=1,ndeg
4857 dsln(n,m,k)=aij(n,m)
4858 130 continue
4859 120 continue
4860 goto 1000
4861 else
4862 do 140 m=1,ndeg
4863 do 150 n=1,ndeg
4864 dsln(n,m,k)=aij(m,n)
4865 150 continue
4866 140 continue
4867 goto 1000
4868 endif
4869 endif
4870 ks=xlnzr(ii)
4871 ke=xlnzr(ii+1)-1
4872 do 200 k=ks,ke
4873 if(colno(k).eq.jj) then
4874 if(itrans.eq.0) then
4875 do 160 m=1,ndeg
4876 do 170 n=1,ndeg
4877 zln(n,m,k)=aij(n,m)
4878 170 continue
4879 160 continue
4880 else
4881 do 180 m=1,ndeg
4882 do 190 n=1,ndeg
4883 zln(n,m,k)=aij(m,n)
4884 190 continue
4885 180 continue
4886 endif
4887 goto 1000
4888 endif
4889 200 continue
4890 ir=20
4891 1000 continue
4892 return
4893 end subroutine addrx
4894
4895 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4896
4897 subroutine dxdot(ndeg,t,a,b,l)
4898
4899 implicit none
4900
4901 real(kind=kreal), intent(in) :: a(ndeg,ndeg,*),b(ndeg,ndeg,*)
4902 real(kind=kreal), intent(out) :: t(ndeg,ndeg)
4903 integer(kind=kint), intent(in) :: ndeg,l
4904
4905 integer(kind=kint) :: k,jj,n,m
4906 !
4907 !----------------------------------------------------------------------
4908 !
4909 ! spdot1 performs inner product of sparse vectors
4910 !
4911 !
4912 ! #coded by t.arakawa
4913 !
4914 !----------------------------------------------------------------------
4915 !
4916 do 221 n=1,ndeg
4917 do 221 m=1,ndeg
4918 t(n,m)=0.0d0
4919 do 221 k=1,ndeg
4920 do 100 jj=1,l
4921 t(n,m)=t(n,m)+a(n,k,jj)*b(m,k,jj)
4922 100 continue
4923 221 continue
4924 return
4925 end subroutine dxdot
4926
4927 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4928
4929 subroutine dxdotl(ndeg,t,a,b,l)
4930
4931 implicit none
4932
4933 real(kind=kreal), intent(in) :: a(ndeg,ndeg,*),b(ndeg,ndeg,*)
4934 real(kind=kreal), intent(out) :: t(ndeg,ndeg)
4935 integer(kind=kint), intent(in) :: ndeg,l
4936
4937 integer(kind=kint) :: n,m,jj,k
4938 !
4939 !----------------------------------------------------------------------
4940 !
4941 ! spdot1 performs inner product of sparse vectors
4942 !
4943 !
4944 ! #coded by t.arakawa
4945 !
4946 !----------------------------------------------------------------------
4947 !
4948 do 221 n=1,ndeg
4949 do 221 m=1,n
4950 t(n,m)=0.0d0
4951 do 221 k=1,ndeg
4952 do 100 jj=1,l
4953 t(n,m)=t(n,m)+a(n,k,jj)*b(m,k,jj)
4954 100 continue
4955 221 continue
4956 return
4957 end subroutine dxdotl
4958
4959 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4960
4961 subroutine dxsdot(ndeg,wi,a,b,n)
4962
4963 implicit none
4964
4965 real(kind=kreal), intent(in) :: a(ndeg,*),b(ndeg,ndeg,*)
4966 real(kind=kreal), intent(out) :: wi(ndeg)
4967 integer(kind=kint), intent(in) :: ndeg, n
4968
4969 integer(kind=kint) :: jj, k, l
4970 !
4971 !----------------------------------------------------------------------
4972 !
4973 ! dxsdot performs inner product of dens vectors
4974 !
4975 !
4976 ! #coded by t.arakawa
4977 ! #reviced by t.kitayama 20071122
4978 !
4979 !----------------------------------------------------------------------
4980 !
4981 do jj=1,n
4982 do k=1,ndeg
4983 do l=1,ndeg
4984 wi(l)=wi(l)-b(l,k,jj)*a(k,jj)
4985 end do
4986 end do
4987 end do
4988 return
4989 end subroutine dxsdot
4990
4991 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4992
4993 subroutine invx(dsln,ndeg,ir)
4994
4995 implicit none
4996
4997 real(kind=kreal), intent(inout) :: dsln(*)
4998 integer(kind=kint), intent(in) :: ndeg
4999 integer(kind=kint), intent(out) :: ir
5000
5001 integer(kind=kint) :: i,j,k,l,ld,l0,k0,ll
5002 real(kind=kreal) :: tem,t
5003
5004 ir=0
5005 l=1
5006 dsln(1)=1.0d0/dsln(1)
5007 do 100 i=2,ndeg
5008 ld=0
5009 l0=l
5010 do 110 j=1,i-1
5011 l=l+1
5012 do 120 k=1,j-1
5013 ld=ld+1
5014 dsln(l)=dsln(l)-dsln(l0+k)*dsln(ld)
5015 120 continue
5016 ld=ld+1
5017 110 continue
5018 t=0.0d0
5019 k0=0
5020 ll=0
5021 do 130 k=l-i+2,l
5022 ll=ll+1
5023 k0=k0+ll
5024 tem=dsln(k)*dsln(k0)
5025 t=t+tem*dsln(k)
5026 dsln(k)=tem
5027 130 continue
5028 l=l+1
5029 dsln(l)=dsln(l)-t
5030 dsln(l)=1.0d0/dsln(l)
5031 100 continue
5032 return
5033 end subroutine invx
5034
5035 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5036
5037 subroutine invxx(zln,zz,diag,ndeg)
5038
5039 implicit none
5040
5041 real(kind=kreal), intent(in) :: zz(ndeg,ndeg),diag(*)
5042 real(kind=kreal), intent(out) :: zln(ndeg,ndeg)
5043 integer(kind=kint), intent(in) :: ndeg
5044
5045 integer(kind=kint) :: i,j,k,l,m,n,loc,loc1
5046
5047 zln=zz
5048 do 100 l=1,ndeg
5049 loc=0
5050 do 120 m=1,ndeg
5051 loc=loc+m
5052 loc1=loc+m
5053 do 120 n=m+1,ndeg
5054 zln(l,n)=zln(l,n)-zln(l,m)*diag(loc1)
5055 loc1=loc1+n
5056 120 continue
5057 loc=0
5058 do 130 m=1,ndeg
5059 loc=loc+m
5060 zln(l,m)=zln(l,m)*diag(loc)
5061 130 continue
5062 do 140 n=ndeg,1,-1
5063 loc=loc-1
5064 do 140 m=n-1,1,-1
5065 zln(l,m)=zln(l,m)-zln(l,n)*diag(loc)
5066 loc=loc-1
5067 140 continue
5068 100 continue
5069 return
5070 end subroutine invxx
5071
5072 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5073
5074 subroutine sxpdot(ndeg,bi,b,zln,colno,ks,ke)
5075
5076 implicit none
5077
5078 integer(kind=kint), intent(in) :: colno(*)
5079 real(kind=kreal), intent(in) :: zln(ndeg,ndeg,*),b(ndeg,*)
5080 real(kind=kreal), intent(out) :: bi(ndeg)
5081 integer(kind=kint),intent(in) :: ndeg,ks,ke
5082
5083 integer(kind=kint) :: j,jj,m,n
5084 !
5085 !----------------------------------------------------------------------
5086 !
5087 ! sxpdot performs inner product of sparse vectors
5088 !
5089 !
5090 ! #coded by t.arakawa
5091 !
5092 !----------------------------------------------------------------------
5093 !
5094 do 100 jj=ks,ke
5095 j=colno(jj)
5096 do 100 m=1,ndeg
5097 do 100 n=1,ndeg
5098 bi(n)=bi(n)-zln(n,m,jj)*b(m,j)
5099 100 continue
5100 return
5101 end subroutine sxpdot
5102
5103 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5104
5105 subroutine sxum(ic,xlnzr,colno,zln,diag,nch,par,neqns,ndeg,ndegl)
5106
5107 implicit none
5108
5109 integer(kind=kint), intent(in) :: xlnzr(*),colno(*),par(*)
5110 integer(kind=kint), intent(out) :: nch(*)
5111 real(kind=kreal), intent(out) :: zln(ndeg,ndeg,*),diag(ndegl,*)
5112 integer(kind=kint), intent(in) :: ic,neqns,ndeg,ndegl
5113
5114 real(kind=kreal) :: zz(ndeg,ndeg),t(ndegl)
5115 integer(kind=kint) :: i,j,k,l,m,n,ndeg22,ks,ke,jc,loc,jj,kk,ir, ierr
5116 real(kind=kreal),allocatable :: temp(:,:,:)
5117 integer(kind=kint),allocatable :: indx(:)
5118
5119 ndeg22=ndeg*ndeg
5120 allocate(temp(ndeg,ndeg,neqns),indx(neqns), stat=ierr)
5121 if(ierr .ne. 0) then
5122 call errtrp('stop due to allocation error.')
5123 end if
5124
5125 ks=xlnzr(ic)
5126 ke=xlnzr(ic+1)
5127 t=0.0
5128
5129 do 200 k=ks,ke-1
5130 jc=colno(k)
5131 indx(jc)=ic
5132 zz=zln(:,:,k)
5133 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
5134 j=colno(jj)
5135 if(indx(j).eq.ic) then
5136 do 311 m=1,ndeg
5137 do 311 n=1,ndeg
5138 do 311 kk=1,ndeg
5139 zz(n,m)=zz(n,m)-temp(n,kk,j)*zln(m,kk,jj)
5140
5141 311 continue
5142 endif
5143 310 continue
5144 call invxx(zln(1,1,k),zz,diag(1,jc),ndeg)
5145 temp(:,:,jc)=zz
5146 loc=0
5147 do 221 n=1,ndeg
5148 do 221 m=1,n
5149 loc=loc+1
5150 do 221 kk=1,ndeg
5151 t(loc)=t(loc)+zz(n,kk)*zln(m,kk,k)
5152 221 continue
5153 200 continue
5154 diag(:,ic)=diag(:,ic)-t
5155 call invx(diag(1,ic),ndeg,ir)
5156 nch(ic)=-1
5157 kk=par(ic)
5158 nch(kk)=nch(kk)-1
5159 return
5160 end subroutine sxum
5161
5162 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5163
5164 subroutine sxum1(ic,xlnzr,colno,zln,diag,nch,par,neqns,ndeg,ndegl)
5165
5166 implicit none
5167
5168 integer(kind=kint), intent(in) :: xlnzr(*),colno(*),nch(*),par(*)
5169 real(kind=kreal), intent(in) :: diag(ndegl,*)
5170 real(kind=kreal), intent(out) :: zln(ndeg,ndeg,*)
5171 integer(kind=kint), intent(in) :: ic,neqns,ndeg,ndegl
5172
5173 real(kind=kreal) :: s(ndeg,ndeg)
5174 integer(kind=kint) :: i,j,k,l,m,n,ks,ke,jc,jj,kk, ierr
5175 real(kind=kreal),allocatable :: temp(:,:,:)
5176 integer(kind=kint),allocatable :: indx(:)
5177
5178 allocate(temp(ndeg,ndeg,neqns),indx(neqns), stat=ierr)
5179 if(ierr .ne. 0) then
5180 call errtrp('stop due to allocation error.')
5181 end if
5182 ks=xlnzr(ic)
5183 ke=xlnzr(ic+1)
5184 do 100 m=1,ndeg
5185 do 100 n=1,ndeg
5186 s(n,m)=0.0d0
5187 100 continue
5188 do 200 k=ks,ke-1
5189 jc=colno(k)
5190 indx(jc)=ic
5191 do 310 jj=xlnzr(jc),xlnzr(jc+1)-1
5192 j=colno(jj)
5193 if(indx(j).eq.ic) then
5194 do 311 m=1,ndeg
5195 do 311 n=1,ndeg
5196 do 311 kk=1,ndeg
5197 s(n,m)=s(n,m)+temp(n,kk,j)*zln(m,kk,jj)
5198 311 continue
5199 endif
5200 310 continue
5201 do 320 m=1,ndeg
5202 do 320 n=1,ndeg
5203 temp(n,m,jc)=zln(n,m,k)-s(n,m)
5204 zln(n,m,k)=temp(n,m,jc)
5205 s(n,m)=0.0d0
5206 320 continue
5207 200 continue
5208 return
5209 end subroutine sxum1
5210
5211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5212
5213 subroutine sxum2_child(neqns,nstop,xlnzr,colno,zln,diag,spdslnidx,spdslnval,nspdsln,ndeg,ndegl)
5214
5215 implicit none
5216
5217 integer(kind=kint), intent(in) :: neqns, nstop
5218 integer(kind=kint), intent(in) :: xlnzr(*),colno(*)
5219 real(kind=kreal), intent(inout) :: zln(ndeg,ndeg,*),diag(ndegl,*)
5220 integer(kind=kint), pointer :: spdslnidx(:)
5221 real(kind=kreal), pointer :: spdslnval(:,:)
5222 integer(kind=kint), intent(out) :: nspdsln
5223 integer(kind=kint), intent(in) :: ndeg, ndegl
5224
5225 integer(kind=kint) :: i,j,k,l,m,n, ic,ks,ke,ii,jj,jc,j1,j2,loc,locd,kk, ierr
5226 integer(kind=kint) :: ispdsln
5227 real(kind=kreal), allocatable :: temp(:,:,:)
5228 integer(kind=kint), allocatable :: indx(:)
5229 logical :: ftflag
5230
5231 allocate(temp(ndeg,ndeg,neqns),indx(neqns), stat=ierr)
5232 if(ierr .ne. 0) then
5233 call errtrp('stop due to allocation error.')
5234 end if
5235
5236 nspdsln=0
5237 do ic=nstop,neqns
5238 ks=xlnzr(ic)
5239 ke=xlnzr(ic+1)-1
5240 do k=ks,ke
5241 jj=colno(k)
5242 indx(jj)=ic
5243 end do
5244 do jc=nstop,ic-1
5245 j1=xlnzr(jc)
5246 j2=xlnzr(jc+1)
5247 do jj=xlnzr(jc),xlnzr(jc+1)-1
5248 j=colno(jj)
5249 if(indx(j).eq.ic) then
5250 nspdsln=nspdsln+1
5251 exit
5252 endif
5253 end do
5254 end do
5255 end do
5256 allocate(spdslnidx(nspdsln),spdslnval(ndeg*ndeg,nspdsln), stat=ierr)
5257 if(ierr .ne. 0) then
5258 call errtrp('stop due to allocation error.')
5259 end if
5260
5261 loc=0
5262 ispdsln=0
5263 spdslnval=0
5264 ftflag = .true.
5265 do 100 ic=nstop,neqns
5266 ks=xlnzr(ic)
5267 ke=xlnzr(ic+1)-1
5268 do 110 k=ks,ke
5269 jj=colno(k)
5270 do 110 m=1,ndeg
5271 do 110 n=1,ndeg
5272 temp(n,m,jj)=zln(n,m,k)
5273 indx(jj)=ic
5274 110 continue
5275 do 111 k=ks,ke
5276 jj=colno(k)
5277 call invxx(zln(1,1,k),temp(1,1,jj),diag(1,jj),ndeg)
5278 111 continue
5279 !
5280 locd=0
5281 do 112 n=1,ndeg
5282 do 112 m=1,n
5283 locd=locd+1
5284 do 112 k=ks,ke
5285 jj=colno(k)
5286 do 112 kk=1,ndeg
5287 diag(locd,ic)=diag(locd,ic)-temp(n,kk,jj)*zln(m,kk,k)
5288 112 continue
5289 do 120 jc=nstop,ic-1
5290 loc=loc+1
5291 j1=xlnzr(jc)
5292 j2=xlnzr(jc+1)
5293 do 220 jj=xlnzr(jc),xlnzr(jc+1)-1
5294 j=colno(jj)
5295 if(indx(j).eq.ic) then
5296 if (ftflag) then
5297 ispdsln=ispdsln+1
5298 ftflag=.false.
5299 end if
5300 spdslnidx(ispdsln)=loc
5301 do 221 m=1,ndeg
5302 do 221 n=1,ndeg
5303 do 221 k=1,ndeg
5304 spdslnval(ndeg*(m-1)+n,ispdsln)=spdslnval(ndeg*(m-1)+n,ispdsln)-temp(n,k,j)*zln(m,k,jj)
5305 221 continue
5306 endif
5307 220 continue
5308 ftflag = .true.
5309 120 continue
5310 100 continue
5311 return
5312 end subroutine sxum2_child
5313
5314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5315
5316 subroutine sxum3(neqns,dsln,diag,ndeg,ndegl)
5317
5318 implicit none
5319
5320 real(kind=kreal), intent(inout):: dsln(ndeg,ndeg,*),diag(ndegl,*)
5321 integer(kind=kint), intent(in) :: neqns, ndeg, ndegl
5322
5323 integer(kind=kint) :: loc, locd, ir, i,j,n,m, ierr
5324 integer(kind=kint), allocatable :: indx(:)
5325 real(kind=kreal), allocatable :: temp(:,:,:)
5326 real(kind=kreal), allocatable :: t(:,:)
5327
5328 allocate(indx(neqns),temp(ndeg,ndeg,neqns),t(ndeg,ndeg), stat=ierr)
5329 if(ierr .ne. 0) then
5330 call errtrp('stop due to allocation error.')
5331 end if
5332
5333 if(neqns.le.0) goto 1000
5334 indx(1)=1 ! it will work...
5335 loc=1
5336 call invx(diag(1,1),ndeg,ir)
5337 do 100 i=2,neqns
5338 indx(i)=loc
5339 do 110 j=1,i-1
5340 call dxdot(ndeg,t,dsln(1,1,indx(i)),dsln(1,1,indx(j)),j-1)
5341 do 111 m=1,ndeg
5342 do 111 n=1,ndeg
5343 dsln(n,m,loc)=dsln(n,m,loc)-t(n,m)
5344 111 continue
5345 loc=loc+1
5346 110 continue
5347 call vxprod(ndeg,ndegl,dsln(1,1,indx(i)),diag,temp,i-1)
5348 call dxdotl(ndeg,t,temp,dsln(1,1,indx(i)),i-1)
5349 locd=0
5350 do 221 n=1,ndeg
5351 do 221 m=1,n
5352 locd=locd+1
5353 diag(locd,i)=diag(locd,i)-t(n,m)
5354 221 continue
5355 call vcopy(temp,dsln(1,1,indx(i)),ndeg*ndeg*(i-1))
5356 call invx(diag(1,i),ndeg,ir)
5357 100 continue
5358 1000 continue
5359 return
5360 end subroutine sxum3
5361
5362 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5363 subroutine vcopy(a,c,n)
5364 implicit none
5365
5366 integer(kind=kint) :: n
5367 real(kind=kreal) :: a(n),c(n)
5368 ! do 100 i=1,n
5369 ! c(i)=a(i)
5370 ! 100 continue
5371 c=a
5372 return
5373 end subroutine vcopy
5374
5375 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5376
5377 subroutine verif0(neqns,ndeg,nttbr,irow,jcol,val,rhs,x)
5378
5379 implicit none
5380
5381 integer(kind=kint), intent(in) :: irow(*),jcol(*)
5382 integer(kind=kint), intent(in) :: neqns,ndeg,nttbr
5383 real(kind=kreal), intent(in) :: val(ndeg,ndeg,*),x(ndeg,*)
5384 real(kind=kreal), intent(out) :: rhs(ndeg,*)
5385
5386 integer(kind=kint) :: i,j,k,l,m
5387 real(kind=kreal) :: rel,err
5388 !
5389 !----------------------------------------------------------------------
5390 !
5391 ! verify the solution(symmetric matrix)
5392 !
5393 !----------------------------------------------------------------------
5394 !
5395 rel=0.0d0
5396 do 10 i=1,neqns
5397 do 10 l=1,ndeg
5398 rel=rel+dabs(rhs(l,i))
5399 10 continue
5400 do 100 k=1,nttbr
5401 i=irow(k)
5402 j=jcol(k)
5403 do 101 l=1,ndeg
5404 do 102 m=1,ndeg
5405 rhs(l,i)=rhs(l,i)-val(l,m,k)*x(m,j)
5406 if(i.ne.j) rhs(l,j)=rhs(l,j)-val(m,l,k)*x(m,i)
5407 102 continue
5408 101 continue
5409 100 continue
5410 err=0.0d0
5411 do 200 i=1,neqns
5412 do 200 l=1,ndeg
5413 err=err+dabs(rhs(l,i))
5414 200 continue
5415 if (m_pds_procinfo%myid .eq. 0) then
5416 write(imsg,6000) err,rel,err/rel
5417 end if
5418 6000 format(' ***verification***(symmetric)'/&
5419 & 'norm(Ax-b) = ',1pd20.10/&
5420 & 'norm(b) = ',1pd20.10/&
5421 & 'norm(Ax-b)/norm(b) = ',1pd20.10)
5422 6010 format(1p4d15.7)
5423 return
5424 end subroutine verif0
5425
5426 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5427
5428 subroutine vxprod(ndeg,ndegl,zln,diag,zz,n)
5429
5430 implicit none
5431
5432 real(kind=kreal), intent(in) :: zln(ndeg*ndeg,n),diag(ndegl,n)
5433 real(kind=kreal), intent(out) :: zz(ndeg*ndeg,n)
5434 integer(kind=kint), intent(in) :: ndeg,ndegl,n
5435
5436 integer(kind=kint) :: i
5437
5438 do 100 i=1,n
5439 call invxx(zz(1,i),zln(1,i),diag(1,i),ndeg)
5440 100 continue
5441 return
5442 end subroutine vxprod
5443
5444 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5445
5446#endif
5447
subroutine verif0(neqns, ndeg, nttbr, irow, jcol, val, rhs, x)
subroutine idntty(neqns, invp, iperm)
subroutine, public hecmw_mat_dump(hecmat, hecmesh)
subroutine, public hecmw_mat_dump_solution(hecmat)
subroutine, public hecmw_solve_direct_parallel(hecmesh, hecmat, ii)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=kint) function hecmw_comm_get_comm()
integer(kind=4), parameter kreal
subroutine hecmw_abort(comm)
subroutine, public symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
Definition: m_elap.F90:7
subroutine, public initelap(t, i)
Definition: m_elap.F90:22
subroutine, public elapout(mes)
Definition: m_elap.F90:35
subroutine, public reovec(r, iperm)
subroutine, public matrix_partition_recursive_bisection(a0, ndiv, pmi)