General linear algebra routines. More...

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)
	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)

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.

***** 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]	g	symmetric matrix, replaced by D,L
[in]	n	size of matrix
[out]	NRANK	rank of matrix g
[out]	EVMAX	largest element in D
[out]	EVMIN	smallest element in D
[in]	MON	flag 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]	n	number of sigmas
[in]	nd	ndf

Returns: Chi2/ndf cut value

Definition at line 2078 of file mpnum.f90.

References chindl(), and mpdef::mps.

Referenced by chindl().

◆ 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]	g	symmetric matrix, replaced by D,L
[in]	n	size 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]	g	decomposed symmetric matrix
[in,out]	v	inverse matrix
[in]	n	size 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]	g	decomposed symmetric matrix
[in,out]	x	r.h.s vector B, replaced by solution vector X
[in]	n	size 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]	g	decomposed symmetric matrix
[in,out]	x	rhs/solution
[in]	n	size 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]	V	symmetric N-by-N matrix
[in]	A	general M-by-N matrix
[in,out]	W	symmetric M-by-M matrix
[in]	N	columns of A
[in]	M	rows 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(2nm+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]	V	symmetric N-by-N matrix
[in]	A	sparse M-by-N matrix, content
[in]	IS	sparse M-by-N matrix, structure
[in,out]	W	symmetric M-by-M matrix
[in]	N	columns of A
[in]	M	rows of A
[in]	iopt	(<>0: don't reset W)
[in]	SC	scratch 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]	n	vector size
[in]	dx	vector
[in,out]	dy	vector
[in]	da	scalar

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]	n	vector size
[in]	dx	vector
[in]	dy	vector

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]	A	general M-by-N matrix (A11 A12 ... A1N A21 A22 ...)
[in]	X	N vector
[out]	Y	= M vector
[in]	M	rows of A
[in]	N	columns 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]	lun	unit number
[in]	x	vector
[in]	v	symmetric matrix
[in]	n	size 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]	lun	unit number
[in]	v	symmetric matrix
[in]	n	size 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]	v	symmetric matrix
[in]	a	vector
[out]	b	vector B = V * A
[in]	n	size 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]	v	symmetric matrix
[in]	a	vector
[out]	b	product vector B = V * A
[in]	n	size 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]	N	size of matrix
[in]	DIAG	diagonal elements
[in]	U	transformation matrix
[out]	V	symmmetric 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]	n	size of matrix
[out]	diag	diagonal elements
[out]	u	transformation matrix
[in]	v	symmetric matrix, unchanged
[out]	work	work array
[out]	iwork	work 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]	N	size of matrix
[in]	DIAG	diagonal elements
[in]	U	transformation matrix
[in]	B	r.h.s. of matrix equation (unchanged)
[out]	X	solution vector
[out]	WORK	work 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)+nm+(mm+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]	n	size of symmetric matrix
[in]	m	number of constrains
[in]	ls	flag for skyline decomposition
[in,out]	c	combined variable-band + constraints matrix, replaced by decomposition
[in,out]	india	pointer array
[out]	nrkd	removed components
[out]	nrkd2	removed 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]	n	size of symmetric matrix
[in]	m	number of constrains
[in]	b	number of constraint blocks
[in]	nm	(integrated) size of constraint blocks (<=N*M)
[in]	ls	flag for skyline decomposition
[in,out]	c	combined variable-band + constraints matrix, replaced by decomposition
[in,out]	india	pointer array
[in]	l	constraint (block matrix row) offsets
[out]	nrkd	removed components
[out]	nrkd2	removed 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)+nm+(mm+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]	n	size of symmetric matrix
[in]	m	number of constrains
[in]	c	decomposed combined variable-band + constraints matrix
[in]	india	pointer array
[in,out]	x	r.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]	n	size of symmetric matrix
[in]	m	number of constrains
[in]	b	number of constraint blocks
[in]	nm	(integrated) size of constraint blocks (<=N*M)
[in]	c	decomposed combined variable-band + constraints matrix
[in]	india	pointer array
[in]	l	constraint (block matrix row) offsets
[in,out]	x	r.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]	a	array of keys
[in]	n	number 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]	n	size of matrix
[in,out]	c	decomposed variable-band matrix
[in]	india	pointer array
[in,out]	x	r.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]	n	size of matrix
[in,out]	c	variable-band matrix, replaced by L
[in]	india	pointer array
[out]	nrkd	removed 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]	n	size of matrix
[in,out]	c	decomposed variable-band matrix
[in]	india	pointer array
[in,out]	x	r.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]	n	size of matrix
[in,out]	c	decomposed variable-band matrix
[in]	india	pointer array
[in,out]	x	r.h.s vector B, replaced by solution vector X
[in]	i0	offset of non zero region in x
[in]	ns	size 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]	p	number of constraints
[in]	n	size of diagonal matrix
[in]	c	diagonal matrix (changed if c=cu as actual parameters)
[out]	cu	1/sqrt(c)
[in,out]	a	constraint matrix (size n*p), modified
[out]	s	Cholesky decomposed symmetric (P,P) matrix
[out]	nrkd	removed 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]	p	number of constraints
[in]	n	size of diagonal matrix
[in]	b	number of constraint blocks
[in]	nm	(integrated) size of constraint blocks (<=N*M)
[in]	c	diagonal matrix (changed if c=cu as actual parameters)
[out]	cu	1/sqrt(c)
[in,out]	a	modified constraint (block) matrix, modified
[in]	l	constraint (block matrix row) offsets
[out]	s	Cholesky decomposed symmetric (P,P) matrix
[out]	nrkd	removed 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]	p	number of constraints
[in]	n	size of diagonal matrix
[in]	cu	1/sqrt(c)
[in]	a	modified constraint matrix (size n*p)
[in]	s	Cholesky decomposed symmetric (P,P) matrix
[out]	x	result vector
[in]	y	rhs 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]	p	number of constraints
[in]	n	size of diagonal matrix
[in]	b	number of constraint blocks
[in]	nm	(integrated) size of constraint blocks (<=N*M)
[in]	cu	1/sqrt(c)
[in]	a	modified constraint (block) matrix
[in]	l	constraint (block matrix row) offsets
[in]	s	Cholesky decomposed symmetric (P,P) matrix
[out]	x	result vector
[in]	y	rhs 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]	a	vector of integers, sorted at return
[in]	n	size 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]	a	vector (quadruplet) of integers, sorted at return and an index
[in,out]	b	vector of long integers, sorted at return and an index
[in]	n	size 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]	a	vector (pair) of integers, sorted at return and an index
[in]	n	size 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]	a	vector (pair) of integers, sorted at return
[in]	n	size 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
	)

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]	v	symmetric 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]	b	N-vector, replaced by solution vector
[in]	n	size of V, B
[in]	nbdr	border size
[in]	nbnd	band width
[in]	inv	=1 calculate band part of inverse (for pulls), >1 calculate complete inverse
[out]	nrank	rank of matrix V
[out]	vbnd	band part of V
[out]	vbdr	border part of V
[out]	aux	solutions for border rows
[out]	vbk	matrix for border solution
[out]	vzru	border solution
[out]	scdiag	workspace (D)
[out]	scflag	workspace (I)

