Dbandmatrix.f90

Go to the documentation of this file.

! Code converted using TO_F90 by Alan Miller
! Date: 2012-04-18  Time: 19:56:08



!     (1) Symmetric matrix W: decomposition, solution, inverse


SUBROUTINE dchdec(w,n, aux)
    implicit none
    integer :: i
    integer :: ii
    integer :: j
    integer :: jj
    integer :: k
    integer :: kk
    integer :: l
    integer :: m


    DOUBLE PRECISION, INTENT(IN OUT)         :: w(*)
    INTEGER, INTENT(IN)                      :: n
    DOUBLE PRECISION, INTENT(OUT)            :: aux(n)

    DOUBLE PRECISION :: b(*),x(*),v(*),suma,ratio
    !     ...
    DO i=1,n
        aux(i)=16.0d0*w((i*i+i)/2) ! save diagonal elements
    END DO
    ii=0
    DO i=1,n
        ii=ii+i
        IF(w(ii)+aux(i) /= aux(i)) THEN     ! GT
            w(ii)=1.0d0/w(ii)                ! (I,I) div !
        ELSE
            w(ii)=0.0d0
        END IF
        jj=ii
        DO j=i+1,n
            ratio=w(i+jj)*w(ii)              ! (I,J) (I,I)
            kk=jj
            DO k=j,n
                w(kk+j)=w(kk+j)-w(kk+i)*ratio   ! (K,J) (K,I)
                kk=kk+k
            END DO ! K
            w(i+jj)=ratio                    ! (I,J)
            jj=jj+j
        END DO ! J
    END DO ! I


    return

    entry dchslv(w,n,b, x)       ! solution B -> X
    WRITE(*,*) 'before copy'
    DO i=1,n
        x(i)=b(i)
    END DO
    WRITE(*,*) 'after copy'
    ii=0
    DO i=1,n
        suma=x(i)
        DO k=1,i-1
            suma=suma-w(k+ii)*x(k)             ! (K,I)
        END DO
        x(i)=suma
        ii=ii+i
    END DO
    WRITE(*,*) 'after forward'
    DO i=n,1,-1
        suma=x(i)*w(ii)                    ! (I,I)
        kk=ii
        DO k=i+1,n
            suma=suma-w(kk+i)*x(k)             ! (K,I)
            kk=kk+k
        END DO
        x(i)=suma
        ii=ii-i
    END DO
    WRITE(*,*) 'after backward'
    return

    entry dchinv(w,n,v)         ! inversion
    ii=(n*n-n)/2
    DO i=n,1,-1
        suma=w(ii+i)                       ! (I,I)
        DO j=i,1,-1
            DO k=j+1,n
                l=min(i,k)
                m=max(i,k)
                suma=suma-w(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K)
            END DO
            v(ii+j)=suma                      ! (I,J)
            suma=0.0d0
        END DO
        ii=ii-i+1
    END DO
END SUBROUTINE dchdec


REAL FUNCTION condes(w,n,aux)
    implicit none
    real :: cond
    integer :: i
    integer :: ir
    integer :: is
    integer :: k
    integer :: kk
    real :: xla1
    real :: xlan


    DOUBLE PRECISION, INTENT(IN)             :: w(*)
    INTEGER, INTENT(IN)                      :: n
    DOUBLE PRECISION, INTENT(IN OUT)         :: aux(n)

    DOUBLE PRECISION :: suma,sumu,sums

    ir=1
    is=1
    DO i=1,n
        IF(w((i*i+i)/2) < w((is*is+is)/2)) is=i ! largest Dii
        IF(w((i*i+i)/2) > w((ir*ir+ir)/2)) ir=i ! smallest Dii
    END DO

    sumu=0.0     ! find smallest eigenvalue
    sums=0.0
    DO i=n,1,-1  ! backward
        suma=0.0
        IF(i == ir) suma=1.0d0
        kk=(i*i+i)/2
        DO k=i+1,n
            suma=suma-w(kk+i)*aux(k)             ! (K,I)
            kk=kk+k
        END DO
        aux(i)=suma
        sumu=sumu+aux(i)*aux(i)
    END DO
    xlan=sngl(w((ir*ir+ir)/2)*dsqrt(sumu))
    IF(xlan /= 0.0) xlan=1.0/xlan

    DO i=1,n
        IF(i == is) THEN
            sums=1.0d0
        ELSE IF(i > is) THEN
            sums=sums+w(is+(i*i-i)/2)**2
        END IF
    END DO       ! is Ws
    xla1=0.0
    IF(w((is*is+is)/2) /= 0.0) xla1=sngl(dsqrt(sums)/w((is*is+is)/2))

    cond=0.0
    IF(xla1 > 0.0.AND.xlan > 0.0) cond=xla1/xlan
    !     estimated condition
    condes=cond
