Solver for block tri-diagonal matrices. [Kalyan S. Perumalla, ORNL, 2009-2011]. More...

Data Types
type	blacsparameters
	BLACS/PBLAS information. More...

type	blacsprocessgrid
	BLACS/PBLAS process grid information. More...

type	levelelement
	Data associated with each row at each level. More...

type	mastertoslavemapping
	Master-to-slave mapping. More...

type	pblaslevelparameters
	Level-specific PBLAS information. More...

type	pblasstats

type	pblastemparray

type	solutionelement
	Solution of selected rows of interest to this rank. More...

type	timecount
	Statistics (timing, etc.) More...

Functions/Subroutines
subroutine	bclockinit ()
	A convenience function to be able to change clock routine(s) easily.

subroutine	bsystemclock (ts)
	A convenience function to be able to change clock routine easily.

subroutine	fl (u)
	A convenience function to be able to change or turn off easily.

subroutine	chargememory (bytes)
	Convenience routine to track allocated memory sizes.

subroutine	chargetime (tot, t2, t1, cnt)
	Convenience routine to accumulate timing values.

subroutine	timecountinit (tc)
	Initialize stats.

subroutine	timecountprint (tc, msg)
	Print stats.

subroutine	plbinitstats ()
	Init statistics.

subroutine	plbprintstats ()
	Print statistics so far.

subroutine	plbinitialize
	Initialize before foward or backward starts.

subroutine	plbfinalize
	Finalize after foward or backward are done.

subroutine	plbforwardinitialize
	Initialize before any forward level starts.

subroutine	plbforwardfinalize
	Finalize after all forward levels are done.

subroutine	plbbackwardinitialize
	Initialize before any backward level starts.

subroutine	plbbackwardfinalize
	Finalize after all backward levels are done.

subroutine	plbforwardinitializelevel (lvl, ammaster)
	Initialize before the forward solve at a given level starts.

subroutine	plbforwardfinalizelevel (lvl, ammaster)
	Finalize after the forward solve at a given level ends.

subroutine	plbbackwardinitializelevel (lvl, ammaster)
	Initialize before the backward solve at a given level starts.

subroutine	plbbackwardfinalizelevel (lvl, ammaster)
	Finalize after the backward solve at a given level ends.

subroutine	determinemasterslaveranks ()
	Determine the ranks of master and slaves at this level.

subroutine	masterbcastnextop (nextop)
	Inform the next operation to be done (calling rank is in master mode)

subroutine	slavegetnextop (nextop)
	Determine the next operation to be done (calling rank is in slave mode)

subroutine	slaveservice ()
	Slave's operation loop.

subroutine	extractsubmatrix (bszr, bszc, pnr, pnc, pi, pj, a, nrows, ncols, suba, subnrows, subncols)
	Extracts from A the submatrix corresponding to process-grid element (pi,pj)

subroutine	mastersendmatrix (a, nrows, ncols, ssuba, ssubnrows, ssubncols)
	Send the matrix by master (caller) to all slaves The submatrix belonging to the caller (master is also one of the slaves) is returned in the ssub* arguments (i.e., sent to self)

subroutine	slavereceivematrix (nrows, ncols, suba, subnrows, subncols)
	Receive the fragment of matrix sent by master to slave (caller)

subroutine	masterbcastvalue (val)
	Send the given value by master (caller) to all slaves.

subroutine	slavereceivevalue (val)
	Receive a value sent by master to slave (caller)

subroutine	injectsubmatrix (bszr, bszc, pnr, pnc, pi, pj, a, nrows, ncols, suba, subnrows, subncols)
	Incorporates into A the submatrix corresp. to process-grid element (pi,pj)

subroutine	injectsubvector (bszr, pnr, pi, v, nrows, subv, subnrows)
	Incorporates into V the subvector corresp. to process-grid element (pi,0)

subroutine	masterrecvmatrix (a, nrows, ncols, ssuba, ssubnrows, ssubncols)
	Receive the submatrices by master (caller) from all slaves The submatrix belonging to the caller (master is also one of the slaves) is given in the ssub* arguments (i.e., received from self)

subroutine	slavesendmatrix (nrows, ncols, suba, subnrows, subncols)
	Send the fragment of matrix by slave (caller) to master.

