Dbandmatrix.f90 File Reference

Symmetric (band) matrix routines. More...

Go to the source code of this file.

Functions/Subroutines
subroutine	dchdec (w, n, aux)
	Decomposition of symmetric matrix.
real(mps) function	condes (w, n, aux)
	Etimate condition.
subroutine	dbcdec (w, mp1, n, aux)
	Decomposition of symmetric band matrix.
subroutine	dbcprv (w, mp1, n, b)
	Print corr band matrix and pars.
subroutine	dbcprb (w, mp1, n)
	Print band matrix.
subroutine	db2dec (w, n, aux)
	Decomposition (M=1).
subroutine	db3dec (w, n, aux)
	Decomposition (M=2).
subroutine	dcfdec (w, n)
	Decomposition of symmetric matrix.
subroutine	dbfdec (w, mp1, n)
	Decomposition of symmetric band matrix.

Detailed Description

Symmetric (band) matrix routines.

Author

Volker Blobel, University Hamburg, 2005-2009 (initial Fortran77 version)

Claus Kleinwort, DESY (maintenance and developement)

Copyright

Copyright (c) 2009 - 2015 Deutsches Elektronen-Synchroton, Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY

This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details.

You should have received a copy of the GNU Library General Public License along with this program (see the file COPYING.LIB for more details); if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