END FUNCTION condes


!     (2) Symmetric band matrix: decomposition, solution, complete
!                                inverse and band part of the inverse

SUBROUTINE dbcdec(w,mp1,n, aux)  ! decomposition, bandwidth M
    implicit none
    integer :: i
    integer :: j
    integer :: k
    !     M=MP1-1                                    N*M(M-1) dot operations

    DOUBLE PRECISION, INTENT(IN OUT)         :: w(mp1,n)
    INTEGER, INTENT(IN)                      :: mp1
    INTEGER, INTENT(IN)                      :: n
    DOUBLE PRECISION, INTENT(OUT)            :: aux(n)
    ! decompos
    DOUBLE PRECISION :: v(mp1,n),b(n),x(n), vs(*),rxw
    !     ...
    DO i=1,n
        aux(i)=16.0d0*w(1,i) ! save diagonal elements
    END DO
    DO i=1,n
        IF(w(1,i)+aux(i) /= aux(i)) THEN
            w(1,i)=1.0/w(1,i)
        ELSE
            w(1,i)=0.0d0  ! singular
        END IF
        DO j=1,min(mp1-1,n-i)
            rxw=w(j+1,i)*w(1,i)
            DO k=1,min(mp1-1,n-i)+1-j
                w(k,i+j)=w(k,i+j)-w(k+j,i)*rxw
            END DO
            w(j+1,i)=rxw
        END DO
    END DO
    return

    entry dbcslv(w,mp1,n,b, x)  ! solution B -> X
    !                                                N*(2M-1) dot operations
    DO i=1,n
        x(i)=b(i)
    END DO
    DO i=1,n ! forward substitution
        DO j=1,min(mp1-1,n-i)
            x(j+i)=x(j+i)-w(j+1,i)*x(i)
        END DO
    END DO
    DO i=n,1,-1 ! backward substitution
        rxw=x(i)*w(1,i)
        DO j=1,min(mp1-1,n-i)
            rxw=rxw-w(j+1,i)*x(j+i)
        END DO
        x(i)=rxw
    END DO
    return

    entry dbciel(w,mp1,n, v)    ! V = inverse band matrix elements
    !                                               N*M*(M-1) dot operations
    DO i=n,1,-1
        rxw=w(1,i)
        DO j=i,max(1,i-mp1+1),-1
            DO k=j+1,min(n,j+mp1-1)
                rxw=rxw-v(1+abs(i-k),min(i,k))*w(1+k-j,j)
            END DO
            v(1+i-j,j)=rxw
            rxw=0.0
        END DO
    END DO
    return

    entry dbcinb(w,mp1,n, vs)   ! VS = band part of inverse symmetric matrix
    !                                             N*M*(M-1) dot operations
    DO i=n,1,-1
        rxw=w(1,i)
        DO j=i,max(1,i-mp1+1),-1
            DO k=j+1,min(n,j+mp1-1)
                rxw=rxw-vs((max(i,k)*(max(i,k)-1))/2+min(i,k))*w(1+k-j,j)
            END DO
            vs((i*i-i)/2+j)=rxw
            rxw=0.0
        END DO
    !       DO J=MAX(1,I-MP1+1)-1,1,-1
    !        VS((I*I-I)/2+J)=0.0
    !       END DO
    END DO
    return

    entry dbcinv(w,mp1,n, vs)   ! V = inverse symmetric matrix
    !                                             N*N/2*(M-1) dot operations
    DO i=n,1,-1
        rxw=w(1,i)
        DO j=i,1,-1
            DO k=j+1,min(n,j+mp1-1)
                rxw=rxw-vs((max(i,k)*(max(i,k)-1))/2+min(i,k))*w(1+k-j,j)
            END DO
            vs((i*i-i)/2+j)=rxw
            rxw=0.0
        END DO
    END DO
    return
END SUBROUTINE dbcdec


