grid_tools  1.13.0
pmat2.f90
Go to the documentation of this file.
1 
4 
21 module pmat2
22 !============================================================================
23 use pkind, only: spi,sp,dp,dpc
24 implicit none
25 private
26 public:: avco
27 public:: dfco
28 public:: dfco2
29 public:: clipb
30 public:: cad1b
31 public:: csb1b
32 public:: cad2b
33 public:: csb2b
34 public:: ldub
35 public:: ldltb
36 public:: udlb
37 public:: l1ubb
38 public:: l1ueb
39 public:: ltdlbv
40 public:: udlbv
41 public:: udlbx
42 public:: udlby
43 public:: udlvb
44 public:: udlxb
45 public:: udlyb
46 public:: u1lbv
47 public:: u1lbx
48 public:: u1lby
49 public:: u1lvb
50 public:: u1lxb
51 public:: u1lyb
52 public:: linbv
53 public:: wrtb
54 real(dp),parameter:: zero=0
55 
56 interface avco; module procedure avco, davco, tavco; end interface
57 interface dfco; module procedure dfco, ddfco, tdfco; end interface
58 interface dfco2; module procedure dfco2, ddfco2, tdfco2; end interface
59 interface clipb; module procedure clib, clib_d, clib_c; end interface
60 interface cad1b; module procedure cad1b; end interface
61 interface csb1b; module procedure csb1b; end interface
62 interface cad2b; module procedure cad2b; end interface
63 interface csb2b; module procedure csb2b; end interface
64 interface ldub; module procedure ldub, dldub; end interface
65 interface ldltb; module procedure ldltb, dldltb; end interface
66 interface l1ubb; module procedure l1ubb, dl1ubb; end interface
67 interface l1ueb; module procedure l1ueb, dl1ueb; end interface
68 interface ltdlbv; module procedure ltdlbv,dltdlbv; end interface
69 interface udlb; module procedure udlb, dudlb; end interface
70 interface udlbv; module procedure udlbv, dudlbv; end interface
71 interface udlbx; module procedure udlbx; end interface
72 interface udlby; module procedure udlby; end interface
73 interface udlvb; module procedure udlvb; end interface
74 interface udlxb; module procedure udlxb; end interface
75 interface udlyb; module procedure udlyb; end interface
76 interface u1lbv; module procedure u1lbv; end interface
77 interface u1lbx; module procedure u1lbx; end interface
78 interface u1lby; module procedure u1lby; end interface
79 interface u1lvb; module procedure u1lvb; end interface
80 interface u1lxb; module procedure u1lxb; end interface
81 interface u1lyb; module procedure u1lyb; end interface
82 interface linbv; module procedure linbv; end interface
83 interface wrtb; module procedure wrtb; end interface
84 contains
85 
99 subroutine avco(na,nb,za,zb,z0,a,b) ! [AVCO]
100 use pietc, only: u0,u1
101 use pmat, only: inv
102 implicit none
103 integer(spi),intent(in ):: na,nb
104 real(sp), intent(in ):: za(na),zb(nb),z0
105 real(sp), intent(out):: a(na),b(nb)
106 !-----------------------------------------------------------------------------
107 integer(spi) :: na1,nab,i
108 real(sp), dimension(na+nb,na+nb):: w
109 real(sp), dimension(na) :: za0,pa
110 real(sp), dimension(nb) :: zb0,pb
111 real(sp), dimension(na+nb) :: ab
112 !=============================================================================
113 na1=na+1; nab=na+nb
114 za0=za-z0; zb0=zb-z0
115 pa=u1; pb=-u1
116 w=u0; ab=u0
117 w(1,1:na)=u1; ab(1)=u1
118 do i=2,nab; w(i,1:na)=pa; pa=pa*za0; w(i,na1:nab)=pb; pb=pb*zb0; enddo
119 call inv(w,ab)
120 a=ab(1:na); b=ab(na1:nab)
121 end subroutine avco
122 
133 subroutine davco(na,nb,za,zb,z0,a,b) ! [AVCO]
134 use pietc, only: u0,u1
135 use pmat, only: inv
136 implicit none
137 integer(spi),intent(IN ):: na,nb
138 real(dp), intent(IN ):: za(na),zb(nb),z0
139 real(dp), intent(OUT):: a(na),b(nb)
140 !-----------------------------------------------------------------------------
141 integer(spi) :: na1,nab,i
142 real(dp),dimension(na+nb,na+nb):: w
143 real(dp),dimension(na) :: za0,pa
144 real(dp),dimension(nb) :: zb0,pb
145 real(dp),dimension(na+nb) :: ab
146 !=============================================================================
147 na1=na+1; nab=na+nb
148 za0=za-z0; zb0=zb-z0
149 pa=u1; pb=-u1
150 w=u0; ab=u0
151 w(1,1:na)=u1; ab(1)=u1
152 do i=2,nab; w(i,1:na)=pa; pa=pa*za0; w(i,na1:nab)=pb; pb=pb*zb0; enddo
153 call inv(w,ab)
154 a=ab(1:na); b=ab(na1:nab)
155 end subroutine davco
156 
164 subroutine tavco(xa,xb,a,b)! [AVCO]
165 implicit none
166 real(dp),dimension(:),intent(IN ):: xa,xb
167 real(dp),dimension(:),intent(OUT):: a,b
168 !-----------------------------------------------------------------------------
169 integer(spi):: na,nb
170 !=============================================================================
171 na=size(xa); if(na /= size(a))stop 'In tavco; sizes of a and xa different'
172 nb=size(xb); if(nb /= size(b))stop 'In tavco; sizes of b and xb different'
173 call davco(na,nb,xa,xb,zero,a,b)
174 end subroutine tavco
175 
189 subroutine dfco(na,nb,za,zb,z0,a,b)! [DFCO]
190 use pietc_s, only: u0,u1
191 use pmat, only: inv
192 implicit none
193 integer(spi),intent(IN ) :: na,nb
194 real(sp), intent(IN ) :: za(na),zb(nb),z0
195 real(sp), intent(OUT) :: a(na),b(nb)
196 !-----------------------------------------------------------------------------
197 integer(spi):: na1,nab,i
198 real(sp), dimension(na+nb,na+nb):: w
199 real(sp), dimension(na) :: za0,pa
200 real(sp), dimension(nb) :: zb0,pb
201 real(sp), dimension(na+nb) :: ab
202 !=============================================================================
203 na1=na+1; nab=na+nb
204 za0=za-z0; zb0=zb-z0
205 pa=u1; pb=-u1
206 w=u0; ab=u0
207 w(1,1:na)=u1; ab(1)=u1
208 do i=3,nab; w(i,1:na) =pa*(i-2); pa=pa*za0; enddo
209 do i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; enddo
210 call inv(w,ab)
211 a=ab(1:na); b=ab(na1:nab)
212 end subroutine dfco
213 
224 subroutine ddfco(na,nb,za,zb,z0,a,b) ! Real(dp) version of [DFCO]
225 use pietc, only: u0,u1
226 use pmat, only: inv
227 implicit none
228 integer(spi),intent(in) :: na,nb
229 real(dp), intent(in) :: za(na),zb(nb),z0
230 real(dp), intent(out):: a(na),b(nb)
231 !-----------------------------------------------------------------------------
232 integer(spi) :: na1,nab,i
233 real(dp), dimension(na+nb,na+nb):: w
234 real(dp), dimension(na) :: za0,pa
235 real(dp), dimension(nb) :: zb0,pb
236 real(dp), dimension(na+nb) :: ab
237 !=============================================================================
238 na1=na+1; nab=na+nb
239 za0=za-z0; zb0=zb-z0
240 pa=u1; pb=-u1
241 w=u0; ab=u0
242 w(1,1:na)=u1; ab(1)=u1
243 do i=3,nab; w(i,1:na) =pa*(i-2); pa=pa*za0; enddo
244 do i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; enddo
245 call inv(w,ab)
246 a=ab(1:na); b=ab(na1:nab)
247 end subroutine ddfco
248 
257 subroutine tdfco(xa,xb,a,b)! [DFCO]
258 implicit none
259 real(dp),dimension(:),intent(IN ):: xa,xb
260 real(dp),dimension(:),intent(OUT):: a,b
261 !-----------------------------------------------------------------------------
262 integer(spi):: na,nb
263 !=============================================================================
264 na=size(xa); if(na /= size(a))stop 'In tdfco; sizes of a and xa different'
265 nb=size(xb); if(nb /= size(b))stop 'In tdfco; sizes of b and xb different'
266 call ddfco(na,nb,xa,xb,zero,a,b)
267 end subroutine tdfco
268 
282 subroutine dfco2(na,nb,za,zb,z0,a,b)! [DFCO2]
283 use pietc_s, only: u0,u1
284 use pmat, only: inv
285 implicit none
286 integer(spi), intent(IN ):: na,nb
287 real(sp), intent(IN ):: za(na),zb(nb),z0
288 real(sp), intent(OUT):: a(na),b(nb)
289 !-----------------------------------------------------------------------------
290 integer(spi) :: na1,nab,i
291 real(sp), dimension(na+nb,na+nb):: w
292 real(sp), dimension(na) :: za0,pa
293 real(sp), dimension(nb) :: zb0,pb
294 real(sp), dimension(na+nb) :: ab
295 !=============================================================================
296 na1=na+1; nab=na+nb
297 za0=za-z0; zb0=zb-z0
298 pa=u1; pb=-u1
299 w=u0; ab=u0
300 w(1,1:na)=u1; ab(1)=u1
301 do i=4,nab; w(i,1:na) =pa*(i-2)*(i-3); pa=pa*za0; enddo
302 do i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; enddo
303 call inv(w,ab)
304 a=ab(1:na); b=ab(na1:nab)
305 end subroutine dfco2
306 
317 subroutine ddfco2(na,nb,za,zb,z0,a,b) ! Real(dp) version of [DFCO2]
318 use pietc, only: u0,u1
319 use pmat, only: inv
320 implicit none
321 integer(spi),intent(IN ) :: na,nb
322 real(dp), intent(IN ) :: za(na),zb(nb),z0
323 real(dp), intent(OUT) :: a(na),b(nb)
324 !-----------------------------------------------------------------------------
325 integer(spi) :: na1,nab,i
326 real(dp), dimension(na+nb,na+nb):: w
327 real(dp), dimension(na) :: za0,pa
328 real(dp), dimension(nb) :: zb0,pb
329 real(dp), dimension(na+nb) :: ab
330 !=============================================================================
331 na1=na+1; nab=na+nb
332 za0=za-z0; zb0=zb-z0
333 pa=u1; pb=-u1
334 w=u0; ab=u0
335 w(1,1:na)=u1; ab(1)=u1
336 do i=4,nab; w(i,1:na) =pa*(i-2)*(i-3); pa=pa*za0; enddo
337 do i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; enddo
338 call inv(w,ab)
339 a=ab(1:na); b=ab(na1:nab)
340 end subroutine ddfco2
341 
349 subroutine tdfco2(xa,xb,a,b)! [DFCO2]
350 !=============================================================================
351 real(dp),dimension(:),intent(IN ):: xa,xb
352 real(dp),dimension(:),intent(OUT):: a,b
353 !-----------------------------------------------------------------------------
354 integer(spi):: na,nb
355 !=============================================================================
356 na=size(xa); if(na /= size(a))stop 'In tdfco2; sizes of a and xa different'
357 nb=size(xb); if(nb /= size(b))stop 'In tdfco2; sizes of b and xb different'
358 call ddfco2(na,nb,xa,xb,zero,a,b)
359 end subroutine tdfco2
360 
369 pure subroutine clib(m1,m2,mah1,mah2,a)! [CLIPB]
370 use pietc_s, only: u0
371 implicit none
372 integer(spi), intent(IN ) :: m1, m2, mah1, mah2
373 real(sp), intent(INOUT) :: a(m1,-mah1:mah2)
374 integer(spi):: j
375 do j=1,mah1; a(1:min(m1,j),-j)=u0; enddo
376 do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=u0; enddo
377 end subroutine clib
378 
387 pure subroutine clib_d(m1,m2,mah1,mah2,a)! [CLIPB]
388 use pietc, only: u0
389 implicit none
390 integer(spi),intent(IN ) :: m1, m2, mah1, mah2
391 real(dp), intent(INOUT) :: a(m1,-mah1:mah2)
392 integer(spi):: j
393 do j=1,mah1; a(1:min(m1,j),-j)=u0; enddo
394 do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=u0; enddo
395 end subroutine clib_d
396 
405 pure subroutine clib_c(m1,m2,mah1,mah2,a)! [CLIPB]
406 use pietc, only: c0
407 implicit none
408 integer(spi), intent(IN ) :: m1, m2, mah1, mah2
409 complex(dpc), intent(INOUT) :: a(m1,-mah1:mah2)
410 integer(spi):: j
411 do j=1,mah1; a(1:min(m1,j),-j)=c0; enddo
412 do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=c0; enddo
413 end subroutine clib_c
414 
425 subroutine cad1b(m1,mah1,mah2,mirror2,a)! [CAD1B]
426 use pietc_s, only: u0
427 implicit none
428 integer(spi),intent(IN ):: m1,mah1,mah2,mirror2
429 real(sp), intent(INOUT):: a(0:m1-1,-mah1:mah2)
430 !-----------------------------------------------------------------------------
431 integer(spi):: i,i2,jm,jp,jpmax
432 !=============================================================================
433 if(mirror2+mah1 > mah2)stop 'In CAD1B; mah2 insufficient'
434 do i=0,m1-1; i2=i*2; jpmax=mirror2+mah1-i2; if(jpmax <= -mah1)exit
435  do jm=-mah1,mah2; jp=mirror2-jm-i2; if(jp <= jm)exit
436  a(i,jp)=a(i,jp)+a(i,jm) ! Reflect and add
437  a(i,jm)=u0 ! zero the exterior part
438  enddo
439 enddo
440 end subroutine cad1b
441 
450 subroutine csb1b(m1,mah1,mah2,mirror2,a)! [CSB1B]
451 use pietc_s, only: u0
452 implicit none
453 integer(spi),intent(IN ):: m1,mah1,mah2,mirror2
454 real(sp), intent(INOUT):: a(0:m1-1,-mah1:mah2)
455 !-----------------------------------------------------------------------------
456 integer(spi):: i,i2,jm,jp,jpmax
457 !=============================================================================
458 if(mirror2+mah1 > mah2)stop 'In CSB1B; mah2 insufficient'
459 do i=0,m1-1; i2=i*2; jpmax=mirror2+mah1-i2; if(jpmax < -mah1)exit
460  do jm=-mah1,mah2; jp=mirror2-jm-i2; if(jp < jm)exit
461  a(i,jp)=a(i,jp)-a(i,jm) ! Reflect and subtract
462  a(i,jm)=u0 ! zero the exterior part
463  enddo
464 enddo
465 end subroutine csb1b
466 
476 subroutine cad2b(m1,m2,mah1,mah2,mirror2,a)! [CAD2B]
477 use pietc_s, only: u0
478 implicit none
479 integer(spi),intent(IN ):: m1,m2,mah1,mah2,mirror2
480 real(sp), intent(INOUT):: a(1-m1:0,m1-m2-mah1:m1-m2+mah2)
481 !-----------------------------------------------------------------------------
482 integer(spi):: i,i2,jm,jp,jmmin,nah1,nah2
483 !=============================================================================
484 nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
485 if(mirror2-nah1 > -nah2)stop 'In CAD2B; mah1 insufficient'
486 do i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; if(jmmin >= nah2)exit
487  do jp=nah2,nah1,-1; jm=mirror2-jp-i2; if(jm >= jp)exit
488  a(i,jm)=a(i,jm)+a(i,jp) ! Reflect and add
489  a(i,jp)=u0 ! zero the exterior part
490  enddo
491 enddo
492 end subroutine cad2b
493 
503 subroutine csb2b(m1,m2,mah1,mah2,mirror2,a)! [CSB2B]
504 use pietc_s, only: u0
505 implicit none
506 integer(spi),intent(IN ):: m1,m2,mah1,mah2,mirror2
507 real(sp), intent(INOUT):: a(1-m1:0,m1-m2-mah1:m1-m2+mah2)
508 !-----------------------------------------------------------------------------
509 integer(spi):: i,i2,jm,jp,jmmin,nah1,nah2
510 !=============================================================================
511 nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
512 if(mirror2-nah1 > -nah2)stop 'In CSB2B; mah1 insufficient'
513 do i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; if(jmmin > nah2)exit
514  do jp=nah2,nah1,-1; jm=mirror2-jp-i2; if(jm > jp)exit
515  a(i,jm)=a(i,jm)-a(i,jp) ! Reflect and subtract
516  a(i,jp)=u0 ! zero the exterior part
517  enddo
518 enddo
519 end subroutine csb2b
520 
531 subroutine ldub(m,mah1,mah2,a)! [LDUB]
532 use pietc_s, only: u0,u1
533 implicit none
534 integer(spi),intent(IN ):: m,mah1, mah2
535 real(sp), intent(INOUT):: a(m,-mah1:mah2)
536 !-----------------------------------------------------------------------------
537 integer(spi):: j, imost, jmost, jp, i
538 real(sp) :: ajj, ajji, aij
539 !=============================================================================
540 do j=1,m
541  imost=min(m,j+mah1)
542  jmost=min(m,j+mah2)
543  jp=j+1
544  ajj=a(j,0)
545  if(ajj == u0)then
546  print '(" Failure in LDUB:"/" Matrix requires pivoting or is singular")'
547  stop
548  endif
549  ajji=u1/ajj
550  a(j,0)=ajji
551  do i=jp,imost
552  aij=ajji*a(i,j-i)
553  a(i,j-i)=aij
554  a(i,jp-i:jmost-i)=a(i,jp-i:jmost-i)-aij*a(j,1:jmost-j)
555  enddo
556  a(j,1:jmost-j)=ajji*a(j,1:jmost-j)
557 enddo
558 end subroutine ldub
559 
570 subroutine dldub(m,mah1,mah2,a) ! Real(dp) version of [LDUB]
571 use pietc, only: u0,u1
572 implicit none
573 integer(spi),intent(IN ):: m,mah1, mah2
574 real(dp), intent(INOUT):: a(m,-mah1:mah2)
575 !-----------------------------------------------------------------------------
576 integer(spi):: j, imost, jmost, jp, i
577 real(dp) :: ajj, ajji, aij
578 !=============================================================================
579 do j=1,m
580  imost=min(m,j+mah1)
581  jmost=min(m,j+mah2)
582  jp=j+1
583  ajj=a(j,0)
584  if(ajj == u0)then
585  print '(" Fails in LDUB_d:"/" Matrix requires pivoting or is singular")'
586  stop
587  endif
588  ajji=u1/ajj
589  a(j,0)=ajji
590  do i=jp,imost
591  aij=ajji*a(i,j-i)
592  a(i,j-i)=aij
593  a(i,jp-i:jmost-i)=a(i,jp-i:jmost-i)-aij*a(j,1:jmost-j)
594  enddo
595  a(j,1:jmost-j)=ajji*a(j,1:jmost-j)
596 enddo
597 end subroutine dldub
598 
608 subroutine ldltb(m,mah1,a) ! Real(sp) version of [LDLTB]
609 use pietc_s, only: u0,u1
610 integer(spi),intent(IN ):: m,mah1
611 real(sp), intent(INOUT):: a(m,-mah1:0)
612 !-----------------------------------------------------------------------------
613 integer(spi):: j, imost, jp, i,k
614 real(sp) :: ajj, ajji, aij
615 !=============================================================================
616 do j=1,m
617  imost=min(m,j+mah1)
618  jp=j+1
619  ajj=a(j,0)
620  if(ajj == u0)then
621  print '(" Fails in LDLTB:"/" Matrix requires pivoting or is singular")'
622  stop
623  endif
624  ajji=u1/ajj
625  a(j,0)=ajji
626  do i=jp,imost
627  aij=a(i,j-i)
628  a(i,j-i)=ajji*aij
629  do k=jp,i
630  a(i,k-i)=a(i,k-i)-aij*a(k,j-k)
631  enddo
632  enddo
633 enddo
634 end subroutine ldltb
635 
644 subroutine dldltb(m,mah1,a) ! Real(dp) version of [LDLTB]
645 use pietc, only: u0,u1
646 integer(spi),intent(IN ) :: m,mah1
647 real(dp), intent(INOUT) :: a(m,-mah1:0)
648 !-----------------------------------------------------------------------------
649 integer(spi):: j, imost, jp, i,k
650 real(dp) :: ajj, ajji, aij
651 !=============================================================================
652 do j=1,m
653  imost=min(m,j+mah1)
654  jp=j+1
655  ajj=a(j,0)
656  if(ajj == u0)then
657  print '(" Fails in LDLTB_d:"/" Matrix requires pivoting or is singular")'
658  stop
659  endif
660  ajji=u1/ajj
661  a(j,0)=ajji
662  do i=jp,imost
663  aij=a(i,j-i)
664  a(i,j-i)=ajji*aij
665  do k=jp,i
666  a(i,k-i)=a(i,k-i)-aij*a(k,j-k)
667  enddo
668  enddo
669 enddo
670 end subroutine dldltb
671 
682 subroutine udlb(m,mah1,mah2,a) ! Reversed-index version of ldub [UDLB]
683 implicit none
684 integer(spi), intent(IN ) :: m,mah1,mah2
685 real(sp),dimension(m,-mah1:mah2),intent(INOUT) :: a(m,-mah1:mah2)
686 !-----------------------------------------------------------------------------
687 real(sp),dimension(m,-mah2:mah1):: at
688 !=============================================================================
689 at=a(m:1:-1,mah2:-mah1:-1); call ldub(m,mah2,mah1,at)
690 a=at(m:1:-1,mah1:-mah2:-1)
691 end subroutine udlb
692 
703 subroutine dudlb(m,mah1,mah2,a) ! real(dp) version of udlb [UDLB]
704 implicit none
705 integer(spi), intent(IN ) :: m,mah1,mah2
706 real(dp),dimension(m,-mah1:mah2),intent(INOUT) :: a(m,-mah1:mah2)
707 !-----------------------------------------------------------------------------
708 real(dp),dimension(m,-mah2:mah1):: at
709 !=============================================================================
710 at=a(m:1:-1,mah2:-mah1:-1); call dldub(m,mah2,mah1,at)
711 a=at(m:1:-1,mah1:-mah2:-1)
712 end subroutine dudlb
713 
732 subroutine l1ubb(m,mah1,mah2,mbh1,mbh2,a,b)! [L1UBB]
733 use pietc_s, only: u0,u1
734 implicit none
735 integer(spi), intent(IN ) :: m,mah1, mah2, mbh1, mbh2
736 real(sp), intent(INOUT) :: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
737 !-----------------------------------------------------------------------------
738 integer(spi):: j, imost, jmost, jleast, jp, i
739 real(sp) :: ajj, ajji, aij
740 !=============================================================================
741 do j=1,m
742  imost=min(m,j+mah1)
743  jmost=min(m,j+mah2)
744  jleast=max(1,j-mah1)
745  jp=j+1
746  ajj=a(j,0)
747  if(ajj == u0)stop 'In L1UBB; zero element found in diagonal factor'
748  ajji=u1/ajj
749  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
750  do i=jp,imost
751  aij=a(i,j-i)
752  a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
753  enddo
754  a(j,0)=u1
755  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
756 enddo
757 end subroutine l1ubb
758 
770 subroutine dl1ubb(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of [L1UBB]
771 use pietc, only: u0,u1
772 implicit none
773 integer(spi),intent(IN ) :: m,mah1, mah2, mbh1, mbh2
774 real(dp), intent(INOUT) :: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
775 !-----------------------------------------------------------------------------
776 integer(spi):: j, imost, jmost, jleast, jp, i
777 real(dp) :: ajj, ajji, aij
778 !=============================================================================
779 do j=1,m
780  imost=min(m,j+mah1)
781  jmost=min(m,j+mah2)
782  jleast=max(1,j-mah1)
783  jp=j+1
784  ajj=a(j,0)
785  if(ajj == u0)stop 'In L1UBB_d; zero element found in diagonal factor'
786  ajji=u1/ajj
787  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
788  do i=jp,imost
789  aij=a(i,j-i)
790  a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
791  enddo
792  a(j,0)=u1
793  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
794 enddo
795 end subroutine dl1ubb
796 
817 subroutine l1ueb(m,mah1,mah2,mbh1,mbh2,a,b)! [L1UEB]
818 use pietc_s, only: u0,u1
819 implicit none
820 integer(spi),intent(IN ) :: m,mah1, mah2, mbh1, mbh2
821 real(sp), intent(INOUT) :: a(0:m,-mah1:mah2), b(m,-mbh1:mbh2)
822 !-----------------------------------------------------------------------------
823 integer(spi):: j, imost, jmost, jleast, jp, i
824 real(sp) :: ajj, ajji, aij
825 !=============================================================================
826 do j=1,m
827  imost=min(m,j+mah1)
828  jmost=min(m,j+mah2)
829  jleast=max(0,j-mah1)
830  jp=j+1
831  ajj=a(j,0)
832  if(ajj == u0)stop 'In L1UEB; zero element found in diagonal factor'
833  ajji=u1/ajj
834  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
835  do i=jp,imost
836  aij=a(i,j-i)
837  a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
838  enddo
839  a(j,0)=u1
840  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
841 enddo
842 end subroutine l1ueb
843 
854 subroutine dl1ueb(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of [L1UEB]
855 use pietc, only: u0,u1
856 implicit none
857 integer(spi),intent(IN ):: m,mah1, mah2, mbh1, mbh2
858 real(dp), intent(INOUT):: a(0:,-mah1:), b(:,-mbh1:)
859 !-----------------------------------------------------------------------------
860 integer(spi):: j, imost, jmost, jleast, jp, i
861 real(dp) :: ajj, ajji, aij
862 !=============================================================================
863 do j=1,m
864  imost=min(m,j+mah1)
865  jmost=min(m,j+mah2)
866  jleast=max(0,j-mah1)
867  jp=j+1
868  ajj=a(j,0)
869  if(ajj == u0)stop 'In L1UEB_D; zero element found in diagonal factor'
870  ajji=u1/ajj
871  a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
872  do i=jp,imost
873  aij=a(i,j-i)
874  a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
875  enddo
876  a(j,0)=u1
877  b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
878 enddo
879 end subroutine dl1ueb
880 
891 subroutine udlbv(m,mah1,mah2,a,v)! [UDLBV]
892 implicit none
893 integer(spi),intent(IN ):: m, mah1, mah2
894 real(sp), intent(IN ):: a(m,-mah1:mah2)
895 real(sp), intent(INOUT):: v(m)
896 !-----------------------------------------------------------------------------
897 integer(spi):: i, j
898 real(sp) :: vj
899 !=============================================================================
900 do j=1,m
901  vj=v(j)
902  do i=j+1,min(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; enddo; v(j)=a(j,0)*vj
903 enddo
904 do j=m,2,-1
905  vj=v(j)
906  do i=max(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; enddo
907 enddo
908 end subroutine udlbv
909 
920 subroutine dudlbv(m,mah1,mah2,a,v)! [udlbv]
921 implicit none
922 integer(spi),intent(IN ) :: m, mah1, mah2
923 real(dp), intent(IN ) :: a(m,-mah1:mah2)
924 real(dp), intent(INOUT) :: v(m)
925 !-----------------------------------------------------------------------------
926 integer(spi):: i, j
927 real(dp) :: vj
928 !=============================================================================
929 do j=1,m
930  vj=v(j)
931  do i=j+1,min(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; enddo; v(j)=a(j,0)*vj
932 enddo
933 do j=m,2,-1
934  vj=v(j)
935  do i=max(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; enddo
936 enddo
937 end subroutine dudlbv
938 
949 subroutine ltdlbv(m,mah1,a,v)! [ltdlbv]
950 implicit none
951 integer(spi),intent(IN ) :: m, mah1
952 real(sp), intent(IN ) :: a(m,-mah1:0)
953 real(sp), intent(INOUT) :: v(m)
954 !-----------------------------------------------------------------------------
955 integer(spi):: i, j
956 real(sp) :: vj
957 !=============================================================================
958 do j=1,m
959  vj=v(j)
960  do i=j+1,min(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; enddo; v(j)=a(j,0)*vj
961 enddo
962 do j=m,2,-1
963  vj=v(j)
964  do i=max(1,j-mah1),j-1; v(i)=v(i)-a(j,i-j)*vj; enddo
965 enddo
966 end subroutine ltdlbv
967 
968 
979 subroutine dltdlbv(m,mah1,a,v)! [ltdlbv]
980 implicit none
981 integer(spi),intent(IN ) :: m, mah1
982 real(dp), intent(IN ) :: a(m,-mah1:0)
983 real(dp), intent(INOUT) :: v(m)
984 !-----------------------------------------------------------------------------
985 integer(spi):: i, j
986 real(dp) :: vj
987 !=============================================================================
988 do j=1,m
989  vj=v(j)
990  do i=j+1,min(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; enddo; v(j)=a(j,0)*vj
991 enddo
992 do j=m,2,-1
993  vj=v(j)
994  do i=max(1,j-mah1),j-1; v(i)=v(i)-a(j,i-j)*vj; enddo
995 enddo
996 end subroutine dltdlbv
997 
1010 subroutine udlbx(mx,mah1,mah2,my,a,v)! [UDLBX]
1011 implicit none
1012 integer(spi),intent(IN ) :: mx, mah1, mah2, my
1013 real(sp), intent(IN ) :: a(mx,-mah1:mah2)
1014 real(sp), intent(INOUT) :: v(mx,my)
1015 !-----------------------------------------------------------------------------
1016 integer(spi):: jx, ix
1017 !=============================================================================
1018 do jx=1,mx
1019  do ix=jx+1,min(mx,jx+mah1); v(ix,:) = v(ix,:) - a(ix,jx-ix)*v(jx,:); enddo
1020  v(jx,:) = a(jx,0) * v(jx,:)
1021 enddo
1022 do jx=mx,2,-1
1023  do ix=max(1,jx-mah2),jx-1; v(ix,:) = v(ix,:) - a(ix,jx-ix)*v(jx,:); enddo
1024 enddo
1025 end subroutine udlbx
1026 
1039 subroutine udlby(my,mah1,mah2,mx,a,v)! [UDLBY]
1040 implicit none
1041 integer(spi),intent(IN ) :: my, mah1, mah2, mx
1042 real(sp), intent(IN ) :: a(my,-mah1:mah2)
1043 real(sp), intent(INOUT) :: v(mx,my)
1044 !-----------------------------------------------------------------------------
1045 integer(spi):: iy, jy
1046 !=============================================================================
1047 do jy=1,my
1048  do iy=jy+1,min(my,jy+mah1); v(:,iy) = v(:,iy)-a(iy,jy-iy)*v(:,jy); enddo
1049  v(:,jy)=a(jy,0)*v(:,jy)
1050 enddo
1051 do jy=my,2,-1
1052  do iy=max(1,jy-mah2),jy-1; v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); enddo
1053 enddo
1054 end subroutine udlby
1055 
1066 subroutine udlvb(m,mah1,mah2,v,a)! [UDLVB]
1067 implicit none
1068 integer(spi), intent(IN ) :: m, mah1, mah2
1069 real(sp), intent(IN ) :: a(m,-mah1:mah2)
1070 real(sp), intent(INOUT) :: v(m)
1071 !-----------------------------------------------------------------------------
1072 integer(spi):: i, j
1073 real(sp) :: vi
1074 !=============================================================================
1075 do i=1,m
1076  vi=v(i)
1077  do j=i+1,min(m,i+mah2); v(j)=v(j)-vi*a(i,j-i); enddo
1078  v(i)=vi*a(i,0)
1079 enddo
1080 do i=m,2,-1
1081  vi=v(i)
1082  do j=max(1,i-mah1),i-1; v(j)=v(j)-vi*a(i,j-i); enddo
1083 enddo
1084 end subroutine udlvb
1085 
1098 subroutine udlxb(mx,mah1,mah2,my,v,a)! [UDLXB]
1099 implicit none
1100 integer(spi),intent(IN ) :: mx, mah1, mah2, my
1101 real(sp), intent(IN ) :: a(mx,-mah1:mah2)
1102 real(sp), intent(INOUT) :: v(mx,my)
1103 !-----------------------------------------------------------------------------
1104 integer(spi):: ix, jx
1105 !=============================================================================
1106 do ix=1,mx
1107  do jx=ix+1,min(mx,ix+mah2); v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); enddo
1108  v(ix,:)=v(ix,:)*a(ix,0)
1109 enddo
1110 do ix=mx,2,-1
1111  do jx=max(1,ix-mah1),ix-1; v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); enddo
1112 enddo
1113 end subroutine udlxb
1114 
1127 subroutine udlyb(my,mah1,mah2,mx,v,a)! [UDLYB]
1128 implicit none
1129 integer(spi),intent(IN ) :: my, mah1, mah2, mx
1130 real(sp), intent(IN ) :: a(my,-mah1:mah2)
1131 real(sp), intent(INOUT) :: v(mx,my)
1132 !-----------------------------------------------------------------------------
1133 integer(spi):: iy, jy
1134 !=============================================================================
1135 do iy=1,my
1136  do jy=iy+1,min(my,iy+mah2); v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); enddo
1137  v(:,iy)=v(:,iy)*a(iy,0)
1138 enddo
1139 do iy=my,2,-1
1140  do jy=max(1,iy-mah1),iy-1; v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); enddo
1141 enddo
1142 end subroutine udlyb
1143 
1154 subroutine u1lbv(m,mah1,mah2,a,v)! [U1LBV]
1155 implicit none
1156 integer(spi),intent(IN ) :: m, mah1, mah2
1157 real(sp), intent(IN ) :: a(m,-mah1:mah2)
1158 real(sp), intent(INOUT) :: v(m)
1159 !-----------------------------------------------------------------------------
1160 integer(spi):: i, j
1161 real(sp) :: vj
1162 !=============================================================================
1163 do j=1,m
1164  vj=v(j)
1165  do i=j+1,min(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; enddo
1166 enddo
1167 do j=m,2,-1
1168  vj=v(j)
1169  do i=max(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; enddo
1170 enddo
1171 end subroutine u1lbv
1172 
1185 subroutine u1lbx(mx,mah1,mah2,my,a,v)! [U1LBX]
1186 implicit none
1187 integer(spi),intent(IN ) :: mx, mah1, mah2, my
1188 real(sp), intent(IN ) :: a(mx,-mah1:mah2)
1189 real(sp), intent(INOUT) :: v(mx,my)
1190 !-----------------------------------------------------------------------------
1191 integer(spi):: ix, jx
1192 !=============================================================================
1193 do jx=1,mx
1194  do ix=jx+1,min(mx,jx+mah1); v(ix,:)=v(ix,:)-a(ix,jx-ix)*v(jx,:); enddo
1195 enddo
1196 do jx=mx,2,-1
1197  do ix=max(1,jx-mah2),jx-1; v(ix,:)=v(ix,:)-a(ix,jx-ix)*v(jx,:); enddo
1198 enddo
1199 end subroutine u1lbx
1200 
1213 subroutine u1lby(my,mah1,mah2,mx,a,v)! [U1LBY]
1214 implicit none
1215 integer(spi),intent(IN ) :: my, mah1, mah2, mx
1216 real(sp), intent(IN ) :: a(my,-mah1:mah2)
1217 real(sp), intent(INOUT) :: v(mx,my)
1218 !-----------------------------------------------------------------------------
1219 integer(spi):: iy, jy
1220 !=============================================================================
1221 do jy=1,my
1222  do iy=jy+1,min(my,jy+mah1); v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); enddo
1223 enddo
1224 do jy=my,2,-1
1225  do iy=max(1,jy-mah2),jy-1; v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); enddo
1226 enddo
1227 end subroutine u1lby
1228 
1239 subroutine u1lvb(m,mah1,mah2,v,a)! [U1LVB]
1240 implicit none
1241 integer(spi),intent(IN ) :: m, mah1, mah2
1242 real(sp), intent(IN ) :: a(m,-mah1:mah2)
1243 real(sp), intent(INOUT) :: v(m)
1244 !-----------------------------------------------------------------------------
1245 integer(spi):: i, j
1246 real(sp) :: vi
1247 !=============================================================================
1248 do i=1,m
1249  vi=v(i)
1250  do j=i+1,min(m,i+mah2); v(j)=v(j)-vi*a(i,j-i); enddo
1251 enddo
1252 do i=m,2,-1
1253  vi=v(i)
1254  do j=max(1,i-mah1),i-1; v(j)=v(j)-vi*a(i,j-i); enddo
1255 enddo
1256 end subroutine u1lvb
1257 
1270 subroutine u1lxb(mx,mah1,mah2,my,v,a)! [U1LXB]
1271 implicit none
1272 integer(spi),intent(IN ) :: mx, mah1, mah2, my
1273 real(sp), intent(IN ) :: a(mx,-mah1:mah2)
1274 real(sp), intent(INOUT) :: v(mx,my)
1275 !-----------------------------------------------------------------------------
1276 integer(spi):: ix, jx
1277 !=============================================================================
1278 do ix=1,mx
1279  do jx=ix+1,min(mx,ix+mah2); v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); enddo
1280 enddo
1281 do ix=mx,2,-1
1282  do jx=max(1,ix-mah1),ix-1; v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); enddo
1283 enddo
1284 end subroutine u1lxb
1285 
1298 subroutine u1lyb(my,mah1,mah2,mx,v,a)! [U1LYB]
1299 implicit none
1300 integer(spi),intent(IN ) :: my, mah1, mah2, mx
1301 real(sp), intent(IN ) :: a(my,-mah1:mah2)
1302 real(sp), intent(INOUT) :: v(mx,my)
1303 !-----------------------------------------------------------------------------
1304 integer(spi):: iy, jy
1305 !=============================================================================
1306 do iy=1,my
1307  do jy=iy+1,min(my,iy+mah2); v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); enddo
1308 enddo
1309 do iy=my,2,-1
1310  do jy=max(1,iy-mah1),iy-1; v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); enddo
1311 enddo
1312 end subroutine u1lyb
1313 
1322 subroutine linbv(m,mah1,mah2,a,v)! [LINBV]
1323 implicit none
1324 integer(spi),intent(IN ) :: m, mah1, mah2
1325 real(sp), intent(INOUT) :: a(m,-mah1:mah2), v(m)
1326 !=============================================================================
1327 call ldub(m,mah1,mah2,a)
1328 call udlbv(m,mah1,mah2,a,v)
1329 end subroutine linbv
1330 
1340 subroutine wrtb(m1,m2,mah1,mah2,a)! [WRTB]
1341 implicit none
1342 integer(spi),intent(IN) :: m1, m2, mah1, mah2
1343 real(sp), intent(IN) :: a(m1,-mah1:mah2)
1344 !-----------------------------------------------------------------------------
1345 integer(spi):: i1, i2, i, j1, j2, j, nj1
1346 !=============================================================================
1347 do i1=1,m1,20
1348  i2=min(i1+19,m1)
1349  print '(7x,6(i2,10x))',(j,j=-mah1,mah2)
1350  do i=i1,i2
1351  j1=max(-mah1,1-i)
1352  j2=min(mah2,m2-i)
1353  nj1=j1+mah1
1354  if(nj1==0)print '(1x,i3,6(1x,e12.5))', i,(a(i,j),j=j1,j2)
1355  if(nj1==1)print '(1x,i3,12x,5(1x,e12.5))',i,(a(i,j),j=j1,j2)
1356  if(nj1==2)print '(1x,i3,24x,4(1x,e12.5))',i,(a(i,j),j=j1,j2)
1357  if(nj1==3)print '(1x,i3,36x,3(1x,e12.5))',i,(a(i,j),j=j1,j2)
1358  if(nj1==4)print '(1x,i3,48x,2(1x,e12.5))',i,(a(i,j),j=j1,j2)
1359  if(nj1==5)print '(1x,i3,60x,1(1x,e12.5))',i,(a(i,j),j=j1,j2)
1360  enddo
1361  read(*,*)
1362 enddo
1363 end subroutine wrtb
1364 
1365 end module pmat2
pkind::sp
integer, parameter sp
Single precision real kind.
Definition: pkind.f90:10
pmat2::dldub
subroutine dldub(m, mah1, mah2, a)
[L]*[D]*[U] factoring of double precision band-matrix.
Definition: pmat2.f90:571
pkind::dp
integer, parameter dp
Double precision real kind.
Definition: pkind.f90:11
pmat2::dudlbv
subroutine dudlbv(m, mah1, mah2, a, v)
Back-substitution step of linear inversion involving banded LDU factored matrix and input vector...
Definition: pmat2.f90:921
pkind
Standard integer, real, and complex single and double precision kinds.
Definition: pkind.f90:7
pietc::u0
real(dp), parameter u0
zero
Definition: pietc.f90:19
pmat2::clib_d
pure subroutine clib_d(m1, m2, mah1, mah2, a)
Clip (set to zero) the unused values in a banded matrix representation.
Definition: pmat2.f90:388
pmat2::zero
real(dp), parameter zero
Double precision real zero.
Definition: pmat2.f90:54
pmat2::dltdlbv
subroutine dltdlbv(m, mah1, a, v)
Like udlbv, except assuming a is the ltdl decomposition of a SYMMETRIC banded matrix, with only the non-upper part provided (to avoid redundancy)
Definition: pmat2.f90:980
pmat2::davco
subroutine davco(na, nb, za, zb, z0, a, b)
Double precision version of subroutine avco for midpoint interpolation.
Definition: pmat2.f90:134
pietc::c0
complex(dpc), parameter c0
complex zero
Definition: pietc.f90:112
pmat2::tdfco2
subroutine tdfco2(xa, xb, a, b)
Simplified computation of compact 2nd-derivative coefficients.
Definition: pmat2.f90:350
pmat2::dl1ubb
subroutine dl1ubb(m, mah1, mah2, mbh1, mbh2, a, b)
Double precision version of L1UBB.
Definition: pmat2.f90:771
pmat2::ddfco2
subroutine ddfco2(na, nb, za, zb, z0, a, b)
Double precision version of DFCO2 to get 2nd-derivative coefficients.
Definition: pmat2.f90:318
pmat2::tavco
subroutine tavco(xa, xb, a, b)
Simplified computation of compact midpoint interpolation coefficients.
Definition: pmat2.f90:165
pmat2::tdfco
subroutine tdfco(xa, xb, a, b)
Simplified computation of compact differencing coefficients to get derivatives d from cumulatives c...
Definition: pmat2.f90:258
pmat2::dudlb
subroutine dudlb(m, mah1, mah2, a)
[U]*[D]*[L] factoring of double precision band matrix A [U] is upper triangular with unit main diagon...
Definition: pmat2.f90:704
pmat2::ddfco
subroutine ddfco(na, nb, za, zb, z0, a, b)
Double precision version of dfco for compact differentiation coefficients.
Definition: pmat2.f90:225
pmat2::clib
pure subroutine clib(m1, m2, mah1, mah2, a)
Clip (set to zero) the unused values in a banded matrix representation.
Definition: pmat2.f90:370
pietc::u1
real(dp), parameter u1
one
Definition: pietc.f90:20
pmat2::clib_c
pure subroutine clib_c(m1, m2, mah1, mah2, a)
Clip (set to zero) the unused values in a banded matrix representation.
Definition: pmat2.f90:406
pietc
Some of the commonly used constants (pi etc) mainly for double-precision subroutines.
Definition: pietc.f90:14
pmat2::dl1ueb
subroutine dl1ueb(m, mah1, mah2, mbh1, mbh2, a, b)
Double precision version of L1UEB.
Definition: pmat2.f90:855
pkind::spi
integer, parameter spi
Single precision integer kind.
Definition: pkind.f90:8
pkind::dpc
integer, parameter dpc
Double precision real kind.
Definition: pkind.f90:13
pmat2
Routines dealing with the operations of banded matrices.
Definition: pmat2.f90:21
pmat2::dldltb
subroutine dldltb(m, mah1, a)
[L]*[D]*[L^T] factoring of symmetric matrix A (root-free Cholesky).
Definition: pmat2.f90:645