Millepede-II V04-17-04
Functions/Subroutines
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. More...
 
subroutine dchslv (w, n, b, x)
 solution B -> X More...
 
subroutine dchinv (w, n, v)
 inversion More...
 
real(mps) function condes (w, n, aux)
 Etimate condition. More...
 
subroutine dbcdec (w, mp1, n, aux)
 Decomposition of symmetric band matrix. More...
 
subroutine dbcslv (w, mp1, n, b, x)
 solution B -> X More...
 
subroutine dbciel (w, mp1, n, v)
 V = inverse band matrix elements. More...
 
subroutine dbcinb (w, mp1, n, vs)
 VS = band part of inverse symmetric matrix. More...
 
subroutine dbcinv (w, mp1, n, vs)
 VS = inverse symmetric matrix. More...
 
subroutine dbcprv (w, mp1, n, b)
 Print corr band matrix and pars. More...
 
subroutine dbcprb (w, mp1, n)
 Print band matrix. More...
 
subroutine db2dec (w, n, aux)
 Decomposition (M=1). More...
 
subroutine db2slv (w, n, b, x)
 solution B -> X More...
 
subroutine db2iel (w, n, v)
 V = inverse band matrix elements. More...
 
subroutine db3dec (w, n, aux)
 Decomposition (M=2). More...
 
subroutine db3slv (w, n, b, x)
 solution B -> X More...
 
subroutine db3iel (w, n, v)
 V = inverse band matrix elements. More...
 
subroutine dcfdec (w, n)
 Decomposition of symmetric matrix. More...
 
subroutine dbfdec (w, mp1, n)
 Decomposition of symmetric band matrix. More...
 

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 - 2022 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
    SUBROUTINE DCHSLV(W,N,B, X)      ! solution B -> X
    SUBROUTINE 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
    SUBROUTINE DBCSLV(W,MP1,N,B, X)  ! solution B -> X
    SUBROUTINE DBCIEL(W,MP1,N, V)    ! V = inverse band matrix elements
    SUBROUTINE 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)
    SUBROUTINE DB2SLV(W,N,B, X)      ! solution B -> X
    SUBROUTINE DB2IEL(W,N, V)        ! V = inverse band matrix elements

    SUBROUTINE DB3DEC(W,N, AUX)      ! decomposition (M=2)
    SUBROUTINE DB3SLV(W,N,B, X)      ! solution B -> X
    SUBROUTINE 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((n*n+n)/2), intent(in)  w,
integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(inout)  aux 
)

Etimate condition.

Parameters
[in]Wsymmetric matrix
[in]Nsize
[in]AUXscratch array
Returns
condition

Definition at line 384 of file Dbandmatrix.f90.

References condes(), and mpdef::mps.

Referenced by condes().

◆ 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.

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

Parameters
[in,out]Wsymmetric band matrix
[in]Nsize
[in]AUXscratch array

Definition at line 743 of file Dbandmatrix.f90.

◆ db2iel()

subroutine db2iel ( real(mpd), dimension(2,n), intent(in)  w,
integer(mpi), intent(in)  n,
real(mpd), dimension(2,n), intent(out)  v 
)

V = inverse band matrix elements.

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]Nsize
[out]Vinverse band matrix elements

Definition at line 815 of file Dbandmatrix.f90.

◆ db2slv()

subroutine db2slv ( real(mpd), dimension(2,n), intent(in)  w,
integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  b,
real(mpd), dimension(n), intent(out)  x 
)

solution B -> X

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]Nsize
[in]Br.h.s.
[out]Xsolution

Definition at line 784 of file Dbandmatrix.f90.

◆ 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.

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

Parameters
[in,out]Wsymmetric band matrix
[in]Nsize
[in]AUXscratch array

Definition at line 865 of file Dbandmatrix.f90.

◆ db3iel()

subroutine db3iel ( real(mpd), dimension(3,n), intent(in)  w,
integer(mpi), intent(in)  n,
real(mpd), dimension(3,n), intent(out)  v 
)

V = inverse band matrix elements.

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]Nsize
[out]Vinverse band matrix elements

Definition at line 953 of file Dbandmatrix.f90.

◆ db3slv()

subroutine db3slv ( real(mpd), dimension(3,n), intent(in)  w,
integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  b,
real(mpd), dimension(n), intent(out)  x 
)

