Commit 33505081 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

Fortran code modernized (EQUIVALENCE and ENTRY statements replaced) and checked

git-svn-id: http://svnsrv.desy.de/public/MillepedeII/trunk@220 3547b9b0-65b8-46d3-b95d-921b3f43af62
parent 0b4e005e
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
!! \author Claus Kleinwort, DESY (maintenance and developement) !! \author Claus Kleinwort, DESY (maintenance and developement)
!! !!
!! \copyright !! \copyright
!! Copyright (c) 2009 - 2015 Deutsches Elektronen-Synchroton, !! Copyright (c) 2009 - 2021 Deutsches Elektronen-Synchroton,
!! Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n !! Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
!! This library is free software; you can redistribute it and/or modify !! This library is free software; you can redistribute it and/or modify
!! it under the terms of the GNU Library General Public License as !! it under the terms of the GNU Library General Public License as
...@@ -199,8 +199,8 @@ ...@@ -199,8 +199,8 @@
!! W(.) and V(.) are symmetric n-by-n matrices with (N*N+N)/2 elements !! W(.) and V(.) are symmetric n-by-n matrices with (N*N+N)/2 elements
!! !!
!! SUBROUTINE DCHDEC(W,N, AUX) ! decomposition, symmetric matrix !! SUBROUTINE DCHDEC(W,N, AUX) ! decomposition, symmetric matrix
!! ENTRY DCHSLV(W,N,B, X) ! solution B -> X !! SUBROUTINE DCHSLV(W,N,B, X) ! solution B -> X
!! ENTRY DCHINV(W,N, V) ! inversion !! SUBROUTINE DCHINV(W,N, V) ! inversion
!! !!
!! SUBROUTINE DCFDEC(W,N) ! modified composition, symmetric !! SUBROUTINE DCFDEC(W,N) ! modified composition, symmetric
!! ! alternative to DCHDEC !! ! alternative to DCHDEC
...@@ -211,9 +211,9 @@ ...@@ -211,9 +211,9 @@
!! The symmetric matrix VS has (N*N+N)/2 elements !! The symmetric matrix VS has (N*N+N)/2 elements
!! !!
!! SUBROUTINE DBCDEC(W,MP1,N, AUX) ! decomposition, bandwidth M !! SUBROUTINE DBCDEC(W,MP1,N, AUX) ! decomposition, bandwidth M
!! ENTRY DBCSLV(W,MP1,N,B, X) ! solution B -> X !! SUBROUTINE DBCSLV(W,MP1,N,B, X) ! solution B -> X
!! ENTRY DBCIEL(W,MP1,N, V) ! V = inverse band matrix elements !! SUBROUTINE DBCIEL(W,MP1,N, V) ! V = inverse band matrix elements
!! ENTRY DBCINV(W,MP1,N, VS) ! V = inverse symmetric matrix !! SUBROUTINE DBCINV(W,MP1,N, VS) ! V = inverse symmetric matrix
!! !!
!! SUBROUTINE DBFDEC(W,MP1,N) ! modified decomposition, bandwidth M !! SUBROUTINE DBFDEC(W,MP1,N) ! modified decomposition, bandwidth M
!! ! alternative to DBCDEC !! ! alternative to DBCDEC
...@@ -222,12 +222,12 @@ ...@@ -222,12 +222,12 @@
!! SUBROUTINE DBCPRV(W,MP1,N,B) ! print corr band matrix and pars !! SUBROUTINE DBCPRV(W,MP1,N,B) ! print corr band matrix and pars
!! !!
!! SUBROUTINE DB2DEC(W,N, AUX) ! decomposition (M=1) !! SUBROUTINE DB2DEC(W,N, AUX) ! decomposition (M=1)
!! ENTRY DB2SLV(W,N,B, X) ! solution B -> X !! SUBROUTINE DB2SLV(W,N,B, X) ! solution B -> X
!! ENTRY DB2IEL(W,N, V) ! V = inverse band matrix elements !! SUBROUTINE DB2IEL(W,N, V) ! V = inverse band matrix elements
!! !!
!! SUBROUTINE DB3DEC(W,N, AUX) ! decomposition (M=2) !! SUBROUTINE DB3DEC(W,N, AUX) ! decomposition (M=2)
!! ENTRY DB3SLV(W,N,B, X) ! solution B -> X !! SUBROUTINE DB3SLV(W,N,B, X) ! solution B -> X
!! ENTRY DB3IEL(W,N, V) ! V = inverse band matrix elements !! SUBROUTINE DB3IEL(W,N, V) ! V = inverse band matrix elements
!!\endverbatim !!\endverbatim
...@@ -235,10 +235,10 @@ ...@@ -235,10 +235,10 @@
!> Decomposition of symmetric matrix. !> Decomposition of symmetric matrix.
!! !!
!! ENTRY DCHSLV(W,N,B, X) for solution B -> X \n !! SUBROUTINE DCHSLV(W,N,B, X) for solution B -> X \n
!! ENTRY DCHINV(W,N,V) for inversion !! SUBROUTINE DCHINV(W,N,V) for inversion
!! !!
!! \param [in,out] W symmetirc matrix !! \param [in,out] W symmetric matrix, replaced by decomposition
!! \param [in] N size !! \param [in] N size
!! \param [in] AUX scratch array !! \param [in] AUX scratch array
...@@ -252,15 +252,12 @@ SUBROUTINE dchdec(w,n, aux) ...@@ -252,15 +252,12 @@ SUBROUTINE dchdec(w,n, aux)
INTEGER(mpi) :: jj INTEGER(mpi) :: jj
INTEGER(mpi) :: k INTEGER(mpi) :: k
INTEGER(mpi) :: kk INTEGER(mpi) :: kk
INTEGER(mpi) :: l REAL(mpd) :: ratio
INTEGER(mpi) :: m
REAL(mpd), INTENT(IN OUT) :: w(*) REAL(mpd), INTENT(IN OUT) :: w(*)
INTEGER(mpi), INTENT(IN) :: n INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(OUT) :: aux(n) REAL(mpd), INTENT(OUT) :: aux(n)
REAL(mpd) :: b(*),x(*),v(*),suma,ratio
! ... ! ...
DO i=1,n DO i=1,n
aux(i)=16.0_mpd*w((i*i+i)/2) ! save diagonal elements aux(i)=16.0_mpd*w((i*i+i)/2) ! save diagonal elements
...@@ -285,11 +282,31 @@ SUBROUTINE dchdec(w,n, aux) ...@@ -285,11 +282,31 @@ SUBROUTINE dchdec(w,n, aux)
jj=jj+j jj=jj+j
END DO ! J END DO ! J
END DO ! I END DO ! I
RETURN
END SUBROUTINE dchdec
!> solution B -> X
!!
!! \param [in,out] W (decomposition of) symmetric matrix
!! \param [in] N size
!! \param [in] B r.h.s.
!! \param [out] X solution
RETURN SUBROUTINE dchslv(w,n,b, x)
USE mpdef
implicit none
INTEGER(mpi) :: i
INTEGER(mpi) :: ii
INTEGER(mpi) :: k
INTEGER(mpi) :: kk
REAL(mpd) :: suma
REAL(mpd), INTENT(IN) :: w(*)
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(IN) :: b(n)
REAL(mpd), INTENT(OUT) :: x(n)
ENTRY dchslv(w,n,b, x) ! solution B -> X
WRITE(*,*) 'before copy' WRITE(*,*) 'before copy'
DO i=1,n DO i=1,n
x(i)=b(i) x(i)=b(i)
...@@ -317,8 +334,30 @@ SUBROUTINE dchdec(w,n, aux) ...@@ -317,8 +334,30 @@ SUBROUTINE dchdec(w,n, aux)
END DO END DO
WRITE(*,*) 'after backward' WRITE(*,*) 'after backward'
RETURN RETURN
END SUBROUTINE dchslv
!> inversion
!!
!! \param [in,out] W (decomposition of) symmetric matrix
!! \param [in] N size
!! \param [out] V inverse symmetric matrix
SUBROUTINE dchinv(w,n,v)
USE mpdef
implicit none
INTEGER(mpi) :: i
INTEGER(mpi) :: ii
INTEGER(mpi) :: j
INTEGER(mpi) :: k
INTEGER(mpi) :: l
INTEGER(mpi) :: m
REAL(mpd) :: suma
REAL(mpd), INTENT(IN) :: w(*)
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(OUT) :: v(*)
ENTRY dchinv(w,n,v) ! inversion
ii=(n*n-n)/2 ii=(n*n-n)/2
DO i=n,1,-1 DO i=n,1,-1
suma=w(ii+i) ! (I,I) suma=w(ii+i) ! (I,I)
...@@ -333,7 +372,7 @@ SUBROUTINE dchdec(w,n, aux) ...@@ -333,7 +372,7 @@ SUBROUTINE dchdec(w,n, aux)
END DO END DO
ii=ii-i+1 ii=ii-i+1
END DO END DO
END SUBROUTINE dchdec END SUBROUTINE dchinv
!> Etimate condition. !> Etimate condition.
!! !!
...@@ -406,23 +445,25 @@ END FUNCTION condes ...@@ -406,23 +445,25 @@ END FUNCTION condes
! inverse and band part of the inverse ! inverse and band part of the inverse
!> Decomposition of symmetric band matrix. !> Decomposition of symmetric band matrix.
!! !!
!! ENTRY DBCSLV(W,MP1,N,B, X) for solution B -> X \n !! SUBROUTINE DBCSLV(W,MP1,N,B, X) for solution B -> X \n
!! ENTRY DBCIEL(W,MP1,N, V), V = inverse band matrix elements \n !! SUBROUTINE DBCIEL(W,MP1,N, V), V = inverse band matrix elements \n
!! ENTRY DBCINB(W,MP1,N, VS), VS = band part of inverse symmetric matrix \n !! SUBROUTINE DBCINB(W,MP1,N, VS), VS = band part of inverse symmetric matrix \n
!! ENTRY DBCINV(W,MP1,N, VS), V = inverse symmetric matrix !! SUBROUTINE DBCINV(W,MP1,N, VS), VS = inverse symmetric matrix
!! !!
!! \param [in,out] W symmetric band matrix !! \param [in,out] W symmetric band matrix, replaced by decomposition
!! \param [in] MP1 band width (M) + 1 !! \param [in] MP1 band width (M) + 1
!! \param [in] N size !! \param [in] N size
!! \param [in] AUX scratch array !! \param [in] AUX scratch array
!!
SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M !! decomposition, bandwidth M
SUBROUTINE dbcdec(w,mp1,n, aux)
USE mpdef USE mpdef
implicit none implicit none
INTEGER(mpi) :: i INTEGER(mpi) :: i
INTEGER(mpi) :: j INTEGER(mpi) :: j
INTEGER(mpi) :: k INTEGER(mpi) :: k
REAL(mpd) :: rxw
! M=MP1-1 N*M(M-1) dot operations ! M=MP1-1 N*M(M-1) dot operations
REAL(mpd), INTENT(IN OUT) :: w(mp1,n) REAL(mpd), INTENT(IN OUT) :: w(mp1,n)
...@@ -430,7 +471,6 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M ...@@ -430,7 +471,6 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M
INTEGER(mpi), INTENT(IN) :: n INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(OUT) :: aux(n) REAL(mpd), INTENT(OUT) :: aux(n)
! decompos ! decompos
REAL(mpd) :: v(mp1,n),b(n),x(n), vs(*),rxw
! ... ! ...
DO i=1,n DO i=1,n
aux(i)=16.0_mpd*w(1,i) ! save diagonal elements aux(i)=16.0_mpd*w(1,i) ! save diagonal elements
...@@ -450,8 +490,29 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M ...@@ -450,8 +490,29 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M
END DO END DO
END DO END DO
RETURN RETURN
END SUBROUTINE dbcdec
!> solution B -> X
!!
!! \param [in,out] W (decomposition of) )symmetric band matrix
!! \param [in] MP1 band width (M) + 1
!! \param [in] N size
!! \param [in] B r.h.s.
!! \param [out] X solution
SUBROUTINE dbcslv(w,mp1,n,b, x)
USE mpdef
ENTRY dbcslv(w,mp1,n,b, x) ! solution B -> X implicit none
INTEGER(mpi) :: i
INTEGER(mpi) :: j
REAL(mpd) :: rxw
REAL(mpd), INTENT(IN) :: w(mp1,n)
INTEGER(mpi), INTENT(IN) :: mp1
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(IN) :: b(n)
REAL(mpd), INTENT(OUT) :: x(n)
! N*(2M-1) dot operations ! N*(2M-1) dot operations
DO i=1,n DO i=1,n
x(i)=b(i) x(i)=b(i)
...@@ -469,8 +530,28 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M ...@@ -469,8 +530,28 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M
x(i)=rxw x(i)=rxw
END DO END DO
RETURN RETURN
END SUBROUTINE dbcslv
!> V = inverse band matrix elements
!!
!! \param [in,out] W (decomposition of) )symmetric band matrix
!! \param [in] MP1 band width (M) + 1
!! \param [in] N size
!! \param [out] V inverse band matrix elements
SUBROUTINE dbciel(w,mp1,n,v)
USE mpdef
ENTRY dbciel(w,mp1,n, v) ! V = inverse band matrix elements implicit none
INTEGER(mpi) :: i
INTEGER(mpi) :: j
INTEGER(mpi) :: k
REAL(mpd) :: rxw
REAL(mpd), INTENT(IN) :: w(mp1,n)
INTEGER(mpi), INTENT(IN) :: mp1
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(OUT) :: v(mp1,n)
! N*M*(M-1) dot operations ! N*M*(M-1) dot operations
DO i=n,1,-1 DO i=n,1,-1
rxw=w(1,i) rxw=w(1,i)
...@@ -483,8 +564,28 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M ...@@ -483,8 +564,28 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M
END DO END DO
END DO END DO
RETURN RETURN
END SUBROUTINE dbciel
!> VS = band part of inverse symmetric matrix
!!
!! \param [in,out] W (decomposition of) )symmetric band matrix
!! \param [in] MP1 band width (M) + 1
!! \param [in] N size
!! \param [out] VS band part of inverse symmetric matrix
ENTRY dbcinb(w,mp1,n, vs) ! VS = band part of inverse symmetric matrix SUBROUTINE dbcinb(w,mp1,n, vs)
USE mpdef
implicit none
INTEGER(mpi) :: i
INTEGER(mpi) :: j
INTEGER(mpi) :: k
REAL(mpd) :: rxw
REAL(mpd), INTENT(IN) :: w(mp1,n)
INTEGER(mpi), INTENT(IN) :: mp1
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(OUT) :: vs(*)
! N*M*(M-1) dot operations ! N*M*(M-1) dot operations
DO i=n,1,-1 DO i=n,1,-1
rxw=w(1,i) rxw=w(1,i)
...@@ -500,8 +601,28 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M ...@@ -500,8 +601,28 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M
! END DO ! END DO
END DO END DO
RETURN RETURN
END SUBROUTINE dbcinb
!> VS = inverse symmetric matrix
!!
!! \param [in,out] W (decomposition of) )symmetric band matrix
!! \param [in] MP1 band width (M) + 1
!! \param [in] N size
!! \param [out] VS inverse symmetric matrix
ENTRY dbcinv(w,mp1,n, vs) ! V = inverse symmetric matrix SUBROUTINE dbcinv(w,mp1,n, vs)
USE mpdef
implicit none
INTEGER(mpi) :: i
INTEGER(mpi) :: j
INTEGER(mpi) :: k
REAL(mpd) :: rxw
REAL(mpd), INTENT(IN) :: w(mp1,n)
INTEGER(mpi), INTENT(IN) :: mp1
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(OUT) :: vs(*)
! N*N/2*(M-1) dot operations ! N*N/2*(M-1) dot operations
DO i=n,1,-1 DO i=n,1,-1
rxw=w(1,i) rxw=w(1,i)
...@@ -514,7 +635,7 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M ...@@ -514,7 +635,7 @@ SUBROUTINE dbcdec(w,mp1,n, aux) ! decomposition, bandwidth M
END DO END DO
END DO END DO
RETURN RETURN
END SUBROUTINE dbcdec END SUBROUTINE dbcinv
!> Print corr band matrix and pars. !> Print corr band matrix and pars.
!! !!
...@@ -612,26 +733,24 @@ END SUBROUTINE dbcprb ...@@ -612,26 +733,24 @@ END SUBROUTINE dbcprb
!! all 1 and not stored. The other elements of L are stored in the !! all 1 and not stored. The other elements of L are stored in the
!! corresponding elements of W. !! corresponding elements of W.
!! !!
!! ENTRY DB2SLV(W,N,B, X), solution B -> X \n !! SUBROUTINE DB2SLV(W,N,B, X), solution B -> X \n
!! ENTRY DB2IEL(W,N, V), V = inverse band matrix elements !! SUBROUTINE DB2IEL(W,N, V), V = inverse band matrix elements
!! !!
!! \param [in,out] W symmetric band matrix !! \param [in,out] W symmetric band matrix
!! \param [in] N size !! \param [in] N size
!! \param [in] AUX scratch array !! \param [in] AUX scratch array
!!
SUBROUTINE db2dec(w,n, aux) SUBROUTINE db2dec(w,n, aux)
USE mpdef USE mpdef
implicit none implicit none
INTEGER(mpi) :: i INTEGER(mpi) :: i
REAL(mpd) :: rxw
REAL(mpd), INTENT(IN OUT) :: w(2,n) REAL(mpd), INTENT(IN OUT) :: w(2,n)
INTEGER(mpi), INTENT(IN OUT) :: n INTEGER(mpi), INTENT(IN OUT) :: n
REAL(mpd), INTENT(OUT) :: aux(n) REAL(mpd), INTENT(OUT) :: aux(n)
REAL(mpd) :: v(2,n),b(n),x(n), rxw
DO i=1,n DO i=1,n
aux(i)=16.0_mpd*w(1,i) ! save diagonal elements aux(i)=16.0_mpd*w(1,i) ! save diagonal elements
END DO END DO
...@@ -652,8 +771,27 @@ SUBROUTINE db2dec(w,n, aux) ...@@ -652,8 +771,27 @@ SUBROUTINE db2dec(w,n, aux)
w(1,n)=0.0_mpd w(1,n)=0.0_mpd
END IF END IF
RETURN RETURN
END SUBROUTINE db2dec
!> solution B -> X
!!
!! \param [in,out] W (decomposition of) )symmetric band matrix
!! \param [in] N size
!! \param [in] B r.h.s.
!! \param [out] X solution
SUBROUTINE db2slv(w,n,b, x)
USE mpdef
implicit none
INTEGER(mpi) :: i
REAL(mpd), INTENT(IN) :: w(2,n)
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(IN) :: b(n)
REAL(mpd), INTENT(OUT) :: x(n)
ENTRY db2slv(w,n,b, x) ! solution B -> X
! The equation W(original)*X=B is solved for X; input is B in vector X. ! The equation W(original)*X=B is solved for X; input is B in vector X.
DO i=1,n DO i=1,n
x(i)=b(i) x(i)=b(i)
...@@ -666,8 +804,24 @@ SUBROUTINE db2dec(w,n, aux) ...@@ -666,8 +804,24 @@ SUBROUTINE db2dec(w,n, aux)
x(i)=x(i)*w(1,i)-w(2,i)*x(i+1) x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)
END DO END DO
RETURN RETURN
END SUBROUTINE db2slv
!> V = inverse band matrix elements
!!
!! \param [in,out] W (decomposition of) )symmetric band matrix
!! \param [in] N size
!! \param [out] V inverse band matrix elements
SUBROUTINE db2iel(w,n, v)
USE mpdef
implicit none
INTEGER(mpi) :: i
REAL(mpd), INTENT(IN) :: w(2,n)
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(OUT) :: v(2,n)
ENTRY db2iel(w,n, v) ! V = inverse band matrix elements
! The band elements of the inverse of W(original) are calculated, ! The band elements of the inverse of W(original) are calculated,
! and stored in V in the same order as in W. ! and stored in V in the same order as in W.
! Remaining elements of the inverse are not calculated. ! Remaining elements of the inverse are not calculated.
...@@ -680,7 +834,8 @@ SUBROUTINE db2dec(w,n, aux) ...@@ -680,7 +834,8 @@ SUBROUTINE db2dec(w,n, aux)
v(1,2)= w(1,2)-v(2,2)*w(2,2) v(1,2)= w(1,2)-v(2,2)*w(2,2)
v(2,1)=-v(1,2)*w(2,1) v(2,1)=-v(1,2)*w(2,1)
v(1,1)= w(1,1)-v(2,1)*w(2,1) v(1,1)= w(1,1)-v(2,1)*w(2,1)
END SUBROUTINE db2dec RETURN
END SUBROUTINE db2iel
! (4) Symmetric band matrix of band width m=2: decomposition, ! (4) Symmetric band matrix of band width m=2: decomposition,
...@@ -700,8 +855,8 @@ END SUBROUTINE db2dec ...@@ -700,8 +855,8 @@ END SUBROUTINE db2dec
!! all 1 and not stored. The other elements of L are stored in the !! all 1 and not stored. The other elements of L are stored in the
!! corresponding elements of W. !! corresponding elements of W.
!! !!
!! ENTRY DB3SLV(W,N,B, X), solution B -> X \n !! SUBROUTINE DB3SLV(W,N,B, X), solution B -> X \n
!! ENTRY DB3IEL(W,N, V), V = inverse band matrix elements !! SUBROUTINE DB3IEL(W,N, V), V = inverse band matrix elements
!! !!
!! \param [in,out] W symmetric band matrix !! \param [in,out] W symmetric band matrix
!! \param [in] N size !! \param [in] N size
...@@ -712,15 +867,13 @@ SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2) ...@@ -712,15 +867,13 @@ SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2)
implicit none implicit none
INTEGER(mpi) :: i INTEGER(mpi) :: i
REAL(mpd) :: rxw
REAL(mpd), INTENT(IN OUT) :: w(3,n) REAL(mpd), INTENT(IN OUT) :: w(3,n)
INTEGER(mpi), INTENT(IN OUT) :: n INTEGER(mpi), INTENT(IN OUT) :: n
REAL(mpd), INTENT(OUT) :: aux(n) REAL(mpd), INTENT(OUT) :: aux(n)
! decompos ! decompos
REAL(mpd) :: v(3,n),b(n),x(n), rxw
DO i=1,n DO i=1,n
aux(i)=16.0_mpd*w(1,i) ! save diagonal elements aux(i)=16.0_mpd*w(1,i) ! save diagonal elements
END DO END DO
...@@ -755,8 +908,26 @@ SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2) ...@@ -755,8 +908,26 @@ SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2)
w(1,n)=0.0_mpd w(1,n)=0.0_mpd
END IF END IF
RETURN RETURN
END SUBROUTINE db3dec
!> solution B -> X
!!
!! \param [in,out] W (decomposition of) )symmetric band matrix
!! \param [in] N size
!! \param [in] B r.h.s.
!! \param [out] X solution
SUBROUTINE db3slv(w,n,b, x)
USE mpdef
implicit none
INTEGER(mpi) :: i
REAL(mpd), INTENT(IN) :: w(3,n)
INTEGER(mpi), INTENT(IN) :: n
REAL(mpd), INTENT(IN) :: b(n)
REAL(mpd), INTENT(OUT) :: x(n)
ENTRY db3slv(w,n,b, x) ! solution B -> X
DO i=1,n DO i=1,n
x(i)=b(i) x(i)=b(i)
END DO END DO
...@@ -771,8 +942,24 @@ SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2) ...@@ -771,8 +942,24 @@ SUBROUTINE db3dec(w,n, aux) ! decomposition (M=2)
x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)-w(3,i)*x(i+2) x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)-w(3,i)*x(i+2)
END DO END DO