For the original broken lines implementation by V. Blobel (University Hamburg).

    *************************************************************
    *                                                           *
    *   Subroutines for symmetric and symmetric band matrices,  *
    *   based on the (square root free) Cholesky decomposition. *
    *                                                           *
    *************************************************************

    All floating point arguments are in DOUBLE PRECISION (and all
    entry names start with a D).

    The Cholesky decomposition transforms a symmetric matrix W
    e.g. the matrix from normal equation of least squares,
    according to
                          W = L D L^         (L^ means L transposed)
    where D is a diagonal matrix and L is a unit triangular matrix
    (diagonal elements all ones, all elements above diagonal zero).

    The above decomposition allows to solve a matrix equation
                        W x = b
    in two steps, using an auxiliary vector y:

                 solve  L y = b  for y, and

               solve D L^ x = y  for x.

    The inverse matrix of W can be calculated from the decomposition.

    In least-squares normal equations the inverse matrix is equal to
    the covariance matrix of the fitted parameters. All diagonal elements
    of the inverse matrix, the parameter variances, are positive, and
    the matrix is positive-definite (all eigenvalues > 0).

    The Cholesky algorithm is stable for a positive-definite matrix.
    The standard form of the Cholesky algorithm includes n square roots
    for a n-by-n matrix, and is possible only for positive-definite
    matrices. The version used here is squareroot-free; this algorithm
    does not necessarily break down in the indefinite case, although
    it is potentially unstable in this case. All decomposition routines
    include a check for singularity, and this check needs an auxiliary
    array AUX of dimension n.

    Method: The Cholesky algorithm for symmetric matrix decomposition
    makes use of the symmetry. The operation count (leading term)
    is n**3/6 (compared to n**3/3 for normal Gaussian elimination).
    The solution of the two triangular systems involves operations
    proportional to n**2.

    The real advantage of the Cholesky algorithm is for band matrices,
    where all matrix elements outside of a band with total width
    (2m+1) around the diagonal are zero. The band structure is kept
    in the decomposition, and allows a fast solution of matrix equations.
    The operation count (leading term) is proportional to m**2*n
    and thus (for fixed m) linear in n. Thus for n=100 and m=2 the
    Cholesky algorithm for the band matrix is 1000 times faster than
    the standard solution method.

    The inverse of a band matrix is a full matrix. It is not necessary
    to calculate the inverse, if only the solution for a matrix equation
    is needed. However the inverse is often needed, because the elements
    of the inverse are the variances and covariances of parameters in
    a least-squares fit. The inverse can be calculated afterwards from
    the decomposition. Since the inverse matrix is a full matrix, this
    has of course an operation count proportional to n**3.

    Usually only the elements of the inverse in and around the diagonal
    are really needed, and this subset of inverse elements, corresponding
    to the original band, can be calculated from the decomposition with
    an operation count, which is linear in n. Thus all variances (the
    diagonal elements) and covariances between neighbour parameters
    are calculated in a short time even for large matrices.

    Matrix storage: the mathematical indexing of matrix elements follows
    the scheme:

                      (  W11   W12   W13 ... W1n  )
                      (  W21   W22   W23 ... W2n  )
                  W = (  ...   ...   ...     ...  )
                      (  ...   ...   ...     ...  )
                      (  Wn1   Wn2   Wn3 ... Wnn  )

    and a storage in an array would require n**2 words, although the
    symmetric matrix has only (n**2+n)/2 different elements, and a band
    matrix has less than (m+1)*n different elements. Therefore the
    following storage schemes are used.

    Symmetric matrix: the elements are in the order
            W11   W12   W22   W13   W23   W33   W14 ... Wnn
    with total (n**2+n)/2 array elements.

    Band matrix: a band matrix of bandwidth m is stored in an array
    of dimension W(m+1,n), according to

                      W(1,.)    W(2,.)    W(3,.)
                     --------------------------------
                       W11       W12       W13
                       W22       W23       W24
                       W33       W34       W35
                       ...
                       Wnn        -         -

    The example is for a bandwidth of m=2; three elements at the end
    are unused. The diagonal elements are in the array elements W(1,.).

    This package includes subroutines for:

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

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

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

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

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

    The subroutines for a fixed band width of m=1 and of m=2 are
    faster than the general routine, because certain loops are avoided
    and replaced by the direct code.

    Historical remark: the square-root algorithm was invented by the
    french Mathematician Andre-Louis Cholesky (1875 - 1918).
    Cholesky's method of computing solutions to the normal equations was
    published 1924, after the death of Cholesky, by Benoit.
    The method received little attention after its publication in 1924.
    In 1948 the method was analysed in a paper by Fox, Huskey and
    Wilkinson, and in the same year Turing published a paper on the
    stability of the method.

    The fast method to calculate the band part of the inverse matrix
    is usually not mentioned in the literature. An exception is:
    I.S.Duff, A.M.Erisman and J.K.Reid, Direct Methods for Sparse
    Matrices, Oxford Science Publications, 1986.
    The following original work is quoted in this book:
    K.Takahashi, J.Fagan and M.Chin, Formation of a sparse bus
    impedance matrix and its application to short circuit study.
    Proceedings 8th PICA Conference, Minneapolis, Minnesota, 1973
    A.M.Erisman and W.F.Tinney, On computing certain elements of the
    inverse of a sparse matrix, CACM 18, 177-179, 1975



    symmetric          decomposit. solution    inv-element    inverse
    ----------------  |-----------|-----------|--------------|-----------|
    n x n matrix        DCHDEC      DCHSLV      -              DCHINV
    band matrix m,n     DBCDEC      DBCSLV      DBCIEL/DBCINB  DBCINV
    bandwidth m=1       DB2DEC      DB2SLV      DB2IEL         -
    bandwidth m=2       DB3DEC      DB3SLV      DB3IEL         -

    The DB2... and DB3... routines are special routines for a fixed bandwidth
    of 1 and 2, they are faster versions of the general DBG... routines.
    The complete inverse matrix can be obtained by DBGINV.
    The routine DBGPRI can be used to print all types of band matrices.

    The decomposition in routines ...DEC replaces (overwrites) the
    original matrix (the number of elements is identical). All other
    routines require W to be the already decomposed matrix.
    The matrix L is a unit lower triangular matrix, with ones on the
    diagonal, which have not be stored. Instead the inverse of the
    diagonal elements of matrix D are stored in those places.

    In the  solution routines ...SLV the array B is the right-hand matrix,
    the array is the resulting solution. The same array can be used
    for B and X.


    W(.) and V(.) are symmetric n-by-n matrices with (N*N+N)/2 elements

    SUBROUTINE DCHDEC(W,N, AUX)      ! decomposition, symmetric matrix
         ENTRY DCHSLV(W,N,B, X)      ! solution B -> X
         ENTRY DCHINV(W,N, V)        ! inversion

    SUBROUTINE DCFDEC(W,N)           ! modified composition, symmetric
                                     ! alternative to DCHDEC

    W(.) and V(.) are band matrices, n rows, band width m (i.e. the total
         width of the band is (2m+1).
         With MP1 = m +1, the array has dimension W(MP1,N).
         The symmetric matrix VS has (N*N+N)/2 elements

    SUBROUTINE DBCDEC(W,MP1,N, AUX)  ! decomposition, bandwidth M
         ENTRY DBCSLV(W,MP1,N,B, X)  ! solution B -> X
         ENTRY DBCIEL(W,MP1,N, V)    ! V = inverse band matrix elements
         ENTRY DBCINV(W,MP1,N, VS)   ! V = inverse symmetric matrix

    SUBROUTINE DBFDEC(W,MP1,N)       ! modified decomposition, bandwidth M
                                     ! alternative to DBCDEC

    SUBROUTINE DBCPRB(W,MP1,N)       ! print band matrix
    SUBROUTINE DBCPRV(W,MP1,N,B)     ! print corr band matrix and pars

    SUBROUTINE DB2DEC(W,N, AUX)      ! decomposition (M=1)
         ENTRY DB2SLV(W,N,B, X)      ! solution B -> X
         ENTRY DB2IEL(W,N, V)        ! V = inverse band matrix elements

    SUBROUTINE DB3DEC(W,N, AUX)      ! decomposition (M=2)
         ENTRY DB3SLV(W,N,B, X)      ! solution B -> X
         ENTRY DB3IEL(W,N, V)        ! V = inverse band matrix elements

Definition in file Dbandmatrix.f90.

Function/Subroutine Documentation

condes()

real(mps) function condes	(	real(mpd), dimension(*), intent(in)	w,
		integer(mpi), intent(in)	n,
		real(mpd), dimension(n), intent(inout)	aux
	)

Etimate condition.

Parameters

[in]	W	symmetric matrix
[in]	N	size
[in]	AUX	scratch array

Returns

condition

Definition at line 345 of file Dbandmatrix.f90.

    USE mpdef
 
    implicit none
    REAL(mps) :: cond
    INTEGER(mpi) :: i
    INTEGER(mpi) :: ir
    INTEGER(mpi) :: is
    INTEGER(mpi) :: k
    INTEGER(mpi) :: kk
    REAL(mps) :: xla1
    REAL(mps) :: xlan
 
 
    REAL(mpd), INTENT(IN)             :: w(*)
    INTEGER(mpi), INTENT(IN)                      :: n
    REAL(mpd), INTENT(IN OUT)         :: aux(n)
 
    REAL(mpd) :: 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.0_mpd
        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=real(w((ir*ir+ir)/2)*sqrt(sumu),mps)
    IF(xlan /= 0.0) xlan=1.0/xlan
 
    DO i=1,n
        IF(i == is) THEN
            sums=1.0_mpd
        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=real(sqrt(sums)/w((is*is+is)/2),mps)
 
    cond=0.0
    IF(xla1 > 0.0.AND.xlan > 0.0) cond=xla1/xlan
    !     estimated condition
    condes=cond

db2dec()

subroutine db2dec	(	real(mpd), dimension(2,n), intent(inout)	w,
		integer(mpi), intent(inout)	n,
		real(mpd), dimension(n), intent(out)	aux
	)

Decomposition (M=1).

W is a symmetrix positive definite band matrix of bandwidth 1(+1). W(1,.) are the diagonal elements, W(2,.) is the next diagonals; W(2,N) is never referenced. AUX is an auxiliary array of length N. W is decomposed to L D Lt, where D = diagonal and L unit triangular. A row is set to zero, if the diagonal element is reduced in previous steps by a word length (i.e. global correlation coefficient large). The resulting L and D replace W: the diagonal elements W(1,...) will contain the inverse of the D-elements; the diagonal elements of L are all 1 and not stored. The other elements of L are stored in the corresponding elements of W.

ENTRY DB2SLV(W,N,B, X), solution B -> X
ENTRY DB2IEL(W,N, V), V = inverse band matrix elements

Parameters

[in,out]	W	symmetric band matrix
[in]	N	size
[in]	AUX	scratch array

Definition at line 622 of file Dbandmatrix.f90.

    USE mpdef
 
    implicit none
    INTEGER(mpi) :: i
 
 
    REAL(mpd), INTENT(IN OUT)         :: w(2,n)
    INTEGER(mpi), INTENT(IN OUT)                  :: n
    REAL(mpd), INTENT(OUT)            :: aux(n)
 
    REAL(mpd) :: v(2,n),b(n),x(n), rxw
 
    DO i=1,n
        aux(i)=16.0_mpd*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.0_mpd/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.0_mpd
            w(2,i)=0.0_mpd
        END IF
    END DO
    IF(w(1,n)+aux(n) > aux(n)) THEN           ! N
        w(1,n)=1.0_mpd/w(1,n)
    ELSE ! singular
        w(1,n)=0.0_mpd
    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)