SUBROUTINE dbcprv(w,mp1,n,b)
    implicit none
    integer :: i
    integer :: j
    integer :: nj
    real :: rho


    DOUBLE PRECISION, INTENT(IN OUT)         :: w(mp1,n)
    INTEGER, INTENT(IN)                      :: mp1
    INTEGER, INTENT(IN)                      :: n
    DOUBLE PRECISION, INTENT(OUT)            :: b(n)

    DOUBLE PRECISION :: err
    INTEGER :: irho(5)
    !     ...
    WRITE(*,*) ' '
    WRITE(*,101)

    DO i=1,n
        err=dsqrt(w(1,i))
        nj=0
        DO j=2,mp1
            IF(i+1-j > 0.AND.nj < 5) THEN
                nj=nj+1
                rho=sngl(w(j,i+1-j)/(err*dsqrt(w(1,i+1-j))))
                irho(nj)=ifix(100.0*abs(rho)+0.5)
                IF(rho < 0.0) irho(nj)=-irho(nj)
            END IF
        END DO
        WRITE(*,102) i,b(i),err,(irho(j),j=1,nj)
    END DO
    WRITE(*,103)
101 format(5x,'i   Param',7x,'error',7x,'  c(i,i-1) c(i,i-2)'/)
102 format(1x,i5,2g12.4,1x,5i9)
103 format(33x,'(correlation coefficients in percent)')
END SUBROUTINE dbcprv


SUBROUTINE dbcprb(w,mp1,n)
    implicit none
    integer :: i
    integer :: j


    DOUBLE PRECISION, INTENT(IN OUT)         :: w(mp1,n)
    INTEGER, INTENT(IN)                      :: mp1
    INTEGER, INTENT(IN)                      :: n


    !     ...
    IF(mp1 > 6) return
    WRITE(*,*) ' '
    WRITE(*,101)
    DO i=1,n
        WRITE(*,102) i,(w(j,i),j=1,mp1)
    END DO
    WRITE(*,*) ' '
101 format(5x,'i   Diag  ')
102 format(1x,i5,6g12.4)
END SUBROUTINE dbcprb


!     (3) Symmetric band matrix of band width m=1: decomposition,
!            solution, complete, inverse and band part of the inverse


SUBROUTINE db2dec(w,n, aux)
    implicit none
    integer :: i


    DOUBLE PRECISION, INTENT(IN OUT)         :: w(2,n)
    INTEGER, INTENT(IN OUT)                  :: n
    DOUBLE PRECISION, INTENT(OUT)            :: aux(n)

    DOUBLE PRECISION :: v(2,n),b(n),x(n), rxw

    DO i=1,n
        aux(i)=16.0d0*w(1,i) ! save diagonal elements
    END DO
    DO i=1,n-1
        IF(w(1,i)+aux(i) /= aux(i)) THEN
            w(1,i)=1.0d0/w(1,i)
            rxw=w(2,i)*w(1,i)
            w(1,i+1)=w(1,i+1)-w(2,i)*rxw
            w(2,i)=rxw
        ELSE ! singular
            w(1,i)=0.0d0
            w(2,i)=0.0d0
        END IF
    END DO
    IF(w(1,n)+aux(n) > aux(n)) THEN           ! N
        w(1,n)=1.0d0/w(1,n)
    ELSE ! singular
        w(1,n)=0.0d0
    END IF
    return

    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.
    DO i=1,n
        x(i)=b(i)
    END DO
    DO i=1,n-1           ! by forward substitution
        x(i+1)=x(i+1)-w(2,i)*x(i)
    END DO
    x(n)=x(n)*w(1,n)     ! and backward substitution
    DO i=n-1,1,-1
        x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)
    END DO
    return

    entry db2iel(w,n, v)        ! V = inverse band matrix elements
    !     The band elements of the inverse of W(original) are calculated,
    !     and stored in V in the same order as in W.
    !     Remaining elements of the inverse are not calculated.
    v(1,n  )= w(1,n)
    v(2,n-1)=-v(1,n  )*w(2,n-1)
    DO i=n-1,3,-1
        v(1,i  )= w(1,i  )-v(2,i  )*w(2,i  )
        v(2,i-1)=-v(1,i  )*w(2,i-1)
    END DO
    v(1,2)= w(1,2)-v(2,2)*w(2,2)
    v(2,1)=-v(1,2)*w(2,1)
    v(1,1)= w(1,1)-v(2,1)*w(2,1)
END SUBROUTINE db2dec


!     (4) Symmetric band matrix of band width m=2: decomposition,
!            solution, complete, inverse and band part of the inverse


