Dbandmatrix.f

Go to the documentation of this file.

 

*     (1) symmetric matrix w: decomposition, solution, inverse

*
      SUBROUTINE dchdec(W,N, AUX)
*
      INTEGER n
      DOUBLE PRECISION w(*),aux(n),b(*),x(*),v(*),sum,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).NE.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
       sum=x(i)
       DO k=1,i-1
        sum=sum-w(k+ii)*x(k)             ! (K,I)
       END DO
       x(i)=sum
       ii=ii+i
      END DO
      WRITE(*,*) 'after forward'
      DO i=n,1,-1
       sum=x(i)*w(ii)                    ! (I,I)
       kk=ii
       DO k=i+1,n
        sum=sum-w(kk+i)*x(k)             ! (K,I)
        kk=kk+k
       END DO
       x(i)=sum
       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
       sum=w(ii+i)                       ! (I,I)
       DO j=i,1,-1
        DO k=j+1,n
         l=min(i,k)
         m=max(i,k)
         sum=sum-w(j+(k*k-k)/2)*v(l+(m*m-m)/2) ! (J,K) (I,K)
        END DO
        v(ii+j)=sum                      ! (I,J)
        sum=0.0d0
       END DO
       ii=ii-i+1
      END DO
      END

*
      FUNCTION condes(W,N,AUX)
*
      INTEGER n
      DOUBLE PRECISION w(*),aux(n),sum,sumu,sums
*
      ir=1
      is=1
      DO i=1,n
       IF(w((i*i+i)/2).LT.w((is*is+is)/2)) is=i ! largest Dii
       IF(w((i*i+i)/2).GT.w((ir*ir+ir)/2)) ir=i ! smallest Dii
      END DO

      sumu=0.0     ! find smallest eigenvalue        
      DO i=n,1,-1  ! backward 
       sum=0.0
       IF(i.EQ.ir) sum=1.0d0
       kk=(i*i+i)/2
       DO k=i+1,n
        sum=sum-w(kk+i)*aux(k)             ! (K,I)
        kk=kk+k
       END DO
       aux(i)=sum
       sumu=sumu+aux(i)*aux(i)
      END DO
      xlan=sngl(w((ir*ir+ir)/2)*dsqrt(sumu))
      IF(xlan.NE.0.0) xlan=1.0/xlan
      
      DO i=1,n
       IF(i.EQ.is) THEN
          sums=1.0d0
       ELSE IF(i.GT.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).NE.0.0) xla1=sngl(dsqrt(sums)/w((is*is+is)/2))
      
      cond=0.0
      IF(xla1.GT.0.0.AND.xlan.GT.0.0) cond=xla1/xlan 
*     estimated condition
      condes=cond
      END 
  

*     (2) symmetric band matrix: decomposition, solution, complete
*                                inverse and band part of the inverse
*
      SUBROUTINE dbcdec(W,MP1,N, AUX)  ! decomposition, bandwidth M
*     m=mp1-1                                    n*m(m-1) dot operations
      INTEGER mp1,n
      DOUBLE PRECISION w(mp1,n),v(mp1,n),b(n),x(n),aux(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).NE.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
c       DO j=max(1,i-mp1+1)-1,1,-1
c        vs((i*i-i)/2+j)=0.0
c       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 dbcprv(W,MP1,N,B)
*
      INTEGER mp1,n
      DOUBLE PRECISION w(mp1,n),b(n),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.GT.0.AND.nj.LT.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.LT.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 dbcprb(W,MP1,N)
*
      INTEGER mp1,n
      DOUBLE PRECISION w(mp1,n)
*     ...
      IF(mp1.GT.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

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

*
      SUBROUTINE db2dec(W,N, AUX)
*
      INTEGER n
      DOUBLE PRECISION w(2,n),v(2,n),b(n),x(n),aux(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).NE.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).GT.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

   
*     (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) 
*
      INTEGER n
      DOUBLE PRECISION w(3,n),v(3,n),b(n),x(n),aux(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).NE.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).GT.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).GT.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


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

c      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
*             
         
c      DOUBLE PRECISION cu(3,n),ci(3,n),ck,bk,au(n),uh(n)
*     ...
c      CALL bdadec(cu,3,n, au)    ! decomposition
c      CALL dbaslv(cu,3,n, ru,uh)  ! solve for zero curvature
c      CALL dbaslv(cu,3,n, ch,au)  ! solve for aux. vector
c      ctz=0.0d0
c      zru=0.0d0
c      DO i=1,n
c       ctz=ctz+ch(i)*au(i)        ! cT z
c       zru=zru+ry(i)*au(i)        ! zT ru       
c      END DO 
c      bk=1.0d0/(ck-ctz)           ! variance of curvature
c      fk=bk   *(rk-zru)           ! curvature
c      DO i=1,n
c       uh(i)=uh(i)-fk*au(i)       ! smoothed data points
c      END DO
c      return

c      entry binv1(n,cu,ci, fk,au)
c      DOUBLE PRECISION ci(3,n)
*     ...
c      CALL dbaibm(cu,3,n, ci)           ! block part of inverse
c      DO i=1,n
c       ci(1,i)=ci(1,i)+fk*au(i)*au(i)   ! diagonal elements
c       IF(i.LT.n) ci(2,i)=ci(2,i)+fk*au(i)*au(i+1) ! next diagonal
c       IF(i.LT.n-1) ci(3,i)=ci(3,i)+fk*au(i)*au(i+2) ! next diagonal 
c      END DO

c      END 

*
      SUBROUTINE dcfdec(W,N)

      IMPLICIT NONE
      INTEGER n,i,j,k
      DOUBLE PRECISION w(*),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 dbfdec(W,MP1,N)
*
      IMPLICIT NONE
      INTEGER mp1,n,i,j,k  
      DOUBLE PRECISION w(mp1,n),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