db3dec()

subroutine db3dec	(	real(mpd), dimension(3,n), intent(inout)	w,
		integer(mpi), intent(inout)	n,
		real(mpd), dimension(n), intent(out)	aux
	)

Decomposition (M=2).

W is a symmetrix positive definite band matrix of bandwidth 2(+1). W(1,.) are the diagonal elements, W(2,.) and W(3,.) are the next diagonals; W(3,N-1), W(2,N) and W(3,N) are never referenced. AUX is an auxiliary array of length N. W is decomposed to L D Lt, where D = diagonal and L unit triangular. A row is set to zero, if the diagonal element is reduced in previous steps by a word length (i.e. global correlation coefficient large). The resulting L and D replace W: the diagonal elements W(1,...) will contain the inverse of the D-elements; the diagonal elements of L are all 1 and not stored. The other elements of L are stored in the corresponding elements of W.

ENTRY DB3SLV(W,N,B, X), solution B -> X
ENTRY DB3IEL(W,N, V), V = inverse band matrix elements

Parameters

[in,out]	W	symmetric band matrix
[in]	N	size
[in]	AUX	scratch array

Definition at line 710 of file Dbandmatrix.f90.

    USE mpdef
 
    implicit none
    INTEGER(mpi) :: i
 
 
    REAL(mpd), INTENT(IN OUT)         :: w(3,n)
    INTEGER(mpi), INTENT(IN OUT)                  :: n
    REAL(mpd), INTENT(OUT)            :: aux(n)
    ! decompos
 
    REAL(mpd) :: v(3,n),b(n),x(n), rxw
 
    DO i=1,n
        aux(i)=16.0_mpd*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.0_mpd/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.0_mpd
            w(2,i)=0.0_mpd
            w(3,i)=0.0_mpd
        END IF
    END DO
    IF(w(1,n-1)+aux(n-1) > aux(n-1)) THEN
        w(1,n-1)=1.0_mpd/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.0_mpd
        w(2,n-1)=0.0_mpd
    END IF
    IF(w(1,n)+aux(n) > aux(n)) THEN
        w(1,n)=1.0_mpd/w(1,n)
    ELSE ! singular
        w(1,n)=0.0_mpd
    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)

