Millepede-II  V04-08-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, ms)
 A V AT product (similarity). More...
 
subroutine dbavats (v, a, is, w, n, ms, 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 sort22 (a, 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 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 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 sqmibb (v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
 Bordered band matrix. More...
 
subroutine sqmibb2 (v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
 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)
   SORT22   sort 2-dim key-array with two additional values (CHK)

Definition in file mpnum.f90.

Function/Subroutine Documentation

◆ chdec2()

subroutine chdec2 ( real(mpd), dimension(*), 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 890 of file mpnum.f90.

References monpgs(), and mpdef::mpl.

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 2071 of file mpnum.f90.

References mpdef::mps.

◆ choldc()

subroutine choldc ( real(mpd), dimension(*), 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 744 of file mpnum.f90.

◆ cholin()

subroutine cholin ( real(mpd), dimension(*), intent(in)  g,
real(mpd), dimension(*), 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 836 of file mpnum.f90.

◆ cholsl()

subroutine cholsl ( real(mpd), dimension(*), 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 790 of file mpnum.f90.

◆ chslv2()

subroutine chslv2 ( real(mpd), dimension(*), 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 952 of file mpnum.f90.

References mpdef::mpl.

Referenced by mchdec().

◆ dbavat()

subroutine dbavat ( real(mpd), dimension(*), intent(in)  v,
real(mpd), dimension(*), intent(in)  a,
real(mpd), dimension(*), intent(inout)  w,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  ms 
)

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]MSrows of A (-rows: don't reset W)
[in]Ncolumns of A

Definition at line 1387 of file mpnum.f90.

Referenced by loopbf().

◆ dbavats()

subroutine dbavats ( real(mpd), dimension(*), intent(in)  v,
real(mpd), dimension(*), intent(in)  a,
integer(mpi), dimension(*), intent(in)  is,
real(mpd), dimension(*), intent(inout)  w,
integer(mpi), intent(in)  n,
integer(mpi), intent(in)  ms,
integer(mpi), dimension(*), 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]MSrows of A (-rows: don't reset W)
[in]Ncolumns of A
[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 1469 of file mpnum.f90.

Referenced by loopbf().

◆ dbaxpy()

subroutine dbaxpy ( integer(mpi), intent(in)  n,
real(mpd), intent(in)  da,
real(mpd), dimension(*), intent(in)  dx,
real(mpd), dimension(*), 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 1232 of file mpnum.f90.

◆ dbdot()

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

Dot product.

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

Definition at line 1201 of file mpnum.f90.

◆ dbgax()

subroutine dbgax ( real(mpd), dimension(*), intent(in)  a,
real(mpd), dimension(*), intent(in)  x,
real(mpd), dimension(*), 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 1350 of file mpnum.f90.

◆ dbmprv()

subroutine dbmprv ( integer(mpi), intent(in)  lun,
real(mpd), dimension(*), intent(in)  x,
real(mpd), dimension(*), 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 1542 of file mpnum.f90.

References mpdef::mpi.

◆ dbprv()

subroutine dbprv ( integer(mpi), intent(in)  lun,
real(mpd), dimension(*), 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 1617 of file mpnum.f90.

◆ dbsvx()

subroutine dbsvx ( real(mpd), dimension(*), intent(in)  v,
real(mpd), dimension(*), intent(in)  a,
real(mpd), dimension(*), 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 1263 of file mpnum.f90.

Referenced by feasib().

◆ dbsvxl()

subroutine dbsvxl ( real(mpd), dimension(*), intent(in)  v,
real(mpd), dimension(*), intent(in)  a,
real(mpd), dimension(*), 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 1307 of file mpnum.f90.

References mpdef::mpl.

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(*), 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]Vsmmmetric matrix

Definition at line 695 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(*), 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 368 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 610 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 648 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(*), 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 2305 of file mpnum.f90.

References lltdec(), and lltfwd().

Referenced by mminrs(), and mminrsqlp().

◆ equslv()

subroutine equslv ( integer(mpi), intent(in)  n,
integer(mpi), intent(in)  m,
real(mpd), dimension(*), intent(in)  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 2371 of file mpnum.f90.

References lltbwd(), and lltfwd().

Referenced by mvsolv().

◆ heapf()

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

Heap sort direct (real).

Real keys A(*), sorted at return.

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

Definition at line 1661 of file mpnum.f90.

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

◆ lltbwd()

subroutine lltbwd ( integer(mpi), intent(in)  n,
real(mpd), dimension(*), intent(in)  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 2265 of file mpnum.f90.

Referenced by equslv().

◆ lltdec()

subroutine lltdec ( integer(mpi), intent(in)  n,
real(mpd), dimension(*), 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 2148 of file mpnum.f90.

Referenced by equdec().

◆ lltfwd()

subroutine lltfwd ( integer(mpi), intent(in)  n,
real(mpd), dimension(*), intent(in)  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 2230 of file mpnum.f90.

Referenced by equdec(), and equslv().

◆ 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 2442 of file mpnum.f90.

Referenced by minresmodule::minres(), 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 2516 of file mpnum.f90.

Referenced by mcsolv().

◆ sort1k()

subroutine sort1k ( integer(mpi), dimension(*), 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 1716 of file mpnum.f90.

References mpdef::parameter, and peend().

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

◆ sort22()

subroutine sort22 ( integer(mpi), dimension(4,*), intent(inout)  a,
integer(mpi), intent(in)  n 
)

Quick sort 2 with index.

Quick sort of A(4,N) integer.

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

Definition at line 1982 of file mpnum.f90.

References mpdef::parameter, and 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 1894 of file mpnum.f90.

References mpdef::parameter, and peend().

Referenced by prpcon().

◆ 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 1801 of file mpnum.f90.

References mpdef::parameter, and peend().

Referenced by peread().

◆ sqmibb()

subroutine sqmibb ( real(mpd), dimension(*), 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 
)

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)

Definition at line 2625 of file mpnum.f90.

References dbcdec(), and sqminv().

Referenced by loopbf().

◆ sqmibb2()

subroutine sqmibb2 ( real(mpd), dimension(*), 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 
)

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)

Definition at line 2881 of file mpnum.f90.

References dbcdec(), and sqminv().

Referenced by loopbf().

◆ sqminl()

subroutine sqminl ( real(mpd), dimension(*), 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 229 of file mpnum.f90.

References monpgs(), and mpdef::mpl.

Referenced by minver().

◆ sqminv()

subroutine sqminv ( real(mpd), dimension(*), 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 96 of file mpnum.f90.

Referenced by 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 1011 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 1123 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 1160 of file mpnum.f90.