Module_Lin_Procedures¶

This file is part of statpack.

statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

statpack 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 Lesser General Public License for more details.

You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)

MODULE EXPORTING SUBROUTINES AND FUNCTIONS FOR THE SOLUTION OF SYSTEMS OF LINEAR EQUATIONS, COMPUTING A TRIANGULAR FACTORIZATION (LU, CHOLESKY), COMPUTING THE INVERSE OF A MATRIX AND COMPUTING THE DETERMINANT OF A MATRIX.

LATEST REVISION : 25/02/2024

`subroutine lu_cmp ( mat, ip, d1, d2, tol, small )`¶

Purpose¶

LU_CMP computes the LU decomposition with partial pivoting and implicit row scaling of a given n-by-n real matrix MAT

P * MAT = L * U

where P is a permutation matrix, L is a n-by-n unit lower triangular matrix and U is a n-by-n upper triangular matrix. P is a permutation matrix, stored in argument IP, such that

P = P(n) * … * P(1)

with P(i) is the identity with row i and IP(i) interchanged.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the coefficient matrix MAT.

On exit, MAT is replaced by the LU decomposition of a rowwise permutation of MAT. The unit diagonal of L is not stored. For solving efficiency, the diagonal reciprocals of the matrix U are saved in the diagonal entries of MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

IP (OUTPUT) integer(i4b), dimension(:)

On exit, IP records the permutations effected by the partial pivoting.

The size of IP must verify: size( IP ) = n .

D1 (OUTPUT) real(stnd)

On exit, if D2 is absent:

D1 = +1, if an even number of interchanges was carried out

D1 = -1, if an odd number of interchanges was carried out

D1 = 0, if MAT is algorithmically singular.

On exit, if D2 is present, D1 is the first component of the determinant of MAT (mantissa of determinant).

D2 (OUTPUT, OPTIONAL) integer(i4b)

If D2 is present, the two components of the determinant of MAT are computed.

On exit, D2 is the second component of the determinant of MAT (exponent of determinant):

determinant(MAT) = scale( D1, D2 )

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4).

If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

SMALL (INPUT, OPTIONAL) real(stnd)

On entry, if the system is singular, replaces a diagonal term of the matrix U if it is smaller in magnitude than the value SMALL using the same sign as the diagonal term. An approximate solution based on this replacement can be obtained if no overflow results.

If SMALL is supplied as less than SAFMIN, the smallest number that can be reciprocated safely, then the value SAFMIN is used in place of SMALL.

Default: SAFMIN, the smallest number that can be reciprocated safely.

Further Details¶

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and D1 is set to zero.

If D1/=0 then the linear system MAT * X = B can be solved with subroutines LU_SOLVE or LU_SOLVE2.

If MAT is algorithmically singular (D1=0), the diagonal terms of U smaller in magnitude than the value SMALL have been replaced by SMALL, using the same sign as the diagonal terms and the decomposition has been completed. An approximate solution based on this replacement can be obtained if no overflow results.

A blocked algorithm is used to compute the factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine lu_cmp2 ( mat, ip, d1, d2, b, matinv, tol, small )`¶

Purpose¶

LU_CMP2 computes the LU decomposition with partial pivoting and implicit row scaling of a given n-by-n real matrix MAT

P * MAT = L * U

where P is a permutation matrix, L is a n-by-n unit lower triangular matrix and U is a n-by-n upper triangular matrix. P is a permutation matrix, stored in argument IP, such that

P = P(n) * … * P(1)

with P(i) is the identity with row i and IP(i) interchanged.

If D2 is present, LU_CMP2 computes the determinant of MAT as

determinant(MAT) = scale( D1, D2 )

If B is present, LU_CMP2 solves the system of linear equations

MAT * X = B

using the LU factorization with scaled partial pivoting of MAT. Here B is a n-vector.

If MATINV is present, LU_CMP2 computes the inverse of MAT

MATINV = MAT**(-1)

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the coefficient matrix MAT.

On exit, MAT is replaced by the LU decomposition of a rowwise permutation of MAT. The unit diagonal of L is not stored. For solving efficiency, the diagonal reciprocals of the matrix U are saved in the diagonal entries of MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

IP (OUTPUT) integer(i4b), dimension(:)

On exit, IP records the permutations effected by the partial pivoting.

The size of IP must verify: size( IP ) = n .

D1 (OUTPUT) real(stnd)

On exit, if D2 is absent:

D1 = +1, if an even number of interchanges was carried out

D1 = -1, if an odd number of interchanges was carried out

D1 = 0, if MAT is algorithmically singular.

On exit, if D2 is present, D1 is the first component of the determinant of MAT (mantissa of determinant).

D2 (OUTPUT, OPTIONAL) integer(i4b)

If D2 is present, the components of the determinant of MAT are computed.

On exit, D2 is the second component of the determinant of MAT (exponent of determinant) :

determinant(MAT) = scale( D1, D2 )

B (INPUT/OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, the solution vector X.

The shape of B must verify: size( B ) = n .

MATINV (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, if MAT is not singular, MATINV contains the inverse of MAT.

The shape of MATINV must verify: size( MATINV, 1 ) = size( MATINV, 2 ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

SMALL (INPUT, OPTIONAL) real(stnd)

On entry, if the system is singular, replaces a diagonal term of the matrix U if it is smaller in magnitude than the value SMALL using the same sign as the diagonal term. An approximate solution for X based on this replacement can be obtained if no overflow results. If SMALL is supplied as less than SAFMIN, the smallest number that can be reciprocated safely, then the value SAFMIN is used in place of SMALL.

Default: SAFMIN, the smallest number that can be reciprocated safely.

Further Details¶

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and D1 is set to zero.

If MAT is algorithmically singular (D1=0), the diagonal terms of U smaller in magnitude than the value SMALL have been replaced by SMALL, using the same sign as the diagonal terms, and the decomposition has been completed. An approximate solution for X based on this replacement is then obtained if no overflow results and MATINV is filled with nan() value.

If D1/=0 then the linear system MAT * Z = D can be solved with subroutines LU_SOLVE or LU_SOLVE2.

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine chol_cmp ( mat, invdiag, d1, d2, upper, tol )`¶

Purpose¶

CHOL_CMP computes the Cholesky factorization of a n-by-n real symmetric positive definite matrix MAT. The factorization has the form

MAT = U’ * U , if UPPER=true or is absent,

and

MAT = L * L’ , if UPPER=false,

where U is an upper triangular matrix and L is a lower triangular matrix.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric positive definite matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

On exit, if D1/=0, the factor U or L from the Cholesky factorization MAT = U’ * U or MAT = L * L’, except for the main diagonal elements which are stored in reciprocal form in INVDIAG. The main diagonal elements of MAT are not modified.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

INVDIAG (OUTPUT) real(stnd), dimension(:)

On exit, INVDIAG contains the reciprocals of the actual diagonal elements of L or U.

The size of INVDIAG must verify: size( INVDIAG ) = n .

D1 (OUTPUT) real(stnd)

On exit, D1 = zero indicates that the matrix MAT is algorithmically not positive definite and that the factorization can not be completed. Any other value indicates successful exit.

On exit, if D2 is present, D1 is the first component of the determinant (mantissa of determinant) of MAT.

D2 (OUTPUT, OPTIONAL) integer(i4b)

If D2 is present, the components of the determinant of MAT are computed.

On exit, D2 is the second component of the determinant of MAT (exponent of determinant):

determinant(MAT) = scale( D1, D2 )

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test the matrix for positive-definiteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

MAT is declared not positive definite if during the j-th stage of the factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= MAT(j,j) * TOL

In this case, the leading minor of order j of MAT is declared not positive definite and on exit of CHOL_CMP:

D1 is set to zero,

INVDIAG(j) = PIV(j),

INVDIAG(j+1_i4b:n) are set to nan() value,

and the Cholesky factorization is not completed.

On the other hand, if MAT is positive definite then

U(j,j) = sqrt(PIV(j)) (if UPPER=true), for j=1 to n,

or

L(j,j) = sqrt(PIV(j)) (if UPPER=false), for j=1 to n,

and on exit of CHOL_CMP:

D1/=0,

INVDIAG(j)=1/sqrt(PIV(j)) for j=1 to n,

and the linear system MAT * Z = D can be solved with subroutine CHOL_SOLVE.

This is a GAxpy version of the Cholesky algorithm, for more details see the reference (1).

A blocked algorithm is used to compute the factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2009:

Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 1, pp 251-254.

`subroutine chol_cmp2 ( mat, invdiag, d1, d2, b, matinv, upper, fill, tol )`¶

Purpose¶

CHOL_CMP2 computes the Cholesky factorization of a n-by-n real symmetric positive definite matrix MAT. The factorization has the form

MAT = U’ * U , if UPPER=true or is absent,