dbcdec()

subroutine dbcdec	(	real(mpd), dimension(mp1,n), intent(inout)	w,
		integer(mpi), intent(in)	mp1,
		integer(mpi), intent(in)	n,
		real(mpd), dimension(n), intent(out)	aux
	)

Decomposition of symmetric band matrix.

ENTRY DBCSLV(W,MP1,N,B, X) for solution B -> X
ENTRY DBCIEL(W,MP1,N, V), V = inverse band matrix elements
ENTRY DBCINB(W,MP1,N, VS), VS = band part of inverse symmetric matrix
ENTRY DBCINV(W,MP1,N, VS), V = inverse symmetric matrix

Parameters

[in,out]	W	symmetric band matrix
[in]	MP1	band width (M) + 1
[in]	N	size
[in]	AUX	scratch array

Definition at line 419 of file Dbandmatrix.f90.

    USE mpdef
 
    implicit none
    INTEGER(mpi) :: i
    INTEGER(mpi) :: j
    INTEGER(mpi) :: k
    !     M=MP1-1                                    N*M(M-1) dot operations
 
    REAL(mpd), INTENT(IN OUT)         :: w(mp1,n)
    INTEGER(mpi), INTENT(IN)                      :: mp1
    INTEGER(mpi), INTENT(IN)                      :: n
    REAL(mpd), INTENT(OUT)            :: aux(n)
    ! decompos
    REAL(mpd) :: v(mp1,n),b(n),x(n), vs(*),rxw
    !     ...
    DO i=1,n
        aux(i)=16.0_mpd*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.0_mpd  ! 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

