Symmetric (band) matrix routines. More...

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]	W	symmetric matrix
[in]	N	size
[in]	AUX	scratch 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]	W	symmetric band matrix
[in]	N	size
[in]	AUX	scratch 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]	N	size
[out]	V	inverse 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]	N	size
[in]	B	r.h.s.
[out]	X	solution

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]	W	symmetric band matrix
[in]	N	size
[in]	AUX	scratch 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]	N	size
[out]	V	inverse 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]	N	size
[in]	B	r.h.s.
[out]	X	solution

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]	W	symmetric band matrix, replaced by decomposition
[in]	MP1	band width (M) + 1
[in]	N	size
[in]	AUX	scratch 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]	MP1	band width (M) + 1
[in]	N	size
[out]	V	inverse 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]	MP1	band width (M) + 1
[in]	N	size
[out]	VS	band 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]	MP1	band width (M) + 1
[in]	N	size
[out]	VS	inverse 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]	W	symmetric band matrix
[in]	MP1	band width (M) + 1
[in]	N	size

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]	W	symmetric band matrix
[in]	MP1	band width (M) + 1
[in]	N	size
[in]	B	vector

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]	MP1	band width (M) + 1
[in]	N	size
[in]	B	r.h.s.
[out]	X	solution

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]	W	symmetric band matrix
[in]	MP1	band width (M) + 1
[in]	N	size

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]	W	symmetric matrix
[in]	N	size

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]	W	symmetric matrix, replaced by decomposition
[in]	N	size
[in]	AUX	scratch 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]	N	size
[out]	V	inverse 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]	N	size
[in]	B	r.h.s.
[out]	X	solution

Definition at line 295 of file Dbandmatrix.f90.

Functions/Subroutines

Detailed Description

Function/Subroutine Documentation

◆ condes()

◆ db2dec()

◆ db2iel()

◆ db2slv()

◆ db3dec()

◆ db3iel()

◆ db3slv()

◆ dbcdec()

◆ dbciel()

◆ dbcinb()

◆ dbcinv()

◆ dbcprb()

◆ dbcprv()

◆ dbcslv()

◆ dbfdec()

◆ dcfdec()

◆ dchdec()

◆ dchinv()

◆ dchslv()