and

MAT = L * L’ , if UPPER=false,

where U is an upper triangular matrix and L is a lower triangular matrix.

If D2 is present, CHOL_CMP2 computes the determinant of MAT as

determinant(MAT) = scale( D1, D2 )

If B is present, CHOL_CMP2 solves the system of linear equations

MAT * X = B

using the Cholesky factorization of MAT. Here B is a n-vector.

If MATINV is present, CHOL_CMP2 computes the inverse of MAT

MATINV = MAT**(-1)

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

On exit, if D1/=0, the factor U or L from the Cholesky factorization MAT = U’ * U or MAT = L * L’, except for the main diagonal elements which are stored in reciprocal form in INVDIAG. The main diagonal elements of MAT are not modified.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

INVDIAG (OUTPUT) real(stnd), dimension(:)

On exit, INVDIAG contains the reciprocals of the actual diagonal elements of L or U.

The size of INVDIAG must verify: size( INVDIAG ) = n .

D1 (OUTPUT) real(stnd)

On exit, D1 = zero indicates that the matrix MAT is algorithmically not positive definite and that the factorization can not be completed. Any other value indicates successful exit.

On exit, if D2 is present, D1 is the first component of the determinant (mantissa of determinant) of MAT.

D2 (OUTPUT, OPTIONAL) integer(i4b)

If D2 is present, the components of the determinant of MAT are computed.

On exit, D2 is the second component of the determinant of MAT (exponent of determinant)

determinant(MAT) = scale( D1, D2 )