subroutine	masterrecvvector (v, nrows, ssubv, ssubnrows)
	Receive the subvector by master (caller) from all slaves in column 1 of pgrid The subvector belonging to the caller (master is also one of the slaves) is given in the ssub* arguments (i.e., received from self)

subroutine	slavesendvector (nrows, subv, subnrows)
	Send the fragment of vector by slave (caller) to master Only slaves who belong to zero'th column of pgrid send to master Other slave simply do no-op for this routine.

subroutine	plbdgemm (alpha, a, b, beta, c)
	Encapsulates BLAS' DGEMM functionality.

subroutine	slavedgemm ()
	DGEMM support from slave.

subroutine	plbdgemv (alpha, a, x, beta, y)
	Encapsulates BLAS' DGEMV functionality.

subroutine	plbdgetrf (a, piv, info)
	Encapsulates BLAS/LAPACK's DGETRF functionality.

subroutine	slavedgetrf ()
	DGETRF support from slave.

subroutine	plbdgetrs (nrhs, a, piv, b, info)
	Encapsulates BLAS/LAPACK's DGETRS functionality.

subroutine	slavedgetrs ()
	DGETRS support from slave.

subroutine	initialize (do_mpiinit, inn, inm)
	To be invoked, before solve, by user.

subroutine	setmatrixrow (globrow, l, d, u)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	setmatrixrowl (globrow, l)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	setmatrixrowd (globrow, d)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	setmatrixrowu (globrow, u)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	setmatrixrowcoll (globrow, lj, j)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	setmatrixrowcold (globrow, dj, j)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	setmatrixrowcolu (globrow, uj, j)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	setmatrixrhs (globrow, b)
	To be invoked, before solve, after initialize, by user, to set up the matrix.

subroutine	getmatrixrowcoll (globrow, lj, j)
	Can be invoked after initialize, by user, to get a copy of the matrix.

subroutine	getmatrixrowcold (globrow, dj, j)
	Can be invoked after initialize, by user, to get a copy of the matrix.

subroutine	getmatrixrowcolu (globrow, uj, j)
	Can be invoked after initialize, by user, to get a copy of the matrix.

subroutine	getmatrixrhs (globrow, b)
	Can be invoked after initialize, by user, to get a copy of the RHS.

subroutine	getsolutionvector (globrow, x)
	To be invoked, after solve, by user to get the solution vector of a local row.

subroutine	setidentitytestcase
	To be invoked, before solve, after initialize, by user, to generate a test.

subroutine	setidentityrhs
	To be invoked, before solve, after SetIdentityCase, by user, to generate a test.

subroutine	setrandomtestcase
	To be invoked, before solve, after initialize, by user, to generate a test.

subroutine	setrandomrhs (randseedoff)
	To be invoked, before solve, after SetRandomTestCase, by user, to generate a an RHS.

subroutine	finalize (do_mpifinalize)
	To be invoked, after solve, by user.

logical function	iseven (num)
	Is given integer even?

logical function	isodd (num)
	Is given integer odd?

integer function	lr2gr (locrow, level)
	Determines the global row number of a given local row number at a given level.

integer function	gr2lr (globrow, level)
	Determines the local row number of "globrow" global row number when globrow participates at a given level; returns zero if this globrow does not operate at the given level.

integer function	gr2rank (globrow)
	Determines the rank of the task holding the given global row (at level 1)

integer function	lr2rank (locrow, level)
	Determines the rank of the task holding the given local row (at a given level)

subroutine	computeforwardoddrowhats (locrow, level, startlocrow, endlocrow, bonly)
	Compute the matrix multiplications in the forward solve at some/any level The row should be odd at the level, The locrow should obey startlocrow <= locrow <= endlocrow.

subroutine	forwardsolve
	BCYCLIC forward phase; to be called after Initialize.

subroutine	forwardupdateb
	BCYCLIC forward phase to deal with a new b; to be called after a SetRHS that may have been invoked after a ForwardSolve, before BackwardSolve.

subroutine	backwardsolve
	BCYCLIC backward phase; to be called after ForwardSolve.

subroutine	verifysolution
	Verify the RMS error of solution after backward solve.

subroutine	checksymmetry (asymindx)

subroutine	storediagonal (blkrow, colnum, buf)

subroutine	writeblocks (flag)

subroutine	displayblocks (flag)

subroutine	parmatvec (invec, outvec, outveclength)