Definition at line 3116 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
	)

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]	v	symmetric 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]	b	N-vector, replaced by solution vector
[in]	n	size of V, B
[in]	nbdr	border size
[in]	nbnd	band width
[in]	inv	=1 calculate band part of inverse (for pulls), >1 calculate complete inverse
[out]	nrank	rank of matrix V
[out]	vbnd	band part of V
[out]	vbdr	border part of V
[out]	aux	solutions for border rows
[out]	vbk	matrix for border solution
[out]	vzru	border solution
[out]	scdiag	workspace (D)
[out]	scflag	workspace (I)

Definition at line 3372 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]	V	symmetric 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]	B	N-vector, replaced by solution vector
[in]	N	size of V, B
[out]	NRANK	rank of matrix V
[out]	DIAG	double precision scratch array
[out]	NEXT	integer aux array
[out]	VK	double precision scratch array (pivot)
[in]	MON	flag 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]	V	symmetric 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]	B	N-vector, replaced by solution vector
[in]	N	size of V, B
[out]	NRANK	rank of matrix V
[out]	DIAG	double precision scratch array
[out]	NEXT	INTEGER(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]	n	size of matrix
[in,out]	val	variable-band matrix, replaced by D,L
[in]	ilptr	pointer 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]	n	size of matrix
[in,out]	val	variable-band matrix, replaced by D,L
[in]	ilptr	pointer 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]	n	size of matrix
[in,out]	val	decomposed variable-band matrix
[in]	ilptr	pointer array
[in,out]	x	r.h.s vector B, replaced by solution vector X

Definition at line 1161 of file mpnum.f90.

Functions/Subroutines

Detailed Description

Function/Subroutine Documentation

◆ chdec2()

◆ chindl()

◆ choldc()

◆ cholin()

◆ cholsl()

◆ chslv2()

◆ dbavat()

◆ dbavats()

◆ dbaxpy()

◆ dbdot()

◆ dbgax()

◆ dbmprv()

◆ dbprv()

◆ dbsvx()

◆ dbsvxl()

◆ devinv()

◆ devrot()

◆ devsig()

◆ devsol()

◆ equdec()

◆ equdecs()

◆ equslv()

◆ equslvs()

◆ heapf()

◆ lltbwd()

◆ lltdec()

◆ lltfwd()

◆ lltfwds()

◆ precon()

◆ precons()

◆ presol()

◆ presols()

◆ sort1k()

◆ sort22l()

◆ sort2i()

◆ sort2k()

◆ sqmibb()

◆ sqmibb2()

◆ sqminl()

◆ sqminv()

◆ vabdec()

◆ vabmmm()

◆ vabslv()