dbcprb()

subroutine dbcprb	(	real(mpd), dimension(mp1,n), intent(inout)	w,
		integer(mpi), intent(in)	mp1,
		integer(mpi), intent(in)	n
	)

Print band matrix.

Parameters

[in]	W	symmetric band matrix
[in]	MP1	band width (M) + 1
[in]	N	size

Definition at line 572 of file Dbandmatrix.f90.

    USE mpdef
 
    implicit none
    INTEGER(mpi) :: i
    INTEGER(mpi) :: j
 
 
    REAL(mpd), INTENT(IN OUT)         :: w(mp1,n)
    INTEGER(mpi), INTENT(IN)                      :: mp1
    INTEGER(mpi), 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)

dbcprv()

subroutine dbcprv	(	real(mpd), dimension(mp1,n), intent(inout)	w,
		integer(mpi), intent(in)	mp1,
		integer(mpi), intent(in)	n,
		real(mpd), dimension(n), intent(out)	b
	)

Print corr band matrix and pars.

Parameters

[in]	W	symmetric band matrix
[in]	MP1	band width (M) + 1
[in]	N	size
[in]	B	vector

Definition at line 526 of file Dbandmatrix.f90.

    USE mpdef
 
    implicit none
    INTEGER(mpi) :: i
    INTEGER(mpi) :: j
    INTEGER(mpi) :: nj
    REAL(mps) :: rho
 
 
    REAL(mpd), INTENT(IN OUT)         :: w(mp1,n)
    INTEGER(mpi), INTENT(IN)                      :: mp1
    INTEGER(mpi), INTENT(IN)                      :: n
    REAL(mpd), INTENT(OUT)            :: b(n)
 
    REAL(mpd) :: ERR
    INTEGER(mpi) :: irho(5)
    !     ...
    WRITE(*,*) ' '
    WRITE(*,101)
 
    DO i=1,n
        err=sqrt(w(1,i))
        nj=0
        DO j=2,mp1
            IF(i+1-j > 0.AND.nj < 5) THEN
                nj=nj+1
                rho=real(w(j,i+1-j)/(err*sqrt(w(1,i+1-j))),mps)
                irho(nj)=nint(100.0*abs(rho),mpi)
                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)')

dbfdec()

subroutine dbfdec	(	real(mpd), dimension(mp1,n), intent(out)	w,
		integer(mpi), intent(inout)	mp1,
		integer(mpi), intent(in)	n
	)