subroutine	setblockrowcol (globrow, colnum, buf, opt)

subroutine	getcolsum (cs)

subroutine	initscalefactors

subroutine	getscalefactors (js, scalevector)

subroutine	finalizescalefactors

subroutine	parallelscaling (lmp, colsum)

subroutine	findminmax_tri (lmp)

subroutine	resetdiagonal (lmp, breset)

Variables
integer, parameter	rprec = SELECTED_REAL_KIND(15, 300)

integer, parameter	iprec = SELECTED_INT_KIND(8)

integer, parameter	cprec = KIND((1.0_rprec, 1.0_rprec))

integer, parameter	dp = rprec

type(levelelement), dimension(:,:), allocatable	lelement
	Data for each row at each level on this rank The first dimension is the level number [1..L], L=#levels of forward solve The 2nd dimension is the row number [1..K+g], K=#rows at level 1 on this rank The +g is for incoming results from neighbors, 0<=g<=2.

type(levelelement), dimension(:), allocatable	orig
	Initial problem specification saved for verification at end of solution The dimension is the row number [1..K], K=#rows at level 1 on this rank.

type(solutionelement), dimension(:), allocatable	selement
	The solution The dimension is the global (level 1) row number [1..N].

integer	p
	Num of "master" tasks, adjusted such that P=Min(nranks,N)

integer	n
	Num of row blocks in input block tri-diag matrix.

integer	m
	Size of each square sub-matrix block.

integer	startglobrow
	Starting row number of this processor at first level.

integer	endglobrow
	Ending row number (inclusive) at first level.

integer	nlevels
	Total number of levels in recursion.

logical	matdirtied
	Has a new matrix been set after ForwardSolve?

logical	rhsdirtied
	Has a new RHS been set after ForwardSolve?

character *100	kenvvar

character *100	kenvval

logical	kpdbg
	Should debugging output be written? : SKS on Nov 9, 2011.

integer	ofu
	Output file unit.

integer	pfu
	Problem file unit.

integer	mpi_err

logical	writeproblemfile
	Should save the random test case to a file?

logical	writesolution
	Should dump solution in output?

logical	usebarriers
	Should barriers be invoked between levels?

real	membytes
	A running count of memory allocated so far.

real	dpsz
	Byte size of a double precision variable.

real	intsz
	Byte size of an integer variable.

real	ptrsz
	Byte size of a pointer variable.

real(dp)	one
	1.0

real(dp)	zero
	0.0

logical(dp)	l_colscale =.FALSE.
	set true if column scaling is applied

logical	use_mpiwtime
	Use MPI's timer function?

double precision	loctimer1

double precision	loctimer2
	Stopwatch snapshots to use locally.

double precision	mattimer1

double precision	mattimer2
	Stopwatch snapshots for matrix ops.

double precision	globtimer1

double precision	globtimer2
	Stopwatch snapshots to use globally.

double precision	timerfreq
	Timer frequency in Hz (count rate)

double precision	tottime
	Total elapsed time.

integer	totcount
	Total count.

double precision	totcommtime
	Total communication time.

integer	totcommcount
	Total communication operation count.

double precision	totinvtime
	Total time to invert matrices.

integer	totinvcount
	Total count of matrix inversions.

double precision	totmatmultime
	Total time for matrix-matrix multiplications.

integer	totmatmulcount
	Total count of matrix-matrix multiplications.

double precision	totmatsoltime
	Total time for matrix-solve multiplications.

integer	totmatsolcount
	Total count of matrix solutions.

real(dp), public	maxeigen_tri
	maximum eigenvalue from Gerschgorin theorem

real(dp), public	dmin_tri
	minimum non-zero diagonal element, used to scale Levenberg parameter (SPH 08-31-16)

real(dp), dimension(:), allocatable	origdiagelement
	Unmodified diagonal element : SKS.

real(dp), dimension(:), allocatable	topscalefac

real(dp), dimension(:), allocatable	botscalefac

logical	doblasonly
	BLACS/PBLAS options.

logical	doblacscomm

type(blacsparameters)	blacs

type(pblaslevelparameters)	pblas

integer, parameter	op_none = 0
	Master-to-slave commands.

integer, parameter	op_done = 1

integer, parameter	op_dgemm = 2

integer, parameter	op_dgemv = 3