SUBROUTINE db3dec(w,n, aux)      ! decomposition (M=2)
    implicit none
    integer :: i


    DOUBLE PRECISION, INTENT(IN OUT)         :: w(3,n)
    INTEGER, INTENT(IN OUT)                  :: n
    DOUBLE PRECISION, INTENT(OUT)            :: aux(n)
    ! decompos

    DOUBLE PRECISION :: v(3,n),b(n),x(n), rxw

    DO i=1,n
        aux(i)=16.0d0*w(1,i) ! save diagonal elements
    END DO
    DO i=1,n-2
        IF(w(1,i)+aux(i) /= aux(i)) THEN
            w(1,i)=1.0d0/w(1,i)
            rxw=w(2,i)*w(1,i)
            w(1,i+1)=w(1,i+1)-w(2,i)*rxw
            w(2,i+1)=w(2,i+1)-w(3,i)*rxw
            w(2,i)=rxw
            rxw=w(3,i)*w(1,i)
            w(1,i+2)=w(1,i+2)-w(3,i)*rxw
            w(3,i)=rxw
        ELSE ! singular
            w(1,i)=0.0d0
            w(2,i)=0.0d0
            w(3,i)=0.0d0
        END IF
    END DO
    IF(w(1,n-1)+aux(n-1) > aux(n-1)) THEN
        w(1,n-1)=1.0d0/w(1,n-1)
        rxw=w(2,n-1)*w(1,n-1)
        w(1,n)=w(1,n)-w(2,n-1)*rxw
        w(2,n-1)=rxw
    ELSE ! singular
        w(1,n-1)=0.0d0
        w(2,n-1)=0.0d0
    END IF
    IF(w(1,n)+aux(n) > aux(n)) THEN
        w(1,n)=1.0d0/w(1,n)
    ELSE ! singular
        w(1,n)=0.0d0
    END IF
    return

    entry db3slv(w,n,b, x)      ! solution B -> X
    DO i=1,n
        x(i)=b(i)
    END DO
    DO i=1,n-2           ! by forward substitution
        x(i+1)=x(i+1)-w(2,i)*x(i)
        x(i+2)=x(i+2)-w(3,i)*x(i)
    END DO
    x(n)=x(n)-w(2,n-1)*x(n-1)
    x(n)=x(n)*w(1,n)     ! and backward substitution
    x(n-1)=x(n-1)*w(1,n-1)-w(2,n-1)*x(n)
    DO i=n-2,1,-1
        x(i)=x(i)*w(1,i)-w(2,i)*x(i+1)-w(3,i)*x(i+2)
    END DO
    return

    entry db3iel(w,n, v)        ! V = inverse band matrix elements
    !     The band elements of the inverse of W(original) are calculated,
    !     and stored in V in the same order as in W.
    !     Remaining elements of the inverse are not calculated.
    v(1,n  )= w(1,n)
    v(2,n-1)=-v(1,n  )*w(2,n-1)
    v(3,n-2)=-v(2,n-1)*w(2,n-2)-v(1,n  )*w(3,n-2)
    v(1,n-1)= w(1,n-1) -v(2,n-1)*w(2,n-1)
    v(2,n-2)=-v(1,n-1)*w(2,n-2)-v(2,n-1)*w(3,n-2)
    v(3,n-3)=-v(2,n-2)*w(2,n-3)-v(1,n-1)*w(3,n-3)
    DO i=n-2,3,-1
        v(1,i  )= w(1,i  ) -v(2,i  )*w(2,i  )-v(3,i)*w(3,i  )
        v(2,i-1)=-v(1,i  )*w(2,i-1)-v(2,i)*w(3,i-1)
        v(3,i-2)=-v(2,i-1)*w(2,i-2)-v(1,i)*w(3,i-2)
    END DO
    v(1,2)= w(1,2) -v(2,2)*w(2,2)-v(3,2)*w(3,2)
    v(2,1)=-v(1,2)*w(2,1)-v(2,2)*w(3,1)
    v(1,1)= w(1,1) -v(2,1)*w(2,1)-v(3,1)*w(3,1)
END SUBROUTINE db3dec


!     (5) Symmetric matrix with band structure, bordered by full row/col
!          - is not yet included -

!      SUBROUTINE BSOLV1(N,CU,RU,CK,RK,CH,     BK,UH,   AU) ! 1
!     Input:  CU = 3*N array         replaced by decomposition
!             RU   N array rhs
!             CK   diagonal element
!             RK   rhs
!             CH   N-vector

!     Aux:    AU   N-vector auxliliary array

!     Result: FK   curvature
!             BK   variance
!             UH   smoothed data points


!      DOUBLE PRECISION CU(3,N),CI(3,N),CK,BK,AU(N),UH(N)
!     ...
!      CALL BDADEC(CU,3,N, AU)    ! decomposition
!      CALL DBASLV(CU,3,N, RU,UH)  ! solve for zero curvature
!      CALL DBASLV(CU,3,N, CH,AU)  ! solve for aux. vector
!      CTZ=0.0D0
!      ZRU=0.0D0
!      DO I=1,N
!       CTZ=CTZ+CH(I)*AU(I)        ! cT z
!       ZRU=ZRU+RY(I)*AU(I)        ! zT ru
!      END DO
!      BK=1.0D0/(CK-CTZ)           ! variance of curvature
!      FK=BK   *(RK-ZRU)           ! curvature
!      DO I=1,N
!       UH(I)=UH(I)-FK*AU(I)       ! smoothed data points
!      END DO
!      RETURN

