Millepede-II V04-17-03
Functions/Subroutines
mpnum.f90 File Reference

General linear algebra routines. More...

Go to the source code of this file.

Functions/Subroutines

subroutine sqminv (v, b, n, nrank, diag, next)
 Matrix inversion and solution. More...
 
subroutine sqminl (v, b, n, nrank, diag, next, vk, mon)
 Matrix inversion for LARGE matrices. More...
 
subroutine devrot (n, diag, u, v, work, iwork)
 Diagonalization. More...
 
subroutine devsig (n, diag, u, b, coef)
 Calculate significances. More...
 
subroutine devsol (n, diag, u, b, x, work)
 Solution by diagonalization. More...
 
subroutine devinv (n, diag, u, v)
 Inversion by diagonalization. More...
 
subroutine choldc (g, n)
 Cholesky decomposition. More...
 
subroutine cholsl (g, x, n)
 Solution after decomposition. More...
 
subroutine cholin (g, v, n)
 Inversion after decomposition. More...
 
subroutine chdec2 (g, n, nrank, evmax, evmin, mon)
 Cholesky decomposition (LARGE pos. More...
 
subroutine chslv2 (g, x, n)
 Solve A*x=b using Cholesky decomposition. More...
 
subroutine vabdec (n, val, ilptr)
 Variable band matrix decomposition. More...
 
subroutine vabmmm (n, val, ilptr)
 Variable band matrix print minimum and maximum. More...
 
subroutine vabslv (n, val, ilptr, x)
 Variable band matrix solution. More...
 
real(mpd) function dbdot (n, dx, dy)
 Dot product. More...
 
subroutine dbaxpy (n, da, dx, dy)
 Multiply, addition. More...
 
subroutine dbsvx (v, a, b, n)
 Product symmetric matrix, vector. More...
 
subroutine dbsvxl (v, a, b, n)
 Product LARGE symmetric matrix, vector. More...
 
subroutine dbgax (a, x, y, m, n)
 Multiply general M-by-N matrix A and N-vector X. More...
 
subroutine dbavat (v, a, w, n, m, iopt)
 A V AT product (similarity). More...
 
subroutine dbavats (v, a, is, w, n, m, iopt, sc)
 A V AT product (similarity, sparse). More...
 
subroutine dbmprv (lun, x, v, n)
 Print symmetric matrix, vector. More...
 
subroutine dbprv (lun, v, n)
 Print symmetric matrix. More...
 
subroutine heapf (a, n)
 Heap sort direct (real). More...
 
subroutine sort1k (a, n)
 Quick sort 1. More...
 
subroutine sort2k (a, n)
 Quick sort 2. More...
 
subroutine sort2i (a, n)
 Quick sort 2 with index. More...
 
subroutine sort22l (a, b, n)
 Quick sort 2 with index. More...
 
real(mps) function chindl (n, nd)
 Chi2/ndf cuts. More...
 
subroutine lltdec (n, c, india, nrkd, iopt)
 LLT decomposition. More...
 
subroutine lltfwd (n, c, india, x)
 Forward solution. More...
 
subroutine lltfwds (n, c, india, x, i0, ns)
 Forward solution (sparse). More...
 
subroutine lltbwd (n, c, india, x)
 Backward solution. More...
 
subroutine equdec (n, m, ls, c, india, nrkd, nrkd2)
 Decomposition of equilibrium systems. More...
 
subroutine equslv (n, m, c, india, x)
 Solution of equilibrium systems (after decomposition). More...
 
subroutine equdecs (n, m, b, nm, ls, c, india, l, nrkd, nrkd2)
 Decomposition of (sparse) equilibrium systems. More...
 
subroutine equslvs (n, m, b, nm, c, india, l, x)
 Solution of (sparse) equilibrium systems (after decomposition). More...
 
subroutine precon (p, n, c, cu, a, s, nrkd)
 Constrained preconditioner, decomposition. More...
 
subroutine presol (p, n, cu, a, s, x, y)
 Constrained preconditioner, solution. More...
 
subroutine precons (p, n, b, nm, c, cu, a, l, s, nrkd)
 Constrained (sparse) preconditioner, decomposition. More...
 
subroutine presols (p, n, b, nm, cu, a, l, s, x, y)
 Constrained (sparse) preconditioner, solution. More...
 
subroutine sqmibb (v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag, evdmin, evdmax)
 Bordered band matrix. More...
 
subroutine sqmibb2 (v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag, evdmin, evdmax)
 Band bordered matrix. More...
 

Detailed Description

General linear algebra routines.

Author
Volker Blobel, University Hamburg, 2005-2009 (initial Fortran77 version)
Claus Kleinwort, DESY (maintenance and developement)

***** Collection of utility routines from V. Blobel *****

V. Blobel, Univ. Hamburg
Numerical subprograms used in MP-II: matrix equations,
   and matrix products, double precision

Solution by inversion
   SQMINV
   SQMINL  for LARGE matrix, use OpenMP (CHK)

Solution by diagonalization
   DEVROT, DEVPRT, DEFSOL,DEVINV

Solution by Cholesky decomposition of symmetric matrix
   CHOLDC
   CHDEC2, CHSLV2 for large (positive definite) matrix, use OpenMP (CHK)

Solution by Cholesky decomposition of variable-band matrix
   VABDEC

Solution by Cholesky decomposition of bordered band matrix
   SQMIBB    upper/left  border (CHK)
   SQMIBB2   lower/right border (CHK)

Matrix/vector products
   DBDOT     dot vector product
   DBAXPY    multiplication and addition
   DBSVX     symmetric matrix vector
   DBSVX     LARGE symmetric matrix vector (CHK)
   DBGAX     general matrix vector
   DBAVAT    AVAT product
   DBAVATS   AVAT product for sparse A (CHK)
   DBMPRV    print parameter and matrix
   DBPRV     print matrix  (CHK)

Chi^2 cut values
   CHINDL

Accurate summation (moved to pede.f90)
   ADDSUM

Sorting
   HEAPF    heap sort reals direct
   SORT1K   sort 1-dim key-array (CHK)
   SORT2K   sort 2-dim key-array
   SORT2I   sort 2-dim key-array with index (CHK)
   SORT22L  sort 2-dim key-array with two additional values (CHK)
                 and one additional "long" value   

Definition in file mpnum.f90.

Function/Subroutine Documentation

◆ chdec2()

subroutine chdec2 ( real(mpd), dimension((int(n,mpl)*int(n,mpl)+int(n,mpl))/2), intent(inout)  g,
integer(mpi), intent(in)  n,
integer(mpi), intent(out)  nrank,
real(mpd), intent(out)  evmax,
real(mpd), intent(out)  evmin,
integer(mpi), intent(in)  mon 
)

Cholesky decomposition (LARGE pos.

def. matrices).

Cholesky decomposition of the matrix G: G = L D L^T

  • G = symmetric matrix, in symmetric storage mode, positive definite
  • L = unit (upper!) triangular matrix (1's on diagonal)
  • D = diagonal matrix (elements store on diagonal of L)

The sqrts of the usual Cholesky decomposition are avoided by D. Matrices L and D are stored in the place of matrix G; after the decomposition, the solution is done by CHSLV2.

Parameters
[in,out]gsymmetric matrix, replaced by D,L
[in]nsize of matrix
[out]NRANKrank of matrix g
[out]EVMAXlargest element in D
[out]EVMINsmallest element in D
[in]MONflag for progress monitoring

Definition at line 891 of file mpnum.f90.

References monpgs().

Referenced by mchdec().

◆ chindl()

real(mps) function chindl ( integer(mpi), intent(in)  n,
integer(mpi), intent(in)  nd 
)

Chi2/ndf cuts.

Return limit in Chi^2/ndf for N sigmas (N=1, 2 or 3).

Parameters
[in]nnumber of sigmas
[in]ndndf
Returns
Chi2/ndf cut value

Definition at line 2078 of file mpnum.f90.

References chindl(), and mpdef::mps.

Referenced by chindl(), fitloc(), and initgl().

◆ choldc()

subroutine choldc ( real(mpd), dimension((n*n+n)/2), intent(inout)  g,
integer(mpi), intent(in)  n 
)

Cholesky decomposition.

Cholesky decomposition of the matrix G: G = L D L^T

  • G = symmetric matrix, in symmetric storage mode
  • L = unit triangular matrix (1's on diagonal)
  • D = diagonal matrix (elements store on diagonal of L)

The sqrts of the usual Cholesky decomposition are avoided by D. Matrices L and D are stored in the place of matrix G; after the decomposition, the solution of matrix equations and the computation of the inverse of the (original) matrix G are done by CHOLSL and CHOLIN.

Parameters
[in,out]gsymmetric matrix, replaced by D,L
[in]nsize of matrix

Definition at line 745 of file mpnum.f90.

◆ cholin()

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

Inversion after decomposition.

The inverse of the (original) matrix G is computed and stored in symmetric storage mode in matrix V. Arrays G and V must be different arrays.

Parameters
[in]gdecomposed symmetric matrix
[in,out]vinverse matrix
[in]nsize of matrix

Definition at line 837 of file mpnum.f90.

◆ cholsl()

subroutine cholsl ( real(mpd), dimension((n*n+n)/2), intent(in)  g,
real(mpd), dimension(n), intent(inout)  x,
integer(mpi), intent(in)  n 
)

Solution after decomposition.

The matrix equation G X = B is solved for X, where the matrix G in the argument is already decomposed by CHOLDC. The vector B is called X in the argument and the content is replaced by the resulting vector X.

Parameters
[in]gdecomposed symmetric matrix
[in,out]xr.h.s vector B, replaced by solution vector X
[in]nsize of matrix

Definition at line 791 of file mpnum.f90.

◆ chslv2()

subroutine chslv2 ( real(mpd), dimension((int(n,mpl)*int(n,mpl)+int(n,mpl))/2), intent(in)  g,
real(mpd), dimension(n), intent(inout)  x,
integer(mpi), intent(in)  n 
)

Solve A*x=b using Cholesky decomposition.

Backward, forward substitution.

Parameters
[in]gdecomposed symmetric matrix
[in,out]xrhs/solution
[in]nsize of matrix

Definition at line 953 of file mpnum.f90.

Referenced by mchdec().

◆ dbavat()

subroutine dbavat ( real(mpd), dimension((n*n+n)/2), intent(in)  v,
real(mpd), dimension(n*m), intent(in)  a,
real(mpd), dimension((m*m+m)/2), intent(inout)  w,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m,
integer(mpi), intent(in)  iopt 
)

A V AT product (similarity).

Multiply symmetric N-by-N matrix from the left with general M-by-N matrix and from the right with the transposed of the same general matrix to form symmetric M-by-M matrix (used for error propagation).

Parameters
[in]Vsymmetric N-by-N matrix
[in]Ageneral M-by-N matrix
[in,out]Wsymmetric M-by-M matrix
[in]Ncolumns of A
[in]Mrows of A
[in]iopt(<>0: don't reset W)

Definition at line 1389 of file mpnum.f90.

Referenced by loopbf().

◆ dbavats()

subroutine dbavats ( real(mpd), dimension((n*n+n)/2), intent(in)  v,
real(mpd), dimension(n*m), intent(in)  a,
integer(mpi), dimension(2*n*m+n+m+1), intent(in)  is,
real(mpd), dimension((m*m+m)/2), intent(inout)  w,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m,
integer(mpi), intent(in)  iopt,
integer(mpi), dimension(n), intent(out)  sc 
)

A V AT product (similarity, sparse).

Multiply symmetric N-by-N matrix from the left with sparse M-by-N matrix and from the right with the transposed of the same general matrix to form symmetric M-by-M matrix (used for error propagation).

Parameters
[in]Vsymmetric N-by-N matrix
[in]Asparse M-by-N matrix, content
[in]ISsparse M-by-N matrix, structure
[in,out]Wsymmetric M-by-M matrix
[in]Ncolumns of A
[in]Mrows of A
[in]iopt(<>0: don't reset W)
[in]SCscratch array

Sparsity structure:

  • IS(1..M): row offsets
  • IS(M+1..N+M+1): column offsets
  • IS(IS(1)+1..IS(M+1)): non-zero columns (column number, index for A)
  • IS(IS(M+1)+1..IS(M+N+1)): non-zero rows (row number, index for A)

Definition at line 1470 of file mpnum.f90.

Referenced by loopbf().

◆ dbaxpy()

subroutine dbaxpy ( integer(mpi), intent(in)  n,
real(mpd), intent(in)  da,
real(mpd), dimension(n), intent(in)  dx,
real(mpd), dimension(n), intent(inout)  dy 
)

Multiply, addition.

Constant times vector added to a vector: DY:=DY+DA*DX

Parameters
[in]nvector size
[in]dxvector
[in,out]dyvector
[in]dascalar

Definition at line 1233 of file mpnum.f90.

◆ dbdot()

real(mpd) function dbdot ( integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  dx,
real(mpd), dimension(n), intent(in)  dy 
)

Dot product.

Parameters
[in]nvector size
[in]dxvector
[in]dyvector
Returns
dot product dx*dy

Definition at line 1202 of file mpnum.f90.

References dbdot().

Referenced by dbdot().

◆ dbgax()

subroutine dbgax ( real(mpd), dimension(n*m), intent(in)  a,
real(mpd), dimension(n), intent(in)  x,
real(mpd), dimension(m), intent(out)  y,
integer(mpi), intent(in)  m,
integer(mpi), intent(in)  n 
)

Multiply general M-by-N matrix A and N-vector X.

Parameters
[in]Ageneral M-by-N matrix (A11 A12 ... A1N A21 A22 ...)
[in]XN vector
[out]Y= M vector
[in]Mrows of A
[in]Ncolumns of A

Definition at line 1351 of file mpnum.f90.

◆ dbmprv()

subroutine dbmprv ( integer(mpi), intent(in)  lun,
real(mpd), dimension(n), intent(in)  x,
real(mpd), dimension((n*n+n)/2), intent(in)  v,
integer(mpi), intent(in)  n 
)

Print symmetric matrix, vector.

Prints the n-vector X and the symmetric N-by-N covariance matrix V, the latter as a correlation matrix.

Parameters
[in]lununit number
[in]xvector
[in]vsymmetric matrix
[in]nsize of matrix, vector

Definition at line 1540 of file mpnum.f90.

References mpdef::mpi.

◆ dbprv()

subroutine dbprv ( integer(mpi), intent(in)  lun,
real(mpd), dimension((n*n+n)/2), intent(in)  v,
integer(mpi), intent(in)  n 
)

Print symmetric matrix.

Prints the symmetric N-by-N matrix V.

Parameters
[in]lununit number
[in]vsymmetric matrix
[in]nsize of matrix,

Definition at line 1615 of file mpnum.f90.

◆ dbsvx()

subroutine dbsvx ( real(mpd), dimension((n*n+n)/2), intent(in)  v,
real(mpd), dimension(n), intent(in)  a,
real(mpd), dimension(n), intent(out)  b,
integer(mpi), intent(in)  n 
)

Product symmetric matrix, vector.

Multiply symmetric N-by-N matrix and N-vector.

Parameters
[in]vsymmetric matrix
[in]avector
[out]bvector B = V * A
[in]nsize of matrix

Definition at line 1264 of file mpnum.f90.

Referenced by feasib().

◆ dbsvxl()

subroutine dbsvxl ( real(mpd), dimension((int(n,mpl)*int(n,mpl)+int(n,mpl))/2), intent(in)  v,
real(mpd), dimension(n), intent(in)  a,
real(mpd), dimension(n), intent(out)  b,
integer(mpi), intent(in)  n 
)

Product LARGE symmetric matrix, vector.

Multiply LARGE symmetric N-by-N matrix and N-vector:

Parameters
[in]vsymmetric matrix
[in]avector
[out]bproduct vector B = V * A
[in]nsize of matrix

Definition at line 1308 of file mpnum.f90.

Referenced by minver().

◆ devinv()

subroutine devinv ( integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  diag,
real(mpd), dimension(n,n), intent(in)  u,
real(mpd), dimension((n*n+n)/2), intent(out)  v 
)

Inversion by diagonalization.

Get inverse matrix V from DIAG and U.

Parameters
[in]Nsize of matrix
[in]DIAGdiagonal elements
[in]Utransformation matrix
[out]Vsymmmetric matrix

Definition at line 696 of file mpnum.f90.

Referenced by zdiags().

◆ devrot()

subroutine devrot ( integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(out)  diag,
real(mpd), dimension(n,n), intent(out)  u,
real(mpd), dimension((n*n+n)/2), intent(in)  v,
real(mpd), dimension(n), intent(out)  work,
integer(mpi), dimension(n), intent(out)  iwork 
)

Diagonalization.

Determination of eigenvalues and eigenvectors of symmetric matrix V by Householder method

Parameters
[in]nsize of matrix
[out]diagdiagonal elements
[out]utransformation matrix
[in]vsymmetric matrix, unchanged
[out]workwork array
[out]iworkwork array

Definition at line 369 of file mpnum.f90.

References peend().

Referenced by mdiags().

◆ devsig()

subroutine devsig ( integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  diag,
real(mpd), dimension(n,n), intent(in)  u,
real(mpd), dimension(n), intent(in)  b,
real(mpd), dimension(n), intent(out)  coef 
)

Calculate significances.

Definition at line 611 of file mpnum.f90.

Referenced by mdiags().

◆ devsol()

subroutine devsol ( integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  diag,
real(mpd), dimension(n,n), intent(in)  u,
real(mpd), dimension(n), intent(in)  b,
real(mpd), dimension(n), intent(out)  x,
real(mpd), dimension(n), intent(out)  work 
)

Solution by diagonalization.

Solution of matrix equation V * X = B after diagonalization of V.

Parameters
[in]Nsize of matrix
[in]DIAGdiagonal elements
[in]Utransformation matrix
[in]Br.h.s. of matrix equation (unchanged)
[out]Xsolution vector
[out]WORKwork array

Definition at line 649 of file mpnum.f90.

Referenced by mdiags().

◆ equdec()

subroutine equdec ( integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m,
integer(mpi), intent(in)  ls,
real(mpd), dimension(india(n)+n*m+(m*m+m)/2), intent(inout)  c,
integer(mpi), dimension(n+m), intent(inout)  india,
integer(mpi), intent(out)  nrkd,
integer(mpi), intent(out)  nrkd2 
)

Decomposition of equilibrium systems.

N x N matrix C:     starting with sym.pos.def. matrix (N)
length  of array C: INDIA(N) + N*M + (M*M+M)/2
Content of array C: band matrix, as described by INDIA(1)...INDIA(N)
       followed by: NxM elements of constraint matrix A
       followed by: (M*M+M)/2 unused elements
                    INDIA(N+1)...INDIA(N+M) defined internally
Parameters
[in]nsize of symmetric matrix
[in]mnumber of constrains
[in]lsflag for skyline decomposition
[in,out]ccombined variable-band + constraints matrix, replaced by decomposition
[in,out]indiapointer array
[out]nrkdremoved components
[out]nrkd2removed components

Definition at line 2355 of file mpnum.f90.

References lltdec(), lltfwd(), and mpdef::mpl.

◆ equdecs()

subroutine equdecs ( integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m,
integer(mpi), intent(in)  b,
integer(mpl), intent(in)  nm,
integer(mpi), intent(in)  ls,
real(mpd), dimension(nm+india(n)+(m*m+m)/2), intent(inout)  c,
integer(mpi), dimension(n+m), intent(inout)  india,
integer(mpi), dimension(4*b), intent(in)  l,
integer(mpi), intent(out)  nrkd,
integer(mpi), intent(out)  nrkd2 
)

Decomposition of (sparse) equilibrium systems.

CHK, June 2021

N x N matrix C:     starting with sym.pos.def. matrix (N)
length  of array C: INDIA(N) + (M*M+M)/2 + NM
Content of array C: band matrix, as described by INDIA(1)...INDIA(N)
       followed by: (M*M+M)/2 unused elements
       followed by: blocks of constraint matrix A
                    INDIA(N+1)...INDIA(N+M) defined internally
Parameters
[in]nsize of symmetric matrix
[in]mnumber of constrains
[in]bnumber of constraint blocks
[in]nm(integrated) size of constraint blocks (<=N*M)
[in]lsflag for skyline decomposition
[in,out]ccombined variable-band + constraints matrix, replaced by decomposition
[in,out]indiapointer array
[in]lconstraint (block matrix row) offsets
[out]nrkdremoved components
[out]nrkd2removed components

Definition at line 2486 of file mpnum.f90.

References lltdec(), and lltfwds().

Referenced by mminrs(), and mminrsqlp().

◆ equslv()

subroutine equslv ( integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m,
real(mpd), dimension(india(n)+n*m+(m*m+m)/2), intent(inout)  c,
integer(mpi), dimension(n+m), intent(in)  india,
real(mpd), dimension(n+m), intent(inout)  x 
)

Solution of equilibrium systems (after decomposition).

N x N matrix C:     starting with sym.pos.def. matrix (N)
length  of array C: INDIA(N) + N*M + (M*M+M)/2
Content of array C: band matrix, as described by INDIA(1)...INDIA(N)
       followed by: NxM elements of constraint matrix A
       followed by: (M*M+M)/2 unused elements
                    INDIA(N+1)...INDIA(N+M) defined internally
Parameters
[in]nsize of symmetric matrix
[in]mnumber of constrains
[in]cdecomposed combined variable-band + constraints matrix
[in]indiapointer array
[in,out]xr.h.s vector B, replaced by solution vector X

Definition at line 2422 of file mpnum.f90.

References lltbwd(), lltfwd(), and mpdef::mpl.

◆ equslvs()

subroutine equslvs ( integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m,
integer(mpi), intent(in)  b,
integer(mpl), intent(in)  nm,
real(mpd), dimension(nm+india(n)+(m*m+m)/2), intent(inout)  c,
integer(mpi), dimension(n+m), intent(in)  india,
integer(mpi), dimension(4*b), intent(in)  l,
real(mpd), dimension(n+m), intent(inout)  x 
)

Solution of (sparse) equilibrium systems (after decomposition).

CHK, June 2021

N x N matrix C:     starting with sym.pos.def. matrix (N)
length  of array C: INDIA(N) + (M*M+M)/2 + NM
Content of array C: band matrix, as described by INDIA(1)...INDIA(N)
       followed by: NxM elements of constraint matrix A
       followed by: (M*M+M)/2 unused elements
                    INDIA(N+1)...INDIA(N+M) defined internally
Parameters
[in]nsize of symmetric matrix
[in]mnumber of constrains
[in]bnumber of constraint blocks
[in]nm(integrated) size of constraint blocks (<=N*M)
[in]cdecomposed combined variable-band + constraints matrix
[in]indiapointer array
[in]lconstraint (block matrix row) offsets
[in,out]xr.h.s vector B, replaced by solution vector X

Definition at line 2613 of file mpnum.f90.

References lltbwd(), and lltfwd().

Referenced by mvsolv().

◆ heapf()

subroutine heapf ( real(mps), dimension(n), intent(inout)  a,
integer(mpi), intent(in)  n 
)

Heap sort direct (real).

Real keys A(1..N), sorted at return.

Parameters
[in,out]aarray of keys
[in]nnumber of keys

Definition at line 1659 of file mpnum.f90.

Referenced by bintab(), hmpmak(), and rmesig().

◆ lltbwd()

subroutine lltbwd ( integer(mpi), intent(in)  n,
real(mpd), dimension(india(n)), intent(inout)  c,
integer(mpi), dimension(n), intent(in)  india,
real(mpd), dimension(n), intent(inout)  x 
)

Backward solution.

The matrix equation A X = B is solved by forward + backward solution. The matrix is assumed to decomposed before using LLTDEC. The array X(N) contains on entry the right-hand-side B(N); at return it contains the solution.

Parameters
[in]nsize of matrix
[in,out]cdecomposed variable-band matrix
[in]indiapointer array
[in,out]xr.h.s vector B, replaced by solution vector X

Definition at line 2315 of file mpnum.f90.

Referenced by equslv(), and equslvs().

◆ lltdec()

subroutine lltdec ( integer(mpi), intent(in)  n,
real(mpd), dimension(india(n)), intent(inout)  c,
integer(mpi), dimension(n), intent(in)  india,
integer(mpi), intent(out)  nrkd,
integer(mpi), intent(in)  iopt 
)

LLT decomposition.

Decomposition: C = L L^T.

Variable-band matrix row-Doolittle decomposition of pos. def. matrix. A variable-band NxN symmetric matrix, is stored row by row in the array C(.). For each row all coefficients from the first non-zero element in the row to the diagonal is stored. The pointer array INDIA(N) contains the indices in C(.) of the diagonal elements. INDIA(1) is always 1, and INDIA(N) is equal to the total number of coefficients stored, called the profile. The form of a variable-band matrix is preserved in the L D L^T decomposition. No fill-in is created ahead in any row or ahead of the first entry in any column, but existing zero-values will become non-zero. The decomposition is done "in-place". (The diagonal will contain the inverse of the diaginal of L).

  • NRKD = 0 no component removed
  • NRKD < 0 1 component removed, negative index
  • NRKD > 1 number of

The matrix C is assumed to be positive definite, e.g. from the normal equations of least squares. The (positive) diagonal elements are reduced during decomposition. If a diagonal element is reduced by about a word length (see line "test for linear dependence"), then the pivot is assumed as zero and the entire row/column is reset to zero, removing the corresponding element from the solution. Optionally use only diagonal element in this case to preserve rank (changing band to skyline matrix).

Parameters
[in]nsize of matrix
[in,out]cvariable-band matrix, replaced by L
[in]indiapointer array
[out]nrkdremoved components
[in]iopt>0: use diagonal to preserve rank ('skyline')

Definition at line 2155 of file mpnum.f90.

Referenced by equdec(), and equdecs().

◆ lltfwd()

subroutine lltfwd ( integer(mpi), intent(in)  n,
real(mpd), dimension(india(n)), intent(inout)  c,
integer(mpi), dimension(n), intent(in)  india,
real(mpd), dimension(n), intent(inout)  x 
)

Forward solution.

The matrix equation A X = B is solved by forward + backward solution. The matrix is assumed to decomposed before using LLTDEC. The array X(N) contains on entry the right-hand-side B(N); at return it contains the solution.

Parameters
[in]nsize of matrix
[in,out]cdecomposed variable-band matrix
[in]indiapointer array
[in,out]xr.h.s vector B, replaced by solution vector X

Definition at line 2237 of file mpnum.f90.

Referenced by equdec(), equslv(), and equslvs().

◆ lltfwds()

subroutine lltfwds ( integer(mpi), intent(in)  n,
real(mpd), dimension(india(n)), intent(inout)  c,
integer(mpi), dimension(n), intent(in)  india,
real(mpd), dimension(n), intent(inout)  x,
integer(mpi), intent(in)  i0,
integer(mpi), intent(in)  ns 
)

Forward solution (sparse).

CHK, June 2021

The matrix equation A X = B is solved by forward + backward solution. The matrix is assumed to decomposed before using LLTDEC. The array X(N) contains on entry the right-hand-side B(N); at return it contains the solution.

Parameters
[in]nsize of matrix
[in,out]cdecomposed variable-band matrix
[in]indiapointer array
[in,out]xr.h.s vector B, replaced by solution vector X
[in]i0offset of non zero region in x
[in]nssize of non zero region in x

Definition at line 2275 of file mpnum.f90.

Referenced by equdecs().

◆ precon()

subroutine precon ( integer(mpi), intent(in)  p,
integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  c,
real(mpd), dimension(n), intent(out)  cu,
real(mpd), dimension(n,p), intent(inout)  a,
real(mpd), dimension((p*p+p)/2), intent(out)  s,
integer(mpi), intent(out)  nrkd 
)

Constrained preconditioner, decomposition.

Constrained preconditioner, e.g for GMRES solution:

                                           intermediate
   (            )  (   )      (   )           (   )
   (   C    A^T )  ( x )   =  ( y )           ( u )
   (            )  (   )      (   )           (   )
   (   A     0  )  ( l )      ( d )           ( v )

input:
   C(N) is diagonal matrix and remains unchanged
        may be identical to CU(N), then it is changed
   A(N,P) is modified
   Y(N+P) is rhs vector, unchanged
        may be identical to X(N), then it is changed

result:
   CU(N) is 1/sqrt of diagonal matrix C(N)
   X(N+P) is result vector
   S((P*P+P)/2) is Cholesky decomposed symmetric (P,P) matrix
Parameters
[in]pnumber of constraints
[in]nsize of diagonal matrix
[in]cdiagonal matrix (changed if c=cu as actual parameters)
[out]cu1/sqrt(c)
[in,out]aconstraint matrix (size n*p), modified
[out]sCholesky decomposed symmetric (P,P) matrix
[out]nrkdremoved components

Definition at line 2710 of file mpnum.f90.

Referenced by minresmodule::minres().

◆ precons()

subroutine precons ( integer(mpi), intent(in)  p,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  b,
integer(mpl), intent(in)  nm,
real(mpd), dimension(n), intent(in)  c,
real(mpd), dimension(n), intent(out)  cu,
real(mpd), dimension(nm), intent(inout)  a,
integer(mpi), dimension(p), intent(in)  l,
real(mpd), dimension((p*p+p)/2), intent(out)  s,
integer(mpi), intent(out)  nrkd 
)

Constrained (sparse) preconditioner, decomposition.

CHK, June 2021

Constrained preconditioner, e.g for GMRES solution:

                                           intermediate
   (            )  (   )      (   )           (   )
   (   C    A^T )  ( x )   =  ( y )           ( u )
   (            )  (   )      (   )           (   )
   (   A     0  )  ( l )      ( d )           ( v )

input:
   C(N) is diagonal matrix and remains unchanged
        may be identical to CU(N), then it is changed
   A(N,P) is represented by nonzero blocks (one per row), modified
   Y(N+P) is rhs vector, unchanged
        may be identical to X(N), then it is changed

result:
   CU(N) is 1/sqrt of diagonal matrix C(N)
   X(N+P) is result vector
   S((P*P+P)/2) is Cholesky decomposed symmetric (P,P) matrix
Parameters
[in]pnumber of constraints
[in]nsize of diagonal matrix
[in]bnumber of constraint blocks
[in]nm(integrated) size of constraint blocks (<=N*M)
[in]cdiagonal matrix (changed if c=cu as actual parameters)
[out]cu1/sqrt(c)
[in,out]amodified constraint (block) matrix, modified
[in]lconstraint (block matrix row) offsets
[out]sCholesky decomposed symmetric (P,P) matrix
[out]nrkdremoved components

Definition at line 2881 of file mpnum.f90.

Referenced by mminrs(), and mminrsqlp().

◆ presol()

subroutine presol ( integer(mpi), intent(in)  p,
integer(mpi), intent(in)  n,
real(mpd), dimension(n), intent(in)  cu,
real(mpd), dimension(n,p), intent(in)  a,
real(mpd), dimension((p*p+p)/2), intent(in)  s,
real(mpd), dimension(n+p), intent(out)  x,
real(mpd), dimension(n+p), intent(in)  y 
)

Constrained preconditioner, solution.

Parameters
[in]pnumber of constraints
[in]nsize of diagonal matrix
[in]cu1/sqrt(c)
[in]amodified constraint matrix (size n*p)
[in]sCholesky decomposed symmetric (P,P) matrix
[out]xresult vector
[in]yrhs vector (changed if x=y as actual parameters)

Definition at line 2784 of file mpnum.f90.

◆ presols()

subroutine presols ( integer(mpi), intent(in)  p,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  b,
integer(mpl), intent(in)  nm,
real(mpd), dimension(n), intent(in)  cu,
real(mpd), dimension(nm), intent(in)  a,
integer(mpi), dimension(p), intent(in)  l,
real(mpd), dimension((p*p+p)/2), intent(in)  s,
real(mpd), dimension(n+p), intent(out)  x,
real(mpd), dimension(n+p), intent(in)  y 
)

Constrained (sparse) preconditioner, solution.

CHK, June 2021

Parameters
[in]pnumber of constraints
[in]nsize of diagonal matrix
[in]bnumber of constraint blocks
[in]nm(integrated) size of constraint blocks (<=N*M)
[in]cu1/sqrt(c)
[in]amodified constraint (block) matrix
[in]lconstraint (block matrix row) offsets
[in]sCholesky decomposed symmetric (P,P) matrix
[out]xresult vector
[in]yrhs vector (changed if x=y as actual parameters)

Definition at line 2980 of file mpnum.f90.

Referenced by mcsolv().

◆ sort1k()

subroutine sort1k ( integer(mpi), dimension(n), intent(inout)  a,
integer(mpi), intent(in)  n 
)

Quick sort 1.

Quick sort of A(1,N) integer.

Parameters
[in,out]avector of integers, sorted at return
[in]nsize of vector

Definition at line 1714 of file mpnum.f90.

References peend().

Referenced by grpcon(), loop2(), loopbf(), and pepgrp().

◆ sort22l()

subroutine sort22l ( integer(mpi), dimension(4,*), intent(inout)  a,
integer(mpl), dimension(*), intent(inout)  b,
integer(mpi), intent(in)  n 
)

Quick sort 2 with index.

Quick sort of A(4,N) integer.

Parameters
[in,out]avector (quadruplet) of integers, sorted at return and an index
[in,out]bvector of long integers, sorted at return and an index
[in]nsize of vector

Definition at line 1981 of file mpnum.f90.

References peend().

Referenced by upone(), and useone().

◆ sort2i()

subroutine sort2i ( integer(mpi), dimension(3,*), intent(inout)  a,
integer(mpi), intent(in)  n 
)

Quick sort 2 with index.

Quick sort of A(3,N) integer.

Parameters
[in,out]avector (pair) of integers, sorted at return and an index
[in]nsize of vector

Definition at line 1892 of file mpnum.f90.

References peend().

Referenced by grpcon().

◆ sort2k()

subroutine sort2k ( integer(mpi), dimension(2,*), intent(inout)  a,
integer(mpi), intent(in)  n 
)

Quick sort 2.

Quick sort of A(2,N) integer.

Parameters
[in,out]avector (pair) of integers, sorted at return
[in]nsize of vector

Definition at line 1799 of file mpnum.f90.

References peend().

Referenced by peread().

◆ sqmibb()

subroutine sqmibb ( real(mpd), dimension((n*n+n)/2), intent(inout)  v,
real(mpd), dimension(n), intent(out)  b,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  nbdr,
integer(mpi), intent(in)  nbnd,
integer(mpi), intent(in)  inv,
integer(mpi), intent(out)  nrank,
real(mpd), dimension(n*(nbnd+1)), intent(out)  vbnd,
real(mpd), dimension(n*nbdr), intent(out)  vbdr,
real(mpd), dimension(n*nbdr), intent(out)  aux,
real(mpd), dimension((nbdr*nbdr+nbdr)/2), intent(out)  vbk,
real(mpd), dimension(nbdr), intent(out)  vzru,
real(mpd), dimension(nbdr), intent(out)  scdiag,
integer(mpi), dimension(nbdr), intent(out)  scflag,
real(mpd), intent(out)  evdmin,
real(mpd), intent(out)  evdmax 
)

Bordered band matrix.

Obtain solution of a system of linear equations with symmetric bordered band matrix (V * X = B), on request inverse is calculated. For band part root-free Cholesky decomposition and forward/backward substitution is used.

Use decomposition in border and band part for block matrix algebra:

| A  Ct |   | x1 |   | b1 |        , A  is the border part
|       | * |    | = |    |        , Ct is the mixed part
| C  D  |   | x2 |   | b2 |        , D  is the band part

Explicit inversion of D is avoided by using solution X of D*X=C (X=D^-1*C, obtained from Cholesky decomposition and forward/backward substitution)

| x1 |   | E*b1 - E*Xt*b2 |        , E^-1 = A-Ct*D^-1*C = A-Ct*X
|    | = |                |
| x2 |   |  x   - X*x1    |        , x is solution of D*x=b2 (x=D^-1*b2)

Inverse matrix is:

|  E   -E*Xt          |
|                     |            , only band part of (D^-1 + X*E*Xt)
| -X*E  D^-1 + X*E*Xt |              is calculated for inv=1
Parameters
[in,out]vsymmetric N-by-N matrix in symmetric storage mode (V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, ...), replaced by inverse matrix
[in,out]bN-vector, replaced by solution vector
[in]nsize of V, B
[in]nbdrborder size
[in]nbndband width
[in]inv=1 calculate band part of inverse (for pulls), >1 calculate complete inverse
[out]nrankrank of matrix V
[out]vbndband part of V
[out]vbdrborder part of V
[out]auxsolutions for border rows
[out]vbkmatrix for border solution
[out]vzruborder solution
[out]scdiagworkspace (D)
[out]scflagworkspace (I)
[out]evdminmin. eigenvalue of diagonal matrix from band decomposition
[out]evdmaxmax. eigenvalue of diagonal matrix from band decomposition

Definition at line 3118 of file mpnum.f90.

References dbcdec(), dbcinb(), dbcinv(), dbcslv(), and sqminv().

Referenced by loopbf().

◆ sqmibb2()

subroutine sqmibb2 ( real(mpd), dimension((n*n+n)/2), intent(inout)  v,
real(mpd), dimension(n), intent(out)  b,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  nbdr,
integer(mpi), intent(in)  nbnd,
integer(mpi), intent(in)  inv,
integer(mpi), intent(out)  nrank,
real(mpd), dimension(n*(nbnd+1)), intent(out)  vbnd,
real(mpd), dimension(n*nbdr), intent(out)  vbdr,
real(mpd), dimension(n*nbdr), intent(out)  aux,
real(mpd), dimension((nbdr*nbdr+nbdr)/2), intent(out)  vbk,
real(mpd), dimension(nbdr), intent(out)  vzru,
real(mpd), dimension(nbdr), intent(out)  scdiag,
integer(mpi), dimension(nbdr), intent(out)  scflag,
real(mpd), intent(out)  evdmin,
real(mpd), intent(out)  evdmax 
)

Band bordered matrix.

Obtain solution of a system of linear equations with symmetric band bordered matrix (V * X = B), on request inverse is calculated. For band part root-free Cholesky decomposition and forward/backward substitution is used.

Use decomposition in band and border part for block matrix algebra:

| A  Ct |   | x1 |   | b1 |        , A  is the band part
|       | * |    | = |    |        , Ct is the mixed part
| C  D  |   | x2 |   | b2 |        , D  is the border part
Parameters
[in,out]vsymmetric N-by-N matrix in symmetric storage mode (V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, ...), replaced by inverse matrix
[in,out]bN-vector, replaced by solution vector
[in]nsize of V, B
[in]nbdrborder size
[in]nbndband width
[in]inv=1 calculate band part of inverse (for pulls), >1 calculate complete inverse
[out]nrankrank of matrix V
[out]vbndband part of V
[out]vbdrborder part of V
[out]auxsolutions for border rows
[out]vbkmatrix for border solution
[out]vzruborder solution
[out]scdiagworkspace (D)
[out]scflagworkspace (I)
[out]evdminmin. eigenvalue of diagonal matrix from band decomposition
[out]evdmaxmax. eigenvalue of diagonal matrix from band decomposition

Definition at line 3390 of file mpnum.f90.

References dbcdec(), dbcinb(), dbcinv(), dbcslv(), and sqminv().

Referenced by loopbf().

◆ sqminl()

subroutine sqminl ( real(mpd), dimension((int(n,mpl)*int(n,mpl)+int(n,mpl))/2), intent(inout)  v,
real(mpd), dimension(n), intent(out)  b,
integer(mpi), intent(in)  n,
integer(mpi), intent(out)  nrank,
real(mpd), dimension(n), intent(out)  diag,
integer(mpi), dimension(n), intent(out)  next,
real(mpd), dimension(n), intent(out)  vk,
integer(mpi), intent(in)  mon 
)

Matrix inversion for LARGE matrices.

Like SQMINV, additional parallelization with OpenMP.

Parameters
[in,out]Vsymmetric N-by-N matrix in symmetric storage mode (V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, ...), replaced by inverse matrix
[in,out]BN-vector, replaced by solution vector
[in]Nsize of V, B
[out]NRANKrank of matrix V
[out]DIAGdouble precision scratch array
[out]NEXTinteger aux array
[out]VKdouble precision scratch array (pivot)
[in]MONflag for progress monitoring

Definition at line 230 of file mpnum.f90.

References monpgs().

Referenced by minver().

◆ sqminv()

subroutine sqminv ( real(mpd), dimension((n*n+n)/2), intent(inout)  v,
real(mpd), dimension(n), intent(out)  b,
integer(mpi), intent(in)  n,
integer(mpi), intent(out)  nrank,
real(mpd), dimension(n), intent(out)  diag,
integer(mpi), dimension(n), intent(out)  next 
)

Matrix inversion and solution.

Obtain solution of a system of linear equations with symmetric matrix (V * X = B) and the inverse.

Method of solution is by elimination selecting the pivot on the diagonal each stage. The rank of the matrix is returned in NRANK. For NRANK ne N, all remaining rows and cols of the resulting matrix V and the corresponding elements of B are set to zero.

Parameters
[in,out]Vsymmetric N-by-N matrix in symmetric storage mode (V(1) = V11, V(2) = V12, V(3) = V22, V(4) = V13, ...), replaced by inverse matrix
[in,out]BN-vector, replaced by solution vector
[in]Nsize of V, B
[out]NRANKrank of matrix V
[out]DIAGdouble precision scratch array
[out]NEXTINTEGER(mpi) aux array

Definition at line 97 of file mpnum.f90.

Referenced by ckpgrp(), feasma(), loopbf(), sqmibb(), and sqmibb2().

◆ vabdec()

subroutine vabdec ( integer(mpi), intent(in)  n,
real(mpd), dimension(*), intent(inout)  val,
integer(mpi), dimension(n), intent(in)  ilptr 
)

Variable band matrix decomposition.

Decomposition: A = L D L^T

Variable-band matrix row Doolittle decomposition. A variable-band NxN symmetric matrix, also called skyline, is stored row by row in the array VAL(.). For each row every coefficient between the first non-zero element in the row and the diagonal is stored. The pointer array ILPTR(N) contains the indices in VAL(.) of the diagonal elements. ILPTR(1) is always 1, and ILPTR(N) is equal to the total number of coefficients stored, called the profile. The form of a variable-band matrix is preserved in the L D L^T decomposition no fill-in is created ahead in any row or ahead of the first entry in any column, but existing zero-values will become non-zero. The decomposition is done "in-place".

Parameters
[in]nsize of matrix
[in,out]valvariable-band matrix, replaced by D,L
[in]ilptrpointer array

Definition at line 1012 of file mpnum.f90.

◆ vabmmm()

subroutine vabmmm ( integer(mpi), intent(in)  n,
real(mpd), dimension(*), intent(inout)  val,
integer(mpi), dimension(n), intent(in)  ilptr 
)

Variable band matrix print minimum and maximum.

Parameters
[in]nsize of matrix
[in,out]valvariable-band matrix, replaced by D,L
[in]ilptrpointer array

Definition at line 1124 of file mpnum.f90.

◆ vabslv()

subroutine vabslv ( integer(mpi), intent(in)  n,
real(mpd), dimension(*), intent(inout)  val,
integer(mpi), dimension(n), intent(in)  ilptr,
real(mpd), dimension(n), intent(inout)  x 
)

Variable band matrix solution.

The matrix equation A X = B is solved. The matrix is assumed to decomposed before using VABDEC. The array X(N) contains on entry the right-hand-side B(N); at return it contains the solution.

Parameters
[in]nsize of matrix
[in,out]valdecomposed variable-band matrix
[in]ilptrpointer array
[in,out]xr.h.s vector B, replaced by solution vector X

Definition at line 1161 of file mpnum.f90.