integer, parameter	op_dgetrf = 4

integer, parameter	op_dgetrs = 5

type(pblasstats)	pstats

Detailed Description

Solver for block tri-diagonal matrices. [Kalyan S. Perumalla, ORNL, 2009-2011].

Function/Subroutine Documentation

◆ bsystemclock()

subroutine blocktridiagonalsolver::bsystemclock ( double precision, intent(inout) ts )

A convenience function to be able to change clock routine easily.

Parameters

[in,out] ts Time snapshot

◆ chargememory()

subroutine blocktridiagonalsolver::chargememory ( real, intent(in) bytes )

Convenience routine to track allocated memory sizes.

Parameters

[in] bytes Number of bytes allocated

◆ chargetime()

subroutine blocktridiagonalsolver::chargetime	(	double precision, intent(inout)	tot,
		double precision, intent(in)	t2,
		double precision, intent(in)	t1,
		integer, intent(inout)	cnt
	)

Convenience routine to accumulate timing values.

Parameters

[in,out]	tot	Total time
[in]	t2	Timer value at end of segment being timed
[in]	t1	Timer value at start of timed segment
[in,out]	cnt	Counter of number of samples

◆ computeforwardoddrowhats()

subroutine blocktridiagonalsolver::computeforwardoddrowhats	(	integer, intent(in)	locrow,
		integer, intent(in)	level,
		integer, intent(in)	startlocrow,
		integer, intent(in)	endlocrow,
		logical, intent(in)	bonly
	)

Compute the matrix multiplications in the forward solve at some/any level The row should be odd at the level, The locrow should obey startlocrow <= locrow <= endlocrow.

Parameters

[in]	locrow	Row number at a given level
[in]	level	Local row number is at this level
[in]	startlocrow	Top-most local row number at given level
[in]	endlocrow	Bot-most local row number at given level
[in]	bonly	Show compute/update only the b, and not L,D,U?

◆ finalize()

subroutine blocktridiagonalsolver::finalize ( logical, intent(in) do_mpifinalize )

To be invoked, after solve, by user.

Parameters

[in] do_mpifinalize Invoke MPI_Finalize?

◆ fl()

subroutine blocktridiagonalsolver::fl ( integer, intent(in) u )

A convenience function to be able to change or turn off easily.

Parameters

[in] u Unit number to flush

◆ getmatrixrhs()

subroutine blocktridiagonalsolver::getmatrixrhs	(	integer, intent(in)	globrow,
		real(dp), dimension(:), intent(out)	b
	)

Can be invoked after initialize, by user, to get a copy of the RHS.

Parameters

[in]	globrow	Original/global block-row num in [1..N]
[out]	b	RHS column corr. to globrow

◆ getmatrixrowcold()

subroutine blocktridiagonalsolver::getmatrixrowcold	(	integer, intent(in)	globrow,
		real(dp), dimension(:), intent(out)	dj,
		integer, intent(in)	j
	)

Can be invoked after initialize, by user, to get a copy of the matrix.

Parameters

[in]	globrow	Original/global block-row num in [1..N]
[out]	dj	j'th colum of D at globrow
[in]	j	column number of L at globrow that is being set

◆ getmatrixrowcoll()

subroutine blocktridiagonalsolver::getmatrixrowcoll	(	integer, intent(in)	globrow,
		real(dp), dimension(:), intent(out)	lj,
		integer, intent(in)	j
	)

Can be invoked after initialize, by user, to get a copy of the matrix.

Parameters

[in]	globrow	Original/global block-row num in [1..N]
[out]	lj	j'th colum of L at globrow; 1st L is always 0
[in]	j	column number of L at globrow that is being set

◆ getmatrixrowcolu()

subroutine blocktridiagonalsolver::getmatrixrowcolu	(	integer, intent(in)	globrow,
		real(dp), dimension(:), intent(out)	uj,
		integer, intent(in)	j
	)

Can be invoked after initialize, by user, to get a copy of the matrix.

Parameters

[in]	globrow	Original/global block-row num in [1..N]
[out]	uj	j'th colum of L at globrow; Nth U is always 0
[in]	j	column number of U at globrow that is being set

◆ getsolutionvector()

subroutine blocktridiagonalsolver::getsolutionvector	(	integer, intent(in)	globrow,
		real(dp), dimension(:), intent(out)	x
	)