!      ENTRY BINV1(N,CU,CI, FK,AU)
!      DOUBLE PRECISION CI(3,N)
!     ...
!      CALL DBAIBM(CU,3,N, CI)           ! block part of inverse
!      DO I=1,N
!       CI(1,I)=CI(1,I)+FK*AU(I)*AU(I)   ! diagonal elements
!       IF(I.LT.N) CI(2,I)=CI(2,I)+FK*AU(I)*AU(I+1) ! next diagonal
!       IF(I.LT.N-1) CI(3,I)=CI(3,I)+FK*AU(I)*AU(I+2) ! next diagonal
!      END DO

!      END


SUBROUTINE dcfdec(w,n)

    IMPLICIT NONE
    DOUBLE PRECISION, INTENT(OUT)            :: w(*)
    INTEGER, INTENT(IN)                      :: n
    INTEGER :: i,j,k
    DOUBLE PRECISION :: epsm,gamm,xchi,beta,delta,theta

    epsm=2.2d-16 ! machine precision
    gamm=0.0d0   ! max diagonal element
    xchi=0.0d0   ! max off-diagonal element
    DO k=1,n
        gamm=max(gamm,abs(w((k*k+k)/2)))
        DO j=k+1,n
            xchi=max(xchi,abs(w((j*j-j)/2+k)))
        END DO
    END DO
    beta=sqrt(max(gamm,xchi/max(1.0d0,sqrt(dfloat(n*n-1))),epsm))
    delta=epsm*max(1.0d0,gamm+xchi)

    DO k=1,n
        DO i=1,k-1
            w((k*k-k)/2+i)=w((k*k-k)/2+i)*w((i*i+i)/2)
        END DO
        DO j=k+1,n
            DO i=1,k-1
                w((j*j-j)/2+k)=w((j*j-j)/2+k)-w((k*k-k)/2+i)*w((j*j-j)/2+i)
            END DO
        END DO
        theta=0.0d0
        DO j=k+1,n
            theta=max(theta,abs(w((j*j-j)/2+k)))
        END DO
        w((k*k+k)/2)=1.0d0/max(abs(w((k*k+k)/2)),(theta/beta)**2,delta)
        DO j=k+1,n
            w((j*j+j)/2)=w((j*j+j)/2)-w((j*j-j)/2+k)**2*w((k*k+k)/2)
        END DO
    END DO ! K

END SUBROUTINE dcfdec


SUBROUTINE dbfdec(w,mp1,n)

    IMPLICIT NONE
    DOUBLE PRECISION, INTENT(OUT)            :: w(mp1,n)
    INTEGER, INTENT(IN OUT)                  :: mp1
    INTEGER, INTENT(IN)                      :: n
    INTEGER :: i,j,k
    DOUBLE PRECISION :: epsm,gamm,xchi,beta,delta,theta

    epsm=2.2d-16 ! machine precision
    gamm=0.0d0   ! max diagonal element
    xchi=0.0d0   ! max off-diagonal element
    DO k=1,n
        gamm=max(gamm,abs(w(1,k)))
        DO j=2,min(mp1,n-k+1)
            xchi=max(xchi,abs(w(j,k)))
        END DO
    END DO
    beta=sqrt(max(gamm,xchi/max(1.0d0,sqrt(dfloat(n*n-1))),epsm))
    delta=epsm*max(1.0d0,gamm+xchi)

    DO k=1,n
        DO i=2,min(mp1,k)
            w(i,k-i+1)=w(i,k-i+1)*w(1,k-i+1)
        END DO
        DO j=2,min(mp1,n-k+1)
            DO i=max(2,j+k+1-mp1),k
                w(j,k)=w(j,k)-w(k-i+2,i-1)*w(j-i+k+1,i-1)
            END DO
        END DO
        theta=0.0d0
        DO j=2,min(mp1,n-k+1)
            theta=max(theta,abs(w(j,k)))
        END DO
        w(1,k)=1.0d0/max(abs(w(1,k)),(theta/beta)**2,delta)
        DO j=2,min(mp1,n-k+1)
            w(1,k+j-1)=w(1,k+j-1)-w(1,k)*w(j,k)**2
        END DO
    END DO ! K

END SUBROUTINE dbfdec