B (INPUT/OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, the right hand side vector B. On exit, the solution vector X.

The shape of B must verify: size( B ) = n .

MATINV (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, if:

FILL = true or is absent: The (symmetric) inverse of MAT.

FILL = false: The upper (if UPPER=true) or lower (if UPPER=false) triangle of the (symmetric) inverse of MAT, is stored in the upper or lower triangular part of the matrix MATINV and the other part of MATINV is not referenced.

The shape of MATINV must verify: size( MATINV, 1 ) = size( MATINV, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of MATINV is filled on exit

FILL= true and UPPER= false, the upper triangle of MATINV is filled on exit

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of MATINV is not filled on exit.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test the matrix for positive-definiteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

MAT is declared not positive definite if during the j-th stage of the factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= MAT(j,j) * TOL

In this case, the leading minor of order j of MAT is declared not positive definite and on exit of CHOL_CMP2:

D1 is set to zero,

INVDIAG(j) = PIV(j),

INVDIAG(j+1_i4b:n) are set to nan() value,

B is filled with nan() value,

the upper or lower triangle of MATINV is filled with nan() value if FILL=false,

the matrix MATINV is filled with nan() value if FILL=true,

and the Cholesky factorization is not completed.

On the other hand, if MAT is positive definite then

U(j,j) = sqrt(PIV(j)) (if UPPER=true), for j=1 to n,

or

L(j,j) = sqrt(PIV(j)) (if UPPER=false), for j=1 to n,

and on exit of CHOL_CMP2:

D1/=0,

INVDIAG(j)=1/sqrt(PIV(j)) for j=1 to n,

B and MATINV are computed, if these arguments are present,

and the linear system MAT * Z = D can be solved with subroutine CHOL_SOLVE.

This is a GAxpy version of the Cholesky algorithm, for more details see the reference (1).

A blocked algorithm is used to compute the factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2009:

Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 1, pp 251-254.

`subroutine gchol_cmp ( mat, invdiag, krank, d1, d2, upper, tol )`¶

Purpose¶

GCHOL_CMP computes the Cholesky factorization of a n-by-n real symmetric positive semidefinite matrix MAT. The factorization has the form

MAT = U’ * U , if UPPER=true or is absent,

and

MAT = L * L’ , if UPPER=false,

where U is an upper triangular matrix and L is a lower triangular matrix.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

On exit, if D1>=0, the factor U or L from the Cholesky factorization MAT = U’ * U or MAT = L * L’, except for the main diagonal elements which are stored in reciprocal form in INVDIAG. The main diagonal elements of MAT are not modified.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

INVDIAG (OUTPUT) real(stnd), dimension(:)

On exit, INVDIAG contains the reciprocals of the actual diagonal elements of L or U, excepted for zeroed elements if MAT is not positive definite.

The shape of INVDIAG must verify: size( INVDIAG ) = n .

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of MAT, which is defined as the number of nonzero elements of the diagonal of L or U. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

D1 (OUTPUT) real(stnd)

On exit, D1 < zero indicates that the matrix MAT is algorithmically not positive semidefinite and that the factorization can not be completed. Any other value indicates successful exit.

On exit, if D2 is present, D1 is the first component of the determinant (mantissa of determinant) of MAT.

D2 (OUTPUT, OPTIONAL) integer(i4b)

If D2 is present, the components of the determinant of MAT are computed.

On exit, D2 is the second component of the determinant of MAT (exponent of determinant)

determinant(MAT) = scale( D1, D2 )

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test MAT for positive-semidefiniteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

MAT is declared not positive semidefinite if during the j-th stage of the factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= - abs( MAT(j,j) * TOL )

In this case, the leading minor of order j of MAT is declared not positive semidefinite and on exit of GCHOL_CMP:

D1 is set to -1,

KRANK is set to -1,

INVDIAG(j) = PIV(j),

INVDIAG(j+1_i4b:n) are set to nan() value,

and the Cholesky factorization is not completed.

On the other hand, if MAT is positive semidefinite (e.g. D1>=0), KRANK is computed as follows:

KRANK is initially set to n. if, during the factorization, a pivot, PIV(j), is such that

abs( PIV(j) ) <= abs( MAT(j,j) * TOL )

KRANK is decreased by 1 and U(j,j:n) (if UPPER=true) or L(j:n,j) (if UPPER=false) is set to zero. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

IF PIV(j) does not satisfy this condition then

U(j,j) = sqrt(PIV(j)) (if UPPER=true),

L(j,j) = sqrt(PIV(j)) (if UPPER=false).

On exit of GCHOL_CMP, if MAT is positive semidefinite, INVDIAG contains the reciprocals of the diagonal elements of U or L, excepted for zeroed elements during the factorization as described above.

If MAT is positive semidefinite (D1>=0), the linear system MAT * Z = D can also be solved with subroutine CHOL_SOLVE.

This is a GAxpy version of the Cholesky algorithm, for more details see the reference (1).

A blocked algorithm is used to compute the factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2009:

Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 1, pp 251-254.

`subroutine gchol_cmp2 ( mat, invdiag, krank, d1, d2, b, matinv, upper, fill, tol )`¶

Purpose¶

GCHOL_CMP2 computes the Cholesky factorization of a n-by-n real symmetric positive semidefinite matrix MAT. The factorization has the form

MAT = U’ * U , if UPPER=true or is absent,

and

MAT = L * L’ , if UPPER=false,

where U is an upper triangular matrix and L is a lower triangular matrix.

If D2 is present, GCHOL_CMP2 computes the determinant of MAT as

determinant(MAT) = scale( D1, D2 )

If B is present, GCHOL_CMP2 solves the system of linear equations

MAT * X = B

using the Cholesky factorization of MAT if B belongs to the range of MAT. Here B is a n-vector. IF B does not belongs to the range of MAT, an approximate solution is computed as

X = MATINV * B

where MATINV is a (generalized) inverse of MAT.

If MATINV is present, GCHOL_CMP2 computes a (generalized) inverse of MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

On exit, if D1>=0, the factor U or L from the Cholesky factorization MAT = U’ * U or MAT = L * L’, except for the main diagonal elements which are stored in reciprocal form in INVDIAG. The main diagonal elements of MAT are not modified.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

INVDIAG (OUTPUT) real(stnd), dimension(:)

On exit, INVDIAG contains the reciprocals of the actual diagonal elements of L or U, excepted for zeroed elements if MAT is not positive definite.

The shape of INVDIAG must verify: size( INVDIAG ) = n .

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of MAT, which is defined as the number of nonzero elements of the diagonal of L or U. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

D1 (OUTPUT) real(stnd)

On exit, D1 < zero indicates that the matrix MAT is algorithmically not positive semidefinite and that the factorization can not be completed. Any other value indicates successful exit.

On exit, if D2 is present, D1 is the first component of the determinant (mantissa of determinant) of MAT.

D2 (OUTPUT, OPTIONAL) integer(i4b)

If D2 is present, the components of the determinant of MAT are computed.

On exit, D2 is the second component of the determinant of MAT (exponent of determinant)

determinant(MAT) = scale( D1, D2 )

B (INPUT/OUTPUT, OPTIONAL) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, one solution vector X if B belongs to the range of MAT, otherwise an approximate solution computed with the help of a symmetric generalized inverse of MAT.

The shape of B must verify: size( B ) = n .

MATINV (OUTPUT) real(stnd), dimension(:,:)

On exit, if:

FILL = true or is absent: The (symmetric) (generalized) inverse of MAT.

FILL = false: The upper (if UPPER=true) or lower (if UPPER=false) triangle of the (symmetric) (generalized) inverse of MAT, is stored in the upper or lower triangular part of the matrix MATINV and the other part of MATINV is not referenced.

The shape of MATINV must verify: size( MATINV, 1 ) = size( MATINV, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of MATINV is filled on exit.

FILL= true and UPPER= false, the upper triangle of MATINV is filled on exit.

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of MATINV is not filled on exit.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test MAT for positive-semidefiniteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default is the machine precision multiplied by n.

Further Details¶

MAT is declared not positive semidefinite if during the j-th stage of the factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= - abs( MAT(j,j) * TOL )

In this case, the leading minor of order j of MAT is declared not positive semidefinite and on exit of GCHOL_CMP2:

D1 is set to -1,

KRANK is set to -1,

INVDIAG(j) = PIV(j),

INVDIAG(j+1_i4b:n) are set to nan() value,

B is filled with nan() value,

the upper or lower triangle of MATINV is filled with nan() value if FILL=false,

the matrix MATINV is filled with nan() value if FILL=true,

and the Cholesky factorization is not completed.

On the other hand, if MAT is positive semidefinite (D1>=0), KRANK is computed as follows:

KRANK is initially set to n. if, during the factorization, a pivot, PIV(j), is such that

abs( PIV(j) ) <= abs( MAT(j,j) * TOL )

KRANK is decreased by 1 and U(j,j:n) (if UPPER=true) or L(j:n,j) (if UPPER=false) is set to zero. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

IF PIV(j) does not satisfy this condition then

U(j,j) = sqrt(PIV(j)) (if UPPER=true),

L(j,j) = sqrt(PIV(j)) (if UPPER=false).

On exit of GCHOL_CMP2, if MAT is positive semidefinite, INVDIAG contains the reciprocals of the diagonal elements of U or L, excepted for zeroed elements during the factorization as described above.

If MAT is positive semidefinite (D1>=0), MATINV is computed as follows. If:

KRANK=n, MATINV is just the inverse of MAT, MATINV = MAT**(-1) ,

KRANK<n, MATINV is a generalized inverse of MAT.

MATINV is a generalized inverse of MAT if

MAT * MATINV * MAT = MAT and MATINV * MAT * MATINV = MATINV

If MAT is positive semidefinite (D1>=0), the linear system MAT * Z = D can also be solved with subroutine CHOL_SOLVE if D belongs to the range of MAT.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2009:

Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 1, pp 251-254.

`subroutine lu_solve ( mat, ip, b )`¶

Purpose¶

LU_SOLVE solves a system of linear equations

MAT * X = B

where MAT is a n-by-n coefficient matrix and B is a n-vector, using the LU factorization with scaled partial pivoting of MAT, P * MAT = L * U, as computed by LU_CMP or LU_CMP2.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the LU factorization of P * MAT for some permutation matrix P specified by argument IP. It is assumed that MAT is as generated by LU_CMP or LU_CMP2.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

IP (INPUT) integer(i4b), dimension(:)

The permutation matrix P as generated by LU_CMP or LU_CMP2.

The shape of IP must verify: size( IP ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, the solution vector X.

The shape of B must verify: size( B ) = n .

Further Details¶

It is assumed that LU_CMP or LU_CMP2 has been used to compute the LU factorization of MAT before LU_SOLVE.

MAT and IP are not modified by this routine and can be left in place for successive calls with different right-hand sides B.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine lu_solve ( mat, ip, b )`¶

Purpose¶

LU_SOLVE solves a system of linear equations with several right hand sides

MAT * X = B

where MAT is a n-by-n coefficient matrix and B is a n-by-nb matrix, using the LU factorization with scaled partial pivoting of MAT, P * MAT = L * U, as computed by LU_CMP or LU_CMP2.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the LU factorization of P * MAT for some permutation matrix P specified by argument IP. It is assumed that MAT is as generated by LU_CMP or LU_CMP2.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

IP (INPUT) integer(i4b), dimension(:)

The permutation matrix P as generated by LU_CMP or LU_CMP2.

The shape of IP must verify: size( IP ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B.

On exit, the solution matrix X.

The shape of B must verify: size( B, 1 ) = n .

Further Details¶

It is assumed that LU_CMP or LU_CMP2 has been used to compute the LU factorization of MAT before LU_SOLVE.

MAT and IP are not modified by this routine and can be left in place for successive calls with different right-hand sides B.

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine lu_solve2 ( mat, ip, b )`¶

Purpose¶

LU_SOLVE2 solves a system of linear equations

MAT * X = B

where MAT is a n-by-n coefficient matrix and B is a n-vector, using the LU factorization with scaled partial pivoting of MAT, P * MAT = L * U, as computed by LU_CMP or LU_CMP2.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the LU factorization of P * MAT for some permutation matrix P specified by argument IP. It is assumed that MAT is as generated by LU_CMP or LU_CMP2.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

IP (INPUT) integer(i4b), dimension(:)

The permutation matrix P as generated by LU_CMP or LU_CMP2.

The shape of IP must verify: size( IP ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, the solution vector X.

The shape of B must verify: size( B ) = n .

Further Details¶

It is assumed that LU_CMP or LU_CMP2 has been used to factor MAT before LU_SOLVE2.

MAT and IP are not modified by this routine and can be left in place for successive calls with different right-hand sides B.

This subroutines takes into account the possibility that B will begin with many zero elements, so it is efficient for use in matrix inversion.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine lu_solve2 ( mat, ip, b )`¶

Purpose¶

LU_SOLVE2 solves a system of linear equations with several right hand sides

MAT * X = B

where MAT is a n-by-n coefficient matrix and B is a n-by-nb matrix, using the LU factorization with scaled partial pivoting of MAT, P * MAT = L * U, as computed by LU_CMP or LU_CMP2.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the LU factorization of P * MAT for some permutation matrix P specified by argument IP. It is assumed that MAT is as generated by LU_CMP or LU_CMP2.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

IP (INPUT) integer(i4b), dimension(:)

The permutation matrix P as generated by LU_CMP or LU_CMP2.

The shape of IP must verify: size( IP ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B.

On exit, the solution matrix X.

The shape of B must verify: size( B, 1 ) = n .

Further Details¶

It is assumed that LU_CMP or LU_CMP2 has been used to factor MAT before LU_SOLVE2.

MAT and IP are not modified by this routine and can be left in place for successive calls with different right-hand sides B.

This subroutines takes into account the possibility that each column of B will begin with many zero elements, so it is efficient for use in matrix inversion.

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`function solve_lin ( mat, b, tol )`¶

Purpose¶

SOLVE_LIN solves a system of linear equations

MAT * X = B

with a n-by-n coefficient matrix MAT. B is a n-vector.

The function returns the solution vector X, if the matrix MAT is not singular.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the coefficient matrix of the equation

MAT * X = B

MAT is not modified by the subroutine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT) real(stnd), dimension(:)

On entry, the right hand side vector B. B is not modified by the subroutine.

The shape of B must verify: size( B ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

Further Details¶

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT is used to solve the linear system.

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular.

On exit, if MAT is singular, the function returns a n-vector filled with nan() value.

A blocked algorithm is used to compute the factorization and solve the triangular systems. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`function solve_lin ( mat, b, tol )`¶

Purpose¶

SOLVE_LIN solves a system of linear equations with several right hand sides

MAT * X = B

with a n-by-n coefficient matrix MAT. B is a n-by-nb matrix.

The function returns the n-by-nb solution matrix X, if MAT is not singular.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the coefficient matrix of the equation

MAT * X = B

MAT is not modified by the subroutine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B. B is not modified by the subroutine.

The shape of B must verify: size( B, 1 ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

Further Details¶

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT is used to solve the linear system.

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular.

On exit, if MAT is algorithmically singular, the function returns a n-by-nb matrix filled with nan() value.

A blocked algorithm is used to compute the factorization and solve the triangular systems. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine lin_lu_solve ( mat, b, failure, tol, small )`¶

Purpose¶

LIN_LU_SOLVE solves a system of linear equations

MAT * X = B

with a n-by-n coefficient matrix MAT. B is a n-vector.

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT

P * MAT = L * U

where P is a permutation matrix, L is a n-by-n unit lower triangular matrix and U is a n-by-n upper triangular matrix, is used to solve the linear system.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, MAT contains the coefficient matrix of the equation

MAT * X = B

On exit, MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, the solution vector X.

The shape of B must verify: size( B ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically singular

FAILURE = false: MAT is not singular.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

SMALL (INPUT, OPTIONAL) real(stnd)

On entry, if the system is singular, replaces a diagonal term of the matrix U if it is smaller in magnitude than the value SMALL using the same sign as the diagonal term. An approximate solution based on this replacement is obtained if no overflow results. If SMALL is supplied as less than SAFMIN, the smallest number that can be reciprocated safely, then the value SAFMIN is used in place of SMALL.

Default: SAFMIN, the smallest number that can be reciprocated safely.

Further Details¶

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and FAILURE is set to true. Otherwise, FAILURE is set to false.

If MAT is algorithmically singular (FAILURE=true), the diagonal terms of U smaller in magnitude than the value SMALL have been replaced by SMALL, using the same sign as the diagonal terms, and the decomposition has been completed. An approximate solution based on this replacement is obtained if no overflow results.

A blocked algorithm is used to compute the factorization and solve the triangular systems. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine lin_lu_solve ( mat, b, failure, tol, small )`¶

Purpose¶

LIN_LU_SOLVE solves a system of linear equations with several right hand sides

MAT * X = B

with a n-by-n coefficient matrix MAT. B is a n-by-nb matrix.

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT

P * MAT = L * U

where P is a permutation matrix, L is a n-by-n unit lower triangular matrix and U is a n-by-n upper triangular matrix, is used to solve the linear systems.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, MAT contains the coefficient matrix of the equation

MAT * X = B

On exit, MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B.

On exit, the solution matrix X.

The shape of B must verify: size( B, 1 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically singular

FAILURE = false: MAT is not singular.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

SMALL (INPUT, OPTIONAL) real(stnd)

On entry, if the system is singular, replaces a diagonal term of the matrix U if it is smaller in magnitude than the value SMALL using the same sign as the diagonal term. Approximate solutions based on this replacement are obtained if no overflow results. If SMALL is supplied as less than SAFMIN, the smallest number that can be reciprocated safely, then the value SAFMIN is used in place of SMALL.

Default: SAFMIN, the smallest number that can be reciprocated safely.

Further Details¶

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and FAILURE is set to true. Otherwise, FAILURE is set to false.

If MAT is algorithmically singular (FAILURE=true), the diagonal terms of U smaller in magnitude than the value SMALL have been replaced by SMALL, using the same sign as the diagonal terms, and the decomposition has been completed. Approximate solutions based on this replacement are obtained if no overflow results.

A blocked algorithm is used to compute the factorization and solve the triangular systems. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine lin_lu_solve ( mat, b, failure, x, tol, small )`¶

Purpose¶

LIN_LU_SOLVE solves a system of linear equations

MAT * X = B

with a n-by-n coefficient matrix MAT. B and X are n-vectors.

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT

P * MAT = L * U

where P is a permutation matrix, L is a n-by-n unit lower triangular matrix and U is a n-by-n upper triangular matrix, is used to solve the linear system.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the coefficient matrix of the equation

MAT * X = B

MAT is not modified by the subroutine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT) real(stnd), dimension(:)

On entry, the right hand side vector B. B is not modified by the subroutine.

The shape of B must verify: size( B ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically singular

FAILURE = false: MAT is not singular.

X (OUTPUT) real(stnd), dimension(:)

On exit, the solution vector X.

The shape of X must verify: size( X ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

SMALL (INPUT, OPTIONAL) real(stnd)

On entry, if the system is singular, replaces a diagonal term of the matrix U if it is smaller in magnitude than the value SMALL using the same sign as the diagonal term. An approximate solution based on this replacement is obtained if no overflow results. If SMALL is supplied as less than SAFMIN, the smallest number that can be reciprocated safely, then the value SAFMIN is used in place of SMALL.

Default: SAFMIN, the smallest number that can be reciprocated safely.

Further Details¶

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and FAILURE is set to true. Otherwise, FAILURE is set to false.

If MAT is algorithmically singular (FAILURE=true), the diagonal terms of U smaller in magnitude than the value SMALL have been replaced by SMALL, using the same sign as the diagonal terms, and the decomposition has been completed. An approximate solution based on this replacement is obtained if no overflow results.

A blocked algorithm is used to compute the factorization and solve the triangular systems. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine lin_lu_solve ( mat, b, failure, x, tol, small )`¶

Purpose¶

LIN_LU_SOLVE solves a system of linear equations with several right hand sides

MAT * X = B

with a n-by-n coefficient matrix MAT. B and X are n-by-nb matrices.

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT

P * MAT = L * U

where P is a permutation matrix, L is a n-by-n unit lower triangular matrix and U is a n-by-n upper triangular matrix, is used to solve the linear systems.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

MAT contains the coefficient matrix of the equation

MAT * X = B

MAT is not modified by the subroutine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B. B is not modified by the subroutine.

The shape of B must verify: size( B, 1 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically singular

FAILURE = false: MAT is not singular.

X (OUTPUT) real(stnd), dimension(:,:)

On exit, the solution matrix X.

The shape of X must verify:

size( X, 1 ) = n ;

size( X, 2 ) = size( B, 2 ) = nb .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

SMALL (INPUT, OPTIONAL) real(stnd)

On entry, if the system is singular, replaces a diagonal term of the matrix U if it is smaller in magnitude than the value SMALL using the same sign as the diagonal term. Approximate solutions based on this replacement are obtained if no overflow results. If SMALL is supplied as less than SAFMIN, the smallest number that can be reciprocated safely, then the value SAFMIN is used in place of SMALL.

Default: SAFMIN, the smallest number that can be reciprocated safely.

Further Details¶

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and FAILURE is set to true. Otherwise, FAILURE is set to false.

If MAT is algorithmically singular (FAILURE=true), the diagonal terms of U smaller in magnitude than the value SMALL have been replaced by SMALL, using the same sign as the diagonal terms, and the decomposition has been completed. Approximate solutions based on this replacement are obtained if no overflow results.

A blocked algorithm is used to compute the factorization and solve the triangular systems. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine chol_solve ( mat, invdiag, b, upper )`¶

Purpose¶

CHOL_SOLVE solves a system of linear equations

MAT * X = B

where MAT is a n-by-n symmetric positive definite matrix and B is a n-vector, using the CHOLESKY factorization MAT = U’ * U or MAT = L * L’, as computed by CHOL_CMP or GCHOL_CMP.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the triangular factor U or L from the Cholesky factorisation, as computed by CHOL_CMP, except for the main diagonal elements which are stored in reciprocal form in INVDIAG.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

INVDIAG (INPUT) real(stnd), dimension(:)

On entry, INVDIAG contains the reciprocals of the actual diagonal elements of L or U, as computed by CHOL_CMP.

The shape of INVDIAG must verify: size( INVDIAG ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, the solution vector X.

The shape of B must verify: size( B ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

This argument must have the same value as used in CHOL_CMP or GCHOL_CMP subroutines for computing the Cholesky factorisation.

The default is true.

`subroutine chol_solve ( mat, invdiag, b, upper )`¶

Purpose¶

CHOL_SOLVE solves a system of linear equations with several right hand sides

MAT * X = B

where MAT is a n-by-n symmetric positive (semi)-definite matrix and B is a n-by-nb matrix, using the CHOLESKY factorization MAT = U’ * U or MAT = L * L’ as computed by CHOL_CMP or GCHOL_CMP.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the triangular factor U or L from the Cholesky factorisation, as computed by CHOL_CMP, except for the main diagonal elements which are stored in reciprocal form in INVDIAG.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

INVDIAG (INPUT) real(stnd), dimension(:)

On entry, INVDIAG contains the reciprocals of the actual diagonal elements of L or U, as computed by CHOL_CMP.

The shape of INVDIAG must verify: size( INVDIAG ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B.

On exit, the solution matrix X.

The shape of B must verify: size( B, 1 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

This argument must have the same value as used in CHOL_CMP or GCHOL_CMP subroutines for computing the Cholesky factorisation.

The default is true.

Further Details¶

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2009:

Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 1, pp 251-254.

`subroutine triang_solve ( mat, b, upper, trans )`¶

Purpose¶

TRIANG_SOLVE solves a triangular system of the form

MAT * X = B or MAT’ * X = B,

where MAT is a triangular matrix of order n, and B is an n-vector.

No test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

The triangular matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of the array MAT contains the upper triangular matrix, and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of the array MAT contains the lower triangular matrix, and the strictly upper triangular part of MAT is not referenced.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, the solution vector X.

The shape of B must verify: size( B ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether MAT is upper or lower triangular. If:

UPPER = true : MAT is upper triangular

UPPER = false: MAT is lower triangular.

The default is true.

TRANS (INPUT, OPTIONAL) logical(lgl)

Specifies the form of the system of equations. If:

TRANS = true : MAT’ * X = B (Transpose)

TRANS = false: MAT * X = B (No transpose)

The default is false.

`subroutine triang_solve ( mat, b, scal, upper, trans )`¶

Purpose¶

TRIANG_SOLVE solves a nonsingular triangular linear system of the form

MAT * X = B or MAT’ * X = B,

where MAT is a triangular matrix of order n, and B is an n-vector.

The matrix MAT is assumed to be ill-conditioned, and frequent rescalings are carried out in order to avoid overflow. However, no test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

The nonsingular triangular matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of the array MAT contains the upper triangular matrix, and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of the array MAT contains the lower triangular matrix, and the strictly upper triangular part of MAT is not referenced.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector B.

On exit, the scaled solution vector X.

The shape of B must verify: size( B ) = n .

SCAL (OUTPUT) real(stnd)

On exit, SCAL is a scaling factor introduced in order to avoid overflow. The solution of the given system of equations is B(:)/SCAL . Note that SCAL may be negative.

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether MAT is upper or lower triangular. If:

UPPER = true : MAT is upper triangular

UPPER = false: MAT is lower triangular.

The default is true.

TRANS (INPUT, OPTIONAL) logical(lgl)

Specifies the form of the system of equations. If:

TRANS = true : MAT’ * X = B (Transpose)

TRANS = false: MAT * X = B (No transpose)

The default is false.

`subroutine triang_solve ( mat, b, upper, trans )`¶

Purpose¶

TRIANG_SOLVE solves a triangular system of the form

MAT * X = B or MAT’ * X = B,

where MAT is a triangular matrix of order n, and B is an n-by-nb matrix.

No test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

The triangular matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of the array MAT contains the upper triangular matrix, and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of the array MAT contains the lower triangular matrix, and the strictly upper triangular part of MAT is not referenced.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B.

On exit, the solution matrix X.

The shape of B must verify: size( B, 1 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether MAT is upper or lower triangular. If:

UPPER = true : MAT is upper triangular

UPPER = false: MAT is lower triangular.

The default is true.

TRANS (INPUT, OPTIONAL) logical(lgl)

Specifies the form of the system of equations. If:

TRANS = true : MAT’ * X = B (Transpose)

TRANS = false: MAT * X = B (No transpose)

The default is false.

Further Details¶

The computations are parallelized if OPENMP is used.

`subroutine triang_solve ( mat, b, scal, upper, trans )`¶

Purpose¶

TRIANG_SOLVE solves a nonsingular triangular linear system of the form

MAT * X = B or MAT’ * X = B,

where MAT is a triangular matrix of order n, and B is an n-by-nb matrix.

The matrix MAT is assumed to be ill-conditioned, and frequent rescalings are carried out in order to avoid overflow. However, no test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

The nonsingular triangular matrix MAT. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of the array MAT contains the upper triangular matrix, and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of the array MAT contains the lower triangular matrix, and the strictly upper triangular part of MAT is not referenced.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

B (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the right hand side matrix B.

On exit, the scaled solution matrix X.

The shape of B must verify:

size( B, 1 ) = n ,

size( B, 2 ) = size( SCAL ) = nb .

SCAL (OUTPUT) real(stnd), dimension(:)

On exit, SCAL is a scaling vector introduced in order to avoid overflow. The solution of the given system of equations is B/spread(SCAL,dim=1, ncopies=n) . Note that elements of SCAL may be negative.

The size of SCAL must verify: size( SCAL ) = size( B, 2 ) = nb .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether MAT is upper or lower triangular. If:

UPPER = true : MAT is upper triangular

UPPER = false: MAT is lower triangular.

The default is true.

TRANS (INPUT, OPTIONAL) logical(lgl)

Specifies the form of the system of equations. If:

TRANS = true : MAT’ * X = B (Transpose)

TRANS = false: MAT * X = B (No transpose)

The default is false.

Further Details¶

The computations are parallelized if OPENMP is used.

`subroutine comp_inv ( mat, failure, tol )`¶

Purpose¶

COMP_INV computes, in place, the inverse of a matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, MAT contains the matrix to be inverted.

On exit, MAT is replaced by its inverse.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically singular.

FAILURE = false: MAT has been inverted.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

Further Details¶

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT is used to compute the inverse.

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and FAILURE is set to true. Otherwise, FAILURE is set to false.

On exit, if FAILURE=true, MAT is filled with nan() value.

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine comp_inv ( mat, failure, matinv, tol )`¶

Purpose¶

COMP_INV computes the inverse of a matrix MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, MAT contains the matrix to be inverted.

On exit, MAT is replaced by the LU decomposition of a rowwise permutation of MAT.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically singular.

FAILURE = false: MAT has been inverted.

MATINV (OUTPUT) real(stnd), dimension(:,:)

On exit, if MAT is not singular, MATINV contains the inverse of MAT.

The shape of MATINV must verify: size( MATINV, 1 ) = size( MATINV, 2 ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

Further Details¶

MAT is modified by COMP_INV.

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT is used to compute the inverse.

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and FAILURE is set to true. Otherwise, FAILURE is set to false.

On exit, if FAILURE=true, MATINV is filled with nan() value.

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`function inv ( mat, tol )`¶

Purpose¶

INV computes the inverse of a real matrix MAT,

MAT * INV(MAT) = I

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, MAT contains the matrix to be inverted. MAT is not modified by the function.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

Further Details¶

MAT is not modified by function INV.

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT is used to compute the inverse.

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular.

On exit, if MAT is algorithmically singular, the function INV returns a n-by-n matrix filled with nan() value.

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine comp_sym_inv ( mat, failure, upper, fill, tol )`¶

Purpose¶

COMP_SYM_INV computes, in place, the inverse of a real symmetric positive definite matrix MAT using the Cholesky factorization MAT = U’ * U or MAT = L * L’.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix to be inverted:

If UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

If UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

On exit:

If FILL = true or is absent: The (symmetric) inverse of MAT overwrites MAT.

If FILL = false: The upper (if UPPER= true) or lower (if UPPER= false) triangle of the (symmetric) inverse of MAT, overwrites the input upper or lower triangular part of MAT and the other part of MAT is not modified.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically not positive definite.

FAILURE = false: MAT has been inverted.

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of MAT is filled on exit.

FILL= true and UPPER= false, the upper triangle of MAT is filled on exit.

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of MAT is not filled and not modified on exit.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test the matrix for positive-definiteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

MAT is declared not positive definite if during the j-th stage of the factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= MAT(j,j) * TOL

In this case, the leading minor of order j of MAT is declared not positive definite, the Cholesky factorization is not completed and, on exit of COMP_SYM_INV, FAILURE is set to true.

On exit, if FAILURE=true:

The upper or lower triangle of MAT is filled with nan() value if FILL=false.

The matrix MAT is filled with nan() value if FILL=true.

A blocked algorithm is used to compute the Cholesky factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine comp_sym_inv ( mat, failure, matinv, upper, fill, tol )`¶

Purpose¶

COMP_SYM_INV computes the inverse of a real symmetric positive definite matrix MAT using the Cholesky factorization MAT = U’ * U or MAT = L * L’.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix to be inverted.

If UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

If UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

MAT is not modified by the subroutine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically not positive definite.

FAILURE = false: MAT has been inverted.

MATINV (OUTPUT) real(stnd), dimension(:,:)

On exit:

If FILL = true or is absent: The (symmetric) inverse of MAT.

If FILL = false: The upper (if UPPER=true) or lower (if UPPER=false) triangle of the (symmetric) inverse of MAT, is stored in the upper or lower triangular part of the matrix MATINV and the other part of MATINV is not modified.

The shape of MATINV must verify: size( MATINV, 1 ) = size( MATINV, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of MATINV is filled on exit.

FILL= true and UPPER= false, the upper triangle of MATINV is filled on exit.

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of MATINV is not filled and not modified on exit.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test the matrix for positive-definiteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

MAT is not modified by COMP_SYM_INV.

MAT is declared not positive definite if during the j-th stage of the factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= MAT(j,j) * TOL

In this case, the leading minor of order j of MAT is declared not positive definite, the Cholesky factorization is not completed and, on exit of COMP_SYM_INV, FAILURE is set to true.

On exit, if FAILURE=true:

The upper or lower triangle of MATINV is filled with nan() value if FILL=false.

The matrix MATINV is filled with nan() value if FILL=true.

A blocked algorithm is used to compute the factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`function sym_inv ( mat, upper, tol )`¶

Purpose¶

SYM_INV computes the inverse of a real symmetric positive definite matrix MAT using the Cholesky factorization MAT = U’ * U or MAT = L * L’.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix to be inverted. If:

If UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

If UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

MAT is not modified by the function.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular

UPPER = false: Lower triangular.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test the matrix for positive-definiteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

MAT is not modified by function SYM_INV.

The (symmetric) inverse of MAT is returned if MAT is positive definite.

MAT is declared not positive definite if during the j-th stage of the Cholesky factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= MAT(j,j) * TOL

In this case, the leading minor of order j of MAT is declared not positive definite.

On exit, if MAT is algorithmically not positive definite, SYM_INV returns a matrix filled with nan() value.

A blocked algorithm is used to compute the factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine comp_sym_ginv ( mat, failure, krank, upper, fill, tol )`¶

Purpose¶

COMP_SYM_GINV computes, in place, the (generalized) inverse of a real symmetric positive semidefinite matrix MAT using the Cholesky factorization MAT = U’ * U or MAT = L * L’.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix to be inverted. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

On exit, if:

FILL = true or is absent: The (symmetric) inverse of MAT overwrites MAT.

FILL = false: The upper or lower triangle of the (symmetric) inverse of MAT, overwrites the input upper or lower triangular part of the matrix MAT and the other part of MAT is not modified.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

FAILURE = true : MAT is algorithmically not positive semidefinite.

FAILURE = false: MAT has been inverted.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of MAT, which is defined as the number of nonzero elements of the diagonal of L or U. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of MAT is filled on exit.

FILL= true and UPPER= false, the upper triangle of MAT is filled on exit.

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of MAT is not filled on exit.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test MAT for positive-semidefiniteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

If MAT is positive semidefinite, the subroutine computes a generalized inverse of MAT. GMAT is a generalized inverse of MAT if

MAT * GMAT * MAT = MAT and GMAT * MAT * GMAT = GMAT

See description of subroutine GCHOL_CMP2 for more details. The subroutine also computes and returns an estimate of the effective rank of MAT in the argument RANK. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

MAT is declared not positive semidefinite if during the j-th stage of the Cholesky factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= - abs( MAT(j,j) * TOL )

In this case, the leading minor of order j of MAT is declared not positive semidefinite, the Cholesky factorization is not completed and on exit of COMP_SYM_GINV, FAILURE is set to true.

On exit, if FAILURE=true:

KRANK is set to -1,

The upper or lower triangle of MAT is filled with nan() value if FILL=false.

The matrix MAT is filled with nan() value if FILL=true.

A blocked algorithm is used to compute the Cholesky factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2009:

Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 1, pp 251-254.

`subroutine comp_sym_ginv ( mat, failure, krank, matinv, upper, fill, tol )`¶

Purpose¶

COMP_SYM_GINV computes the (generalized) inverse of a real symmetric positive semidefinite matrix MAT using the Cholesky factorization MAT = U’ * U or MAT = L * L’.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the symmetric matrix to be inverted. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular part of the matrix and the strictly lower triangular part of MAT is not referenced.

If UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular part of the matrix and the strictly upper triangular part of MAT is not referenced.

MAT is not modified by the subroutine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

FAILURE (OUTPUT) logical(lgl)

On exit, if:

-FAILURE = true : MAT is algorithmically not positive semidefinite. -FAILURE = false: MAT has been inverted.

KRANK (OUTPUT) integer(i4b)

On exit, KRANK contains the effective rank of MAT, which is defined as the number of nonzero elements of the diagonal of L or U. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

MATINV (OUTPUT) real(stnd), dimension(:,:)

On exit, if:

FILL = true or is absent: The (symmetric) (generalized) inverse of MAT.

If FILL = false: The upper (if UPPER=true) or lower (if UPPER=false) triangle of the (symmetric) (generalized) inverse of MAT, is stored in the upper or lower triangular part of the matrix MATINV and the other part of MATINV is not modified.

The shape of MATINV must verify: size( MATINV, 1 ) = size( MATINV, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the upper or lower triangular part of the symmetric matrix MAT is stored. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of MATINV is filled on exit.

FILL= true and UPPER= false, the upper triangle of MATINV is filled on exit.

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of MATINV is not filled on exit.

The default is true.

TOL (INPUT, OPTIONAL) real(stnd)

Tolerance used to test MAT for positive-semidefiniteness. TOL is used as a multiplying factor for determining effective zero for pivots. TOL must be greater or equal to zero, otherwise the default value is used.

The default value is the machine precision multiplied by n.

Further Details¶

MAT is not modified by COMP_SYM_GINV.

If MAT is positive semidefinite, the subroutine computes a generalized inverse of MAT. MATINV is a generalized inverse of MAT if

MAT * MATINV * MAT = MAT and MATINV * MAT * MATINV = MATINV

See description of subroutine GCHOL_CMP2 for more details. The subroutine also computes and returns an estimate of the rank of MAT in the argument RANK. Note that KRANK may be different from the true rank of MAT. See the reference (2) for details.

MAT is declared not positive semidefinite if during the j-th stage of the Cholesky factorization of MAT, a pivot, PIV(j), is such that

PIV(j) <= - abs( MAT(j,j) * TOL )

In this case, the leading minor of order j of MAT is declared not positive semidefinite and on exit of COMP_SYM_GINV, FAILURE is set to true.

On exit, if FAILURE=true:

KRANK is set to -1,

The upper or lower triangle of MATINV is filled with nan() value if FILL=false,

The matrix MATINV is filled with nan() value if FILL=true.

A blocked algorithm is used to compute the Cholesky factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2009:

Cholesky factorization. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 1, pp 251-254.

`subroutine comp_triang_inv ( mat, upper )`¶

Purpose¶

COMP_TRIANG_INV computes, in place, the inverse of a real upper or lower triangular matrix MAT.

On entry, if MAT is algorithmically singular, diagonal terms smaller in magnitude than the value SAFMIN (the smallest number that can be reciprocated safely) are replaced by SAFMIN using the same sign as the diagonal term. An approximate solution based on this replacement is then obtained.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the triangular matrix to be inverted.

If UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular matrix and the strictly lower triangular part of MAT is not referenced.

If UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular matrix and the strictly upper triangular part of MAT is not referenced.

On exit, the (triangular) inverse of the original matrix in the same storage format.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the matrix MAT is upper or lower triangular. If:

UPPER = true : MAT is Upper triangular,

UPPER = false: MAT is Lower triangular.

The default is true.

Further Details¶

The computations are not parallelized even if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine comp_triang_inv ( mat, matinv, upper )`¶

Purpose¶

COMP_TRIANG_INV computes the inverse of a real upper or lower triangular matrix MAT.

On entry, if MAT is algorithmically singular, diagonal terms smaller in magnitude than the value SAFMIN (the smallest number that can be reciprocated safely) are replaced by SAFMIN using the same sign as the diagonal term. An approximate solution based on this replacement is then obtained.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the triangular matrix to be inverted. If:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular matrix and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular matrix and the strictly upper triangular part of MAT is not referenced.

MAT is not modified by the subroutine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

MATINV (OUTPUT) real(stnd), dimension(:,:)

On exit, the (triangular) inverse of the original matrix in the same storage format.

The shape of MATINV must verify: size( MATINV, 1 ) = size( MATINV, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the matrix MAT is upper or lower triangular. If:

UPPER = true : MAT is upper triangular,

UPPER = false: MAT is lower triangular.

The default is true.

Further Details¶

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine comp_uut_ltl ( mat, upper, fill )`¶

Purpose¶

COMP_UUT_LTL computes, in place, the product U * U’ or L’ * L, where the triangular factor U or L is stored in the upper or lower triangular part of MAT.

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, if:

UPPER = true or is absent: The leading n-by-n upper triangular part of MAT contains the upper triangular factor U and the strictly lower triangular part of MAT is not referenced.

UPPER = false: The leading n-by-n lower triangular part of MAT contains the lower triangular factor L and the strictly upper triangular part of MAT is not referenced.

On exit, if:

If FILL = true or is absent: The product U * U’ or L’ * L overwrites MAT.

If FILL = false and UPPER = true or is absent: The leading n-by-n upper triangular part of MAT is overwritten with the upper triangular part of the product U * U’ and the strictly lower triangular part of MAT is not referenced.

If FILL = false and UPPER = false: The leading n-by-n lower triangular part of MAT is overwritten with the lower triangular part of the product L’ * L and the strictly upper triangular part of MAT is not referenced.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the triangular factor stored in matrix MAT is upper or lower triangular. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangularis stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of MAT is filled on exit.

FILL= true and UPPER= false, the upper triangle of MAT is filled on exit.

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of MAT is not filled on exit.

The default is true.

Further Details¶

The computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

`subroutine comp_uut_ltl ( mat, prod, upper, fill )`¶

Purpose¶

COMP_UUT_LTL computes the product U * U’ or L’ * L, where the triangular factor U or L is stored in the upper or lower triangular part of MAT.

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the triangular factor U or L.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

PROD (OUTPUT) real(stnd), dimension(:,:)

On exit, if:

FILL = true or is absent: The product U * U’ or L’ * L .

FILL = false and UPPER = true or is absent: The leading n-by-n upper triangular part of PROD contains the upper triangular part of the product U*Ut and the strictly lower triangular part of PROD is not referenced.

FILL = false and UPPER = false: The leading n-by-n lower triangular part of PROD contains the lower triangular part of the product Lt*L and the strictly upper triangular part of PROD is not referenced.

The shape of PROD must verify: size( PROD, 1 ) = size( PROD, 2 ) = n .

UPPER (INPUT, OPTIONAL) logical(lgl)

Specifies whether the triangular factor stored in matrix MAT is the upper or lower triangle. If:

UPPER = true : Upper triangular is stored

UPPER = false: Lower triangular is stored.

The default is true.

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows. If:

FILL= true and UPPER= true, the lower triangle of PROD is filled on exit.

FILL= true and UPPER= false, the upper triangle of PROD is filled on exit.

FILL= false, the lower (UPPER= true) or upper (UPPER= false) triangle of PROD is not filled on exit.

The default is true.

Further Details¶

The computations are parallelized if OPENMP is used.

`subroutine comp_det ( mat, det, tol, man_det, exp_det )`¶

Purpose¶

COMP_DET computes the determinant of a real matrix MAT

DET = determinant( MAT )

Arguments¶

MAT (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, the real matrix MAT.

On exit, MAT is destroyed.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

DET (OUTPUT) real(stnd)

On exit, the determinant of MAT.

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

MAN_DET (OUTPUT, OPTIONAL) real(stnd)

On exit, the mantissa of the determinant of MAT.

EXP_DET (OUTPUT, OPTIONAL) integer(i4b)

On exit, the exponent of the determinant of MAT.

Further Details¶

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT is used to compute the determinant.

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and DET is set to zero.

On exit:

DET = nan() if MAT is a zero sized matrix.

DET = scale( MAN_DET, EXP_DET ) if minexponent(DET) <= EXP_DET <= maxexponent(DET)

DET = sign( 0, MAN_DET ) if EXP_DET < minexponent(DET)

DET = sign( huge(DET), MAN_DET ) if maxexponent(DET) < EXP_DET

A blocked algorithm is used to compute the LU factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`function det ( mat, tol )`¶

Purpose¶

DET computes the determinant of a real matrix MAT

determinant( MAT ) = DET( MAT )

Arguments¶

MAT (INPUT) real(stnd), dimension(:,:)

On entry, the real matrix MAT. MAT is not modified by the routine.

The shape of MAT must verify: size( MAT, 1 ) = size( MAT, 2 ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not MAT is nearly singular. Tol should normally be choose as approximately the largest relative error in the elements of MAT. For example, if the elements of MAT are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Default: EPS, the relative machine precision.

Further Details¶

The LU decomposition with partial pivoting and implicit row scaling of the matrix MAT is used to compute the determinant.

MAT is declared singular if a diagonal element of U is such that

abs( U(j,j) ) <= n * norm( MAT(j,:) ) * TOL

where norm( MAT(j,:) ) denotes the maximum of the absolute values of the jth row of the matrix MAT. In this case, a diagonal element of U is small, indicating that MAT is singular or nearly singular and DET returns the value zero.

If MAT is a zero sized matrix, DET( MAT ) = nan().

A blocked algorithm is used to compute the LU factorization. Furthermore, the computations are parallelized if OPENMP is used.

For further details, see:

Golub, G.H., and Van Loan, C.F., 1996:

Matrix Computations. 3rd ed. The Johns Hopkins University Press, Baltimore.

Higham, N.J., 2011:

Gaussian elimination. Wiley Interdisciplinary Reviews: Computational Statistics, Vol. 3, Issue 3, pp 230-238.

`subroutine sym_trid_cmp ( d, e, sub, diag, sup1, sup2, perm, tol )`¶

Purpose¶

SYM_TRID_CMP factorizes an n by n symmetric tridiagonal matrix T as

T = P * L * U

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The factorization is obtained by Gaussian elimination with partial pivoting and implicit row scaling.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary .

The size of E must be size( E ) = size( D ) = n .

SUB (OUTPUT) real(stnd), dimension(:)

On exit, SUB(:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L of the factorization of T, SUB(n) is arbitrary .

The size of SUB must verify: size( SUB ) = size( D ) = n .

DIAG (OUTPUT) real(stnd), dimension(:)

On exit, DIAG(:) contains the n diagonal elements of the upper triangular matrix U of the factorization of T.

The size of DIAG must verify: size( DIAG ) = size( D ) = n .

SUP1 (OUTPUT) real(stnd), dimension(:)

On exit, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U of the factorization of T, SUP1(n) is arbitrary .

The size of SUP1 must verify: size( SUP1 ) = size( D ) = n .

SUP2 (OUTPUT) real(stnd), dimension(:)

On exit, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U of the factorization of T, SUP2(n-1:n) is arbitrary .

The size of SUP2 must verify: size( SUP2 ) = size( D ) = n .

PERM (OUTPUT) logical(lgl), dimension(:)

On exit, PERM(:n-1) contains details of the permutation matrix P(j):

if, an interchange occurred at the kth step of the elimination in the factorization of T, then PERM(k) = TRUE

otherwise, PERM(k) = FALSE.

The element PERM(n) is set to TRUE, if there is an integer l such that

abs( U(l,l) ) <= norm( T(l) ) * TOL,

where norm( T(l) ) denotes the sum of the absolute values of the lth row of the matrix T. If no such l exists then PERM(n) is returned as FALSE.

If PERM(n) is returned as TRUE, then a diagonal element of U is small, indicating that T is singular or nearly singular.

The size of PERM must verify: size( PERM ) = size( D ) = n .

TOL (INPUT, OPTIONAL) real(stnd)

On entry, a relative tolerance used to indicate whether or not the matrix T is nearly singular. TOL should normally be choose as approximately the largest relative error in the elements of T. For example, if the elements of T are correct to about 4 significant figures, then TOL should be set to about 5 * 10**(-4). If TOL is supplied as less than EPS, where EPS is the relative machine precision, then the value EPS is used in place of TOL.

Further Details¶

This subroutine is adapted from the subroutine DLAGTF in LAPACK.

`subroutine sym_trid_cmp2 ( d, e, sub, diag, sup1, sup2, perm )`¶

Purpose¶

SYM_TRID_CMP2 factorizes an n by n symmetric tridiagonal matrix T, as

T = P * L * U

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The factorization is obtained by Gaussian elimination with partial pivoting and row interchanges.

Arguments¶

D (INPUT) real(stnd), dimension(:)

On entry, the diagonal elements of the symmetric tridiagonal matrix T.

E (INPUT) real(stnd), dimension(:)

On entry, the n-1 subdiagonal elements of the symmetric tridiagonal matrix T and E(n) is arbitrary .

The size of E must be size( E ) = size( D ) = n .

SUB (OUTPUT) real(stnd), dimension(:)

On exit, SUB(:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L of the factorization of T, SUB(n) is arbitrary .

The size of SUB must verify: size( SUB ) = size( D ) = n .

DIAG (OUTPUT) real(stnd), dimension(:)

On exit, DIAG(:) contains the n diagonal elements of the upper triangular matrix U of the factorization of T.

The size of DIAG must verify: size( DIAG ) = size( D ) = n .

SUP1 (OUTPUT) real(stnd), dimension(:)

On exit, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U of the factorization of T, SUP1(n) is arbitrary .

The size of SUP1 must verify: size( SUP1 ) = size( D ) = n .

SUP2 (OUTPUT) real(stnd), dimension(:)

On exit, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U of the factorization of T, SUP2(n-1:n) is arbitrary .

The size of SUP2 must verify: size( SUP2 ) = size( D ) = n .

PERM (OUTPUT) logical(lgl), dimension(:)

On exit, PERM(:n-1) contains details of the permutation matrix P(j):

If an interchange occurred at the kth step of the elimination in the factorization of T, then PERM(k) = TRUE

Otherwise PERM(k) = FALSE.

PERM(n) is arbitrary .

The size of PERM must verify: size( PERM ) = size( D ) = n .

Further Details¶

SYM_TRID_CMP2 is a simplified version of SYM_TRID_CMP. This subroutine is adapted from the subroutine DGTTRF in LAPACK.

`subroutine sym_trid_solve ( sub, diag, sup1, sup2, perm, y, scale )`¶

Purpose¶

SYM_TRID_SOLVE may be used to solve the system of equations

x(:) * T = scale * y(:)

, where T is an n by n symmetric tridiagonal matrix and scale is a scalar for x(:), following the factorization of T by SYM_TRID_CMP or SYM_TRID_CMP2 as

T = P * L * U

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The matrix T is assumed to be ill-conditioned, and frequent rescalings are carried out in order to avoid overflow. However, no test for singularity or near-singularity is included in this routine. Such tests must be performed before calling this routine.

Arguments¶

SUB (INPUT) real(stnd), dimension(:)

On entry, SUB(:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L of the factorization of T,

SUB(n) is arbitrary.

The size of SUB must verify: size( SUB ) = size( Y ) = n .

DIAG (INPUT) real(stnd), dimension(:)

On entry, DIAG(:) contains the n diagonal elements of the upper triangular matrix U of the factorization of T.

The shape of DIAG must verify: size( DIAG ) = size( Y ) = n .

SUP1 (INPUT) real(stnd), dimension(:)

On entry, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U of the factorization of T,

SUP1(n) is arbitrary.

The shape of SUP1 must verify: size( SUP1 ) = size( Y ) = n .

SUP2 (INPUT) real(stnd), dimension(:)

On entry, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U of the factorization of T, SUP2(n-1:n) is arbitrary.

The shape of SUP2 must verify: size( SUP2 ) = size( Y ) = n .

PERM (INPUT) logical(lgl), dimension(:)

On entry, PERM(:n-1) contains details of the permutation matrix P:

if an interchange occurred at the kth step of the elimination in the factorization of T, then PERM(k) = TRUE,

otherwise PERM(k) = FALSE.

PERM(n) is arbitrary .

The shape of PERM must verify: size( PERM ) = size( Y ) = n .

Y (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector y.

On exit, Y is overwritten the solution vector x.

The shape of Y must verify: size( Y ) = n .

SCALE (OUTPUT) real(stnd)

On exit, the scalar SCALE.

Further Details¶

This subroutine is adapted from the subroutine DLAGTS in LAPACK.

`subroutine sym_trid_solve ( sub, diag, sup1, sup2, perm, y )`¶

Purpose¶

SYM_TRID_SOLVE may be used to solve the system of equations

x(:) * T = y(:)

, where T is an n by n symmetric tridiagonal matrix, following the factorization of T by SYM_TRID_CMP or SYM_TRID_CMP2 as

T = P * L * U

where P is a permutation matrix, L is a unit lower tridiagonal matrix with at most one non-zero sub-diagonal elements per column and U is an upper triangular matrix with at most two non-zero super-diagonal elements per column.

The matrix T is assumed to be no singular and well-conditioned.

Arguments¶

SUB (INPUT) real(stnd), dimension(:)

On entry, SUB(:n-1) contains the n-1 subdiagonal elements of the lower triangular matrix L of the factorization of T, SUB(n) is arbitrary .

The size of SUB must verify: size( SUB ) = size( Y ) = n .

DIAG (INPUT) real(stnd), dimension(:)

On entry, DIAG(:) contains the n diagonal elements of the upper triangular matrix U of the factorization of T.

The shape of DIAG must verify: size( DIAG ) = size( Y ) = n .

SUP1 (INPUT) real(stnd), dimension(:)

On entry, SUP1(:n-1) contains the n-1 superdiagonal elements of the upper triangular matrix U of the factorization of T, SUP1(n) is arbitrary.

The shape of SUP1 must verify: size( SUP1 ) = size( Y ) = n .

SUP2 (INPUT) real(stnd), dimension(:)

On entry, SUP2(:n-2) contains the n-2 second superdiagonal elements of the upper triangular matrix U of the factorization of T, SUP2(n-1:n) is arbitrary.

The shape of SUP2 must verify: size( SUP2 ) = size( Y ) = n .

PERM (INPUT) logical(lgl), dimension(:)

On entry, PERM(:n-1) contains details of the permutation matrix P:

if an interchange occurred at the kth step of the elimination in the factorization of T, then PERM(k) = TRUE,

otherwise PERM(k) = FALSE.

PERM(n) is arbitrary .

The shape of PERM must verify: size( PERM ) = size( Y ) = n .

Y (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the right hand side vector y.

On exit, Y is overwritten the solution vector x.

The shape of Y must verify: size( Y ) = n .

Further Details¶

This subroutine is adapted from the routine DLAGTS in LAPACK.

Module_Lin_Procedures¶

subroutine lu_cmp ( mat, ip, d1, d2, tol, small )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lu_cmp2 ( mat, ip, d1, d2, b, matinv, tol, small )¶

Purpose¶

Arguments¶

Further Details¶

subroutine chol_cmp ( mat, invdiag, d1, d2, upper, tol )¶

Purpose¶

Arguments¶

Further Details¶

subroutine chol_cmp2 ( mat, invdiag, d1, d2, b, matinv, upper, fill, tol )¶

Purpose¶

Arguments¶

Further Details¶

subroutine gchol_cmp ( mat, invdiag, krank, d1, d2, upper, tol )¶

Purpose¶

Arguments¶

Further Details¶

subroutine gchol_cmp2 ( mat, invdiag, krank, d1, d2, b, matinv, upper, fill, tol )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lu_solve ( mat, ip, b )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lu_solve ( mat, ip, b )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lu_solve2 ( mat, ip, b )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lu_solve2 ( mat, ip, b )¶

Purpose¶

Arguments¶

Further Details¶

function solve_lin ( mat, b, tol )¶

Purpose¶

Arguments¶

Further Details¶

function solve_lin ( mat, b, tol )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lin_lu_solve ( mat, b, failure, tol, small )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lin_lu_solve ( mat, b, failure, tol, small )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lin_lu_solve ( mat, b, failure, x, tol, small )¶

Purpose¶

Arguments¶

Further Details¶

subroutine lin_lu_solve ( mat, b, failure, x, tol, small )¶

Purpose¶

Arguments¶

Further Details¶

subroutine chol_solve ( mat, invdiag, b, upper )¶

Purpose¶

Arguments¶

subroutine chol_solve ( mat, invdiag, b, upper )¶

Purpose¶

Arguments¶

Further Details¶

subroutine triang_solve ( mat, b, upper, trans )¶

Purpose¶

Arguments¶

subroutine triang_solve ( mat, b, scal, upper, trans )¶

Purpose¶

Arguments¶

subroutine triang_solve ( mat, b, upper, trans )¶

Purpose¶

`subroutine lu_cmp ( mat, ip, d1, d2, tol, small )`¶

`subroutine lu_cmp2 ( mat, ip, d1, d2, b, matinv, tol, small )`¶

`subroutine chol_cmp ( mat, invdiag, d1, d2, upper, tol )`¶

`subroutine chol_cmp2 ( mat, invdiag, d1, d2, b, matinv, upper, fill, tol )`¶

`subroutine gchol_cmp ( mat, invdiag, krank, d1, d2, upper, tol )`¶

`subroutine gchol_cmp2 ( mat, invdiag, krank, d1, d2, b, matinv, upper, fill, tol )`¶

`subroutine lu_solve ( mat, ip, b )`¶

`subroutine lu_solve ( mat, ip, b )`¶

`subroutine lu_solve2 ( mat, ip, b )`¶

`subroutine lu_solve2 ( mat, ip, b )`¶

`function solve_lin ( mat, b, tol )`¶

`function solve_lin ( mat, b, tol )`¶

`subroutine lin_lu_solve ( mat, b, failure, tol, small )`¶

`subroutine lin_lu_solve ( mat, b, failure, tol, small )`¶

`subroutine lin_lu_solve ( mat, b, failure, x, tol, small )`¶

`subroutine lin_lu_solve ( mat, b, failure, x, tol, small )`¶

`subroutine chol_solve ( mat, invdiag, b, upper )`¶

`subroutine chol_solve ( mat, invdiag, b, upper )`¶

`subroutine triang_solve ( mat, b, upper, trans )`¶

`subroutine triang_solve ( mat, b, scal, upper, trans )`¶

`subroutine triang_solve ( mat, b, upper, trans )`¶

`subroutine triang_solve ( mat, b, scal, upper, trans )`¶

`subroutine comp_inv ( mat, failure, tol )`¶

`subroutine comp_inv ( mat, failure, matinv, tol )`¶

`function inv ( mat, tol )`¶

`subroutine comp_sym_inv ( mat, failure, upper, fill, tol )`¶

`subroutine comp_sym_inv ( mat, failure, matinv, upper, fill, tol )`¶

`function sym_inv ( mat, upper, tol )`¶

`subroutine comp_sym_ginv ( mat, failure, krank, upper, fill, tol )`¶

`subroutine comp_sym_ginv ( mat, failure, krank, matinv, upper, fill, tol )`¶

`subroutine comp_triang_inv ( mat, upper )`¶

`subroutine comp_triang_inv ( mat, matinv, upper )`¶

`subroutine comp_uut_ltl ( mat, upper, fill )`¶

`subroutine comp_uut_ltl ( mat, prod, upper, fill )`¶

`subroutine comp_det ( mat, det, tol, man_det, exp_det )`¶

`function det ( mat, tol )`¶

`subroutine sym_trid_cmp ( d, e, sub, diag, sup1, sup2, perm, tol )`¶