solution B -> X

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]Nsize
[in]Br.h.s.
[out]Xsolution

Definition at line 920 of file Dbandmatrix.f90.

◆ 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.

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

Parameters
[in,out]Wsymmetric band matrix, replaced by decomposition
[in]MP1band width (M) + 1
[in]Nsize
[in]AUXscratch array

decomposition, bandwidth M

Definition at line 459 of file Dbandmatrix.f90.

Referenced by sqmibb(), and sqmibb2().

◆ dbciel()

subroutine dbciel ( real(mpd), dimension(mp1,n), intent(in)  w,
integer(mpi), intent(in)  mp1,
integer(mpi), intent(in)  n,
real(mpd), dimension(mp1,n), intent(out)  v 
)

V = inverse band matrix elements.

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]MP1band width (M) + 1
[in]Nsize
[out]Vinverse band matrix elements

Definition at line 542 of file Dbandmatrix.f90.

◆ dbcinb()

subroutine dbcinb ( real(mpd), dimension(mp1,n), intent(in)  w,
integer(mpi), intent(in)  mp1,
integer(mpi), intent(in)  n,
real(mpd), dimension((n*n+n)/2), intent(out)  vs 
)

VS = band part of inverse symmetric matrix.

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]MP1band width (M) + 1
[in]Nsize
[out]VSband part of inverse symmetric matrix

Definition at line 576 of file Dbandmatrix.f90.

Referenced by sqmibb(), and sqmibb2().

◆ dbcinv()

subroutine dbcinv ( real(mpd), dimension(mp1,n), intent(in)  w,
integer(mpi), intent(in)  mp1,
integer(mpi), intent(in)  n,
real(mpd), dimension((n*n+n)/2), intent(out)  vs 
)

VS = inverse symmetric matrix.

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]MP1band width (M) + 1
[in]Nsize
[out]VSinverse symmetric matrix

Definition at line 613 of file Dbandmatrix.f90.

Referenced by sqmibb(), and sqmibb2().

◆ 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]Wsymmetric band matrix
[in]MP1band width (M) + 1
[in]Nsize

Definition at line 693 of file Dbandmatrix.f90.

◆ 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]Wsymmetric band matrix
[in]MP1band width (M) + 1
[in]Nsize
[in]Bvector

Definition at line 647 of file Dbandmatrix.f90.

References mpdef::mpi.

◆ dbcslv()

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

solution B -> X

Parameters
[in,out]W(decomposition of) )symmetric band matrix
[in]MP1band width (M) + 1
[in]Nsize
[in]Br.h.s.
[out]Xsolution

Definition at line 503 of file Dbandmatrix.f90.

Referenced by sqmibb(), and sqmibb2().

◆ 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]Wsymmetric band matrix
[in]MP1band width (M) + 1
[in]Nsize

Definition at line 1092 of file Dbandmatrix.f90.

◆ dcfdec()

subroutine dcfdec ( real(mpd), dimension((n*n+n)/2), 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]Wsymmetric matrix
[in]Nsize

Definition at line 1040 of file Dbandmatrix.f90.

◆ dchdec()

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

Decomposition of symmetric matrix.

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

Parameters
[in,out]Wsymmetric matrix, replaced by decomposition
[in]Nsize
[in]AUXscratch array

Definition at line 245 of file Dbandmatrix.f90.

◆ dchinv()

subroutine dchinv ( real(mpd), dimension((n*n+n)/2), intent(in)  w,
integer(mpi), intent(in)  n,
real(mpd), dimension((n*n+n)/2), intent(out)  v 
)

inversion

Parameters
[in,out]W(decomposition of) symmetric matrix
[in]Nsize
[out]Vinverse symmetric matrix

Definition at line 345 of file Dbandmatrix.f90.

◆ dchslv()

subroutine dchslv ( real(mpd), dimension((n*n+n)/2), intent(in)  w,
integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  b,
real(mpd), dimension(n), intent(out)  x 
)

solution B -> X

Parameters
[in,out]W(decomposition of) symmetric matrix
[in]Nsize
[in]Br.h.s.
[out]Xsolution

Definition at line 295 of file Dbandmatrix.f90.

Generated on Thu Dec 5 2024 14:18:59 for Millepede-II by doxygen 1.9.3