Decomposition of symmetric band matrix.

Band matrix modified Cholesky decomposition, Philip E.Gill, Walter Murray and Margarete H.Wright: Practical Optimization, Academic Press, 1981

Parameters

[in,out]	W	symmetric band matrix
[in]	MP1	band width (M) + 1
[in]	N	size

Definition at line 904 of file Dbandmatrix.f90.

    USE mpdef
 
    IMPLICIT NONE
    REAL(mpd), INTENT(OUT)            :: w(mp1,n)
    INTEGER(mpi), INTENT(IN OUT)                  :: mp1
    INTEGER(mpi), INTENT(IN)                      :: n
    INTEGER(mpi) :: i,j,k
    REAL(mpd) :: epsm,gamm,xchi,beta,delta,theta
 
    epsm=epsilon(epsm) ! machine precision
    gamm=0.0_mpd   ! max diagonal element
    xchi=0.0_mpd   ! 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.0_mpd,sqrt(real(n*n-1,mpd))),epsm))
    delta=epsm*max(1.0_mpd,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.0_mpd
        DO j=2,min(mp1,n-k+1)
            theta=max(theta,abs(w(j,k)))
        END DO
        w(1,k)=1.0_mpd/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
 

dcfdec()

subroutine dcfdec	(	real(mpd), dimension(*), intent(out)	w,
		integer(mpi), intent(in)	n
	)

Decomposition of symmetric matrix.

Modified Cholesky decomposition, Philip E.Gill, Walter Murray and Margarete H.Wright: Practical Optimization, Academic Press, 1981

Parameters

[in,out]	W	symmetirc matrix
[in]	N	size

Definition at line 852 of file Dbandmatrix.f90.

    USE mpdef
 
    IMPLICIT NONE
    REAL(mpd), INTENT(OUT)            :: w(*)
    INTEGER(mpi), INTENT(IN)                      :: n
    INTEGER(mpi) :: i,j,k
    REAL(mpd) :: epsm,gamm,xchi,beta,delta,theta
 
    epsm=epsilon(epsm) ! machine precision
    gamm=0.0_mpd   ! max diagonal element
    xchi=0.0_mpd   ! 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.0_mpd,sqrt(real(n*n-1,mpd))),epsm))
    delta=epsm*max(1.0_mpd,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.0_mpd
        DO j=k+1,n
            theta=max(theta,abs(w((j*j-j)/2+k)))
        END DO
        w((k*k+k)/2)=1.0_mpd/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
 

dchdec()

subroutine dchdec	(	real(mpd), dimension(*), intent(inout)	w,
		integer(mpi), intent(in)	n,
		real(mpd), dimension(n), intent(out)	aux
	)

Decomposition of symmetric matrix.

ENTRY DCHSLV(W,N,B, X) for solution B -> X
ENTRY DCHINV(W,N,V) for inversion

Parameters

[in,out]	W	symmetirc matrix
[in]	N	size
[in]	AUX	scratch array

Definition at line 245 of file Dbandmatrix.f90.

    USE mpdef
 
    implicit none
    INTEGER(mpi) :: i
    INTEGER(mpi) :: ii
    INTEGER(mpi) :: j
    INTEGER(mpi) :: jj
    INTEGER(mpi) :: k
    INTEGER(mpi) :: kk
    INTEGER(mpi) :: l
    INTEGER(mpi) :: m
 
 
    REAL(mpd), INTENT(IN OUT)         :: w(*)
    INTEGER(mpi), INTENT(IN)                      :: n
    REAL(mpd), INTENT(OUT)            :: aux(n)
 
    REAL(mpd) :: b(*),x(*),v(*),suma,ratio
    !     ...
    DO i=1,n
        aux(i)=16.0_mpd*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.0_mpd/w(ii)                ! (I,I) div !
        ELSE
            w(ii)=0.0_mpd
        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.0_mpd
        END DO
        ii=ii-i+1
    END DO