To be invoked, after solve, by user to get the solution vector of a local row.

Parameters

[in]	globrow	Original/global block-row num in [1..N]
[out]	x	Solution column corr. to globrow

◆ gr2lr()

integer function blocktridiagonalsolver::gr2lr	(	integer, intent(in)	globrow,
		integer, intent(in)	level
	)

Determines the local row number of "globrow" global row number when globrow participates at a given level; returns zero if this globrow does not operate at the given level.

Parameters

[in]	globrow	Row number in original input (level 1)
[in]	level	Level at which globrow's position is needed

Returns: Returned: local row number of global row, at given level

◆ gr2rank()

integer function blocktridiagonalsolver::gr2rank ( integer, intent(in) globrow )

Determines the rank of the task holding the given global row (at level 1)

Parameters

[in] globrow Row number in original input (level 1)

Returns: Returned: Rank holding the given global row

◆ initialize()

subroutine blocktridiagonalsolver::initialize	(	logical, intent(in)	do_mpiinit,
		integer, intent(in)	inn,
		integer, intent(in)	inm
	)

To be invoked, before solve, by user.

Parameters

[in]	do_mpiinit	Invoke MPI_Init?
[in]	inn	Num of row blocks in input block tri-diag matrix
[in]	inm	Size of each square sub-matrix block

◆ iseven()

logical function blocktridiagonalsolver::iseven ( integer, intent(in) num )

Is given integer even?

Parameters

[in] num A number

◆ isodd()

logical function blocktridiagonalsolver::isodd ( integer, intent(in) num )

Is given integer odd?

Parameters

[in] num A number

◆ lr2gr()

integer function blocktridiagonalsolver::lr2gr	(	integer, intent(in)	locrow,
		integer, intent(in)	level
	)

Determines the global row number of a given local row number at a given level.

Parameters

[in]	locrow	local row number of global row, at given level
[in]	level	Level at which locrow's position is given

Returns: Returned: Row number in original input (level 1)

◆ lr2rank()

integer function blocktridiagonalsolver::lr2rank	(	integer, intent(in)	locrow,
		integer, intent(in)	level
	)

Determines the rank of the task holding the given local row (at a given level)

Parameters

[in]	locrow	Row number at a given level
[in]	level	Local row number is at this level

Returns: Returned: Rank holding the given global row

◆ plbdgemm()

subroutine blocktridiagonalsolver::plbdgemm	(	real(dp), intent(in)	alpha,
		real(dp), dimension(:,:), intent(in)	a,
		real(dp), dimension(:,:), intent(in)	b,
		real(dp), intent(in)	beta,
		real(dp), dimension(:,:), intent(inout)	c
	)

Encapsulates BLAS' DGEMM functionality.

Parameters

[in]	b	Each matrix of size MxM
[in,out]	c	Matrix of size MxM

◆ plbdgemv()

subroutine blocktridiagonalsolver::plbdgemv	(	real(dp), intent(in)	alpha,
		real(dp), dimension(:,:), intent(in)	a,
		real(dp), dimension(:), intent(in)	x,
		real(dp), intent(in)	beta,
		real(dp), dimension(:), intent(inout)	y
	)

Encapsulates BLAS' DGEMV functionality.

Parameters

[in]	a	Matrix of size MxM
[in]	x	Vector of size Mx1
[in,out]	y	Vector of size Mx1

◆ plbdgetrf()

subroutine blocktridiagonalsolver::plbdgetrf	(	real(dp), dimension(:,:), intent(inout)	a,
		integer, dimension(:), intent(out)	piv,
		integer, intent(inout)	info
	)

Encapsulates BLAS/LAPACK's DGETRF functionality.

Parameters

[in,out]	a	Matrix of size MxM
[out]	piv	Mx1

◆ plbdgetrs()

subroutine blocktridiagonalsolver::plbdgetrs	(	integer, intent(in)	nrhs,
		real(dp), dimension(:,:), intent(in)	a,
		integer, dimension(:), intent(in)	piv,
		real(dp), dimension(:,:), intent(inout)	b,
		integer	info
	)

Encapsulates BLAS/LAPACK's DGETRS functionality.

Parameters

[in]	a	Matrix of size MxM
[in]	piv	Mx1
[in,out]	b	Matrix of size Mxnrhs