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