mirror of https://gitlab.com/QEF/q-e.git
3212 lines
92 KiB
Fortran
3212 lines
92 KiB
Fortran
c
|
|
c This file contains LAPACK routines used in quantum-espresso
|
|
c that are part of ATLAS - from www.netlib.org
|
|
c These are:
|
|
* [S,D,C,Z]GESV
|
|
* [S,D,C,Z]GETRF
|
|
* [S,D,C,Z]GETRS
|
|
* [S,D,C,Z]GETRI
|
|
* [S,D,C,Z]TRTRI
|
|
* [S,D,C,Z]POSV
|
|
* [S,D,C,Z]POTRF
|
|
* [S,D,C,Z]POTRS
|
|
* [S,D,C,Z]POTRI
|
|
* [S,D,C,Z]LAUUM
|
|
c
|
|
SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
|
|
*
|
|
* -- LAPACK driver routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* March 31, 1993
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, LDB, N, NRHS
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DGESV computes the solution to a real system of linear equations
|
|
* A * X = B,
|
|
* where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
|
|
*
|
|
* The LU decomposition with partial pivoting and row interchanges is
|
|
* used to factor A as
|
|
* A = P * L * U,
|
|
* where P is a permutation matrix, L is unit lower triangular, and U is
|
|
* upper triangular. The factored form of A is then used to solve the
|
|
* system of equations A * X = B.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* N (input) INTEGER
|
|
* The number of linear equations, i.e., the order of the
|
|
* matrix A. N >= 0.
|
|
*
|
|
* NRHS (input) INTEGER
|
|
* The number of right hand sides, i.e., the number of columns
|
|
* of the matrix B. NRHS >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the N-by-N coefficient matrix A.
|
|
* On exit, the factors L and U from the factorization
|
|
* A = P*L*U; the unit diagonal elements of L are not stored.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* IPIV (output) INTEGER array, dimension (N)
|
|
* The pivot indices that define the permutation matrix P;
|
|
* row i of the matrix was interchanged with row IPIV(i).
|
|
*
|
|
* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
|
|
* On entry, the N-by-NRHS matrix of right hand side matrix B.
|
|
* On exit, if INFO = 0, the N-by-NRHS solution matrix X.
|
|
*
|
|
* LDB (input) INTEGER
|
|
* The leading dimension of the array B. LDB >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, U(i,i) is exactly zero. The factorization
|
|
* has been completed, but the factor U is exactly
|
|
* singular, so the solution could not be computed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. External Subroutines ..
|
|
EXTERNAL DGETRF, DGETRS, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
IF( N.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( NRHS.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -4
|
|
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
|
|
INFO = -7
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DGESV ', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Compute the LU factorization of A.
|
|
*
|
|
CALL DGETRF( N, N, A, LDA, IPIV, INFO )
|
|
IF( INFO.EQ.0 ) THEN
|
|
*
|
|
* Solve the system A*X = B, overwriting B with X.
|
|
*
|
|
CALL DGETRS( 'No transpose', N, NRHS, A, LDA, IPIV, B, LDB,
|
|
$ INFO )
|
|
END IF
|
|
RETURN
|
|
*
|
|
* End of DGESV
|
|
*
|
|
END
|
|
SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* June 30, 1992
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, M, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
DOUBLE PRECISION A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DGETF2 computes an LU factorization of a general m-by-n matrix A
|
|
* using partial pivoting with row interchanges.
|
|
*
|
|
* The factorization has the form
|
|
* A = P * L * U
|
|
* where P is a permutation matrix, L is lower triangular with unit
|
|
* diagonal elements (lower trapezoidal if m > n), and U is upper
|
|
* triangular (upper trapezoidal if m < n).
|
|
*
|
|
* This is the right-looking Level 2 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* M (input) INTEGER
|
|
* The number of rows of the matrix A. M >= 0.
|
|
*
|
|
* N (input) INTEGER
|
|
* The number of columns of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the m by n matrix to be factored.
|
|
* On exit, the factors L and U from the factorization
|
|
* A = P*L*U; the unit diagonal elements of L are not stored.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,M).
|
|
*
|
|
* IPIV (output) INTEGER array, dimension (min(M,N))
|
|
* The pivot indices; for 1 <= i <= min(M,N), row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -k, the k-th argument had an illegal value
|
|
* > 0: if INFO = k, U(k,k) is exactly zero. The factorization
|
|
* has been completed, but the factor U is exactly
|
|
* singular, and division by zero will occur if it is used
|
|
* to solve a system of equations.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE, ZERO
|
|
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER J, JP
|
|
* ..
|
|
* .. External Functions ..
|
|
INTEGER IDAMAX
|
|
EXTERNAL IDAMAX
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DGER, DSCAL, DSWAP, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
IF( M.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DGETF2', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( M.EQ.0 .OR. N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
DO 10 J = 1, MIN( M, N )
|
|
*
|
|
* Find pivot and test for singularity.
|
|
*
|
|
JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 )
|
|
IPIV( J ) = JP
|
|
IF( A( JP, J ).NE.ZERO ) THEN
|
|
*
|
|
* Apply the interchange to columns 1:N.
|
|
*
|
|
IF( JP.NE.J )
|
|
$ CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
|
|
*
|
|
* Compute elements J+1:M of J-th column.
|
|
*
|
|
IF( J.LT.M )
|
|
$ CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
|
|
*
|
|
ELSE IF( INFO.EQ.0 ) THEN
|
|
*
|
|
INFO = J
|
|
END IF
|
|
*
|
|
IF( J.LT.MIN( M, N ) ) THEN
|
|
*
|
|
* Update trailing submatrix.
|
|
*
|
|
CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA,
|
|
$ A( J+1, J+1 ), LDA )
|
|
END IF
|
|
10 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of DGETF2
|
|
*
|
|
END
|
|
SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* March 31, 1993
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, M, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
DOUBLE PRECISION A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DGETRF computes an LU factorization of a general M-by-N matrix A
|
|
* using partial pivoting with row interchanges.
|
|
*
|
|
* The factorization has the form
|
|
* A = P * L * U
|
|
* where P is a permutation matrix, L is lower triangular with unit
|
|
* diagonal elements (lower trapezoidal if m > n), and U is upper
|
|
* triangular (upper trapezoidal if m < n).
|
|
*
|
|
* This is the right-looking Level 3 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* M (input) INTEGER
|
|
* The number of rows of the matrix A. M >= 0.
|
|
*
|
|
* N (input) INTEGER
|
|
* The number of columns of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the M-by-N matrix to be factored.
|
|
* On exit, the factors L and U from the factorization
|
|
* A = P*L*U; the unit diagonal elements of L are not stored.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,M).
|
|
*
|
|
* IPIV (output) INTEGER array, dimension (min(M,N))
|
|
* The pivot indices; for 1 <= i <= min(M,N), row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, U(i,i) is exactly zero. The factorization
|
|
* has been completed, but the factor U is exactly
|
|
* singular, and division by zero will occur if it is used
|
|
* to solve a system of equations.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE
|
|
PARAMETER ( ONE = 1.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER I, IINFO, J, JB, NB
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DGEMM, DGETF2, DLASWP, DTRSM, XERBLA
|
|
* ..
|
|
* .. External Functions ..
|
|
INTEGER ILAENV
|
|
EXTERNAL ILAENV
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
IF( M.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DGETRF', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( M.EQ.0 .OR. N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Determine the block size for this environment.
|
|
*
|
|
NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 )
|
|
IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
|
|
*
|
|
* Use unblocked code.
|
|
*
|
|
CALL DGETF2( M, N, A, LDA, IPIV, INFO )
|
|
ELSE
|
|
*
|
|
* Use blocked code.
|
|
*
|
|
DO 20 J = 1, MIN( M, N ), NB
|
|
JB = MIN( MIN( M, N )-J+1, NB )
|
|
*
|
|
* Factor diagonal and subdiagonal blocks and test for exact
|
|
* singularity.
|
|
*
|
|
CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO )
|
|
*
|
|
* Adjust INFO and the pivot indices.
|
|
*
|
|
IF( INFO.EQ.0 .AND. IINFO.GT.0 )
|
|
$ INFO = IINFO + J - 1
|
|
DO 10 I = J, MIN( M, J+JB-1 )
|
|
IPIV( I ) = J - 1 + IPIV( I )
|
|
10 CONTINUE
|
|
*
|
|
* Apply interchanges to columns 1:J-1.
|
|
*
|
|
CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
|
|
*
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Apply interchanges to columns J+JB:N.
|
|
*
|
|
CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1,
|
|
$ IPIV, 1 )
|
|
*
|
|
* Compute block row of U.
|
|
*
|
|
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
|
|
$ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
|
|
$ LDA )
|
|
IF( J+JB.LE.M ) THEN
|
|
*
|
|
* Update trailing submatrix.
|
|
*
|
|
CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1,
|
|
$ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
|
|
$ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
|
|
$ LDA )
|
|
END IF
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
RETURN
|
|
*
|
|
* End of DGETRF
|
|
*
|
|
END
|
|
SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* June 30, 1999
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, LWORK, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
DOUBLE PRECISION A( LDA, * ), WORK( * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DGETRI computes the inverse of a matrix using the LU factorization
|
|
* computed by DGETRF.
|
|
*
|
|
* This method inverts U and then computes inv(A) by solving the system
|
|
* inv(A)*L = inv(U) for inv(A).
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the factors L and U from the factorization
|
|
* A = P*L*U as computed by DGETRF.
|
|
* On exit, if INFO = 0, the inverse of the original matrix A.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* IPIV (input) INTEGER array, dimension (N)
|
|
* The pivot indices from DGETRF; for 1<=i<=N, row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
|
|
* On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
|
|
*
|
|
* LWORK (input) INTEGER
|
|
* The dimension of the array WORK. LWORK >= max(1,N).
|
|
* For optimal performance LWORK >= N*NB, where NB is
|
|
* the optimal blocksize returned by ILAENV.
|
|
*
|
|
* If LWORK = -1, then a workspace query is assumed; the routine
|
|
* only calculates the optimal size of the WORK array, returns
|
|
* this value as the first entry of the WORK array, and no error
|
|
* message related to LWORK is issued by XERBLA.
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
|
|
* singular and its inverse could not be computed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ZERO, ONE
|
|
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL LQUERY
|
|
INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
|
|
$ NBMIN, NN
|
|
* ..
|
|
* .. External Functions ..
|
|
INTEGER ILAENV
|
|
EXTERNAL ILAENV
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 )
|
|
LWKOPT = N*NB
|
|
WORK( 1 ) = LWKOPT
|
|
LQUERY = ( LWORK.EQ.-1 )
|
|
IF( N.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -3
|
|
ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
|
|
INFO = -6
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DGETRI', -INFO )
|
|
RETURN
|
|
ELSE IF( LQUERY ) THEN
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Form inv(U). If INFO > 0 from DTRTRI, then U is singular,
|
|
* and the inverse is not computed.
|
|
*
|
|
CALL DTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
|
|
IF( INFO.GT.0 )
|
|
$ RETURN
|
|
*
|
|
NBMIN = 2
|
|
LDWORK = N
|
|
IF( NB.GT.1 .AND. NB.LT.N ) THEN
|
|
IWS = MAX( LDWORK*NB, 1 )
|
|
IF( LWORK.LT.IWS ) THEN
|
|
NB = LWORK / LDWORK
|
|
NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) )
|
|
END IF
|
|
ELSE
|
|
IWS = N
|
|
END IF
|
|
*
|
|
* Solve the equation inv(A)*L = inv(U) for inv(A).
|
|
*
|
|
IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
|
|
*
|
|
* Use unblocked code.
|
|
*
|
|
DO 20 J = N, 1, -1
|
|
*
|
|
* Copy current column of L to WORK and replace with zeros.
|
|
*
|
|
DO 10 I = J + 1, N
|
|
WORK( I ) = A( I, J )
|
|
A( I, J ) = ZERO
|
|
10 CONTINUE
|
|
*
|
|
* Compute current column of inv(A).
|
|
*
|
|
IF( J.LT.N )
|
|
$ CALL DGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
|
|
$ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
|
|
20 CONTINUE
|
|
ELSE
|
|
*
|
|
* Use blocked code.
|
|
*
|
|
NN = ( ( N-1 ) / NB )*NB + 1
|
|
DO 50 J = NN, 1, -NB
|
|
JB = MIN( NB, N-J+1 )
|
|
*
|
|
* Copy current block column of L to WORK and replace with
|
|
* zeros.
|
|
*
|
|
DO 40 JJ = J, J + JB - 1
|
|
DO 30 I = JJ + 1, N
|
|
WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
|
|
A( I, JJ ) = ZERO
|
|
30 CONTINUE
|
|
40 CONTINUE
|
|
*
|
|
* Compute current block column of inv(A).
|
|
*
|
|
IF( J+JB.LE.N )
|
|
$ CALL DGEMM( 'No transpose', 'No transpose', N, JB,
|
|
$ N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
|
|
$ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
|
|
CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
|
|
$ ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
|
|
50 CONTINUE
|
|
END IF
|
|
*
|
|
* Apply column interchanges.
|
|
*
|
|
DO 60 J = N - 1, 1, -1
|
|
JP = IPIV( J )
|
|
IF( JP.NE.J )
|
|
$ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
|
|
60 CONTINUE
|
|
*
|
|
WORK( 1 ) = IWS
|
|
RETURN
|
|
*
|
|
* End of DGETRI
|
|
*
|
|
END
|
|
SUBROUTINE DGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* March 31, 1993
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER TRANS
|
|
INTEGER INFO, LDA, LDB, N, NRHS
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DGETRS solves a system of linear equations
|
|
* A * X = B or A' * X = B
|
|
* with a general N-by-N matrix A using the LU factorization computed
|
|
* by DGETRF.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* TRANS (input) CHARACTER*1
|
|
* Specifies the form of the system of equations:
|
|
* = 'N': A * X = B (No transpose)
|
|
* = 'T': A'* X = B (Transpose)
|
|
* = 'C': A'* X = B (Conjugate transpose = Transpose)
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* NRHS (input) INTEGER
|
|
* The number of right hand sides, i.e., the number of columns
|
|
* of the matrix B. NRHS >= 0.
|
|
*
|
|
* A (input) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* The factors L and U from the factorization A = P*L*U
|
|
* as computed by DGETRF.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* IPIV (input) INTEGER array, dimension (N)
|
|
* The pivot indices from DGETRF; for 1<=i<=N, row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
|
|
* On entry, the right hand side matrix B.
|
|
* On exit, the solution matrix X.
|
|
*
|
|
* LDB (input) INTEGER
|
|
* The leading dimension of the array B. LDB >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE
|
|
PARAMETER ( ONE = 1.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL NOTRAN
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
EXTERNAL LSAME
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DLASWP, DTRSM, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
NOTRAN = LSAME( TRANS, 'N' )
|
|
IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
|
|
$ LSAME( TRANS, 'C' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( NRHS.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
|
|
INFO = -8
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DGETRS', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 .OR. NRHS.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
IF( NOTRAN ) THEN
|
|
*
|
|
* Solve A * X = B.
|
|
*
|
|
* Apply row interchanges to the right hand sides.
|
|
*
|
|
CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, 1 )
|
|
*
|
|
* Solve L*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS,
|
|
$ ONE, A, LDA, B, LDB )
|
|
*
|
|
* Solve U*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
|
|
$ NRHS, ONE, A, LDA, B, LDB )
|
|
ELSE
|
|
*
|
|
* Solve A' * X = B.
|
|
*
|
|
* Solve U'*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS,
|
|
$ ONE, A, LDA, B, LDB )
|
|
*
|
|
* Solve L'*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Unit', N, NRHS, ONE,
|
|
$ A, LDA, B, LDB )
|
|
*
|
|
* Apply row interchanges to the solution vectors.
|
|
*
|
|
CALL DLASWP( NRHS, B, LDB, 1, N, IPIV, -1 )
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of DGETRS
|
|
*
|
|
END
|
|
|
|
SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX )
|
|
*
|
|
* -- LAPACK auxiliary routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* June 30, 1999
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INCX, K1, K2, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
DOUBLE PRECISION A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DLASWP performs a series of row interchanges on the matrix A.
|
|
* One row interchange is initiated for each of rows K1 through K2 of A.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* N (input) INTEGER
|
|
* The number of columns of the matrix A.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the matrix of column dimension N to which the row
|
|
* interchanges will be applied.
|
|
* On exit, the permuted matrix.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A.
|
|
*
|
|
* K1 (input) INTEGER
|
|
* The first element of IPIV for which a row interchange will
|
|
* be done.
|
|
*
|
|
* K2 (input) INTEGER
|
|
* The last element of IPIV for which a row interchange will
|
|
* be done.
|
|
*
|
|
* IPIV (input) INTEGER array, dimension (M*abs(INCX))
|
|
* The vector of pivot indices. Only the elements in positions
|
|
* K1 through K2 of IPIV are accessed.
|
|
* IPIV(K) = L implies rows K and L are to be interchanged.
|
|
*
|
|
* INCX (input) INTEGER
|
|
* The increment between successive values of IPIV. If IPIV
|
|
* is negative, the pivots are applied in reverse order.
|
|
*
|
|
* Further Details
|
|
* ===============
|
|
*
|
|
* Modified by
|
|
* R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Local Scalars ..
|
|
INTEGER I, I1, I2, INC, IP, IX, IX0, J, K, N32
|
|
DOUBLE PRECISION TEMP
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Interchange row I with row IPIV(I) for each of rows K1 through K2.
|
|
*
|
|
IF( INCX.GT.0 ) THEN
|
|
IX0 = K1
|
|
I1 = K1
|
|
I2 = K2
|
|
INC = 1
|
|
ELSE IF( INCX.LT.0 ) THEN
|
|
IX0 = 1 + ( 1-K2 )*INCX
|
|
I1 = K2
|
|
I2 = K1
|
|
INC = -1
|
|
ELSE
|
|
RETURN
|
|
END IF
|
|
*
|
|
N32 = ( N / 32 )*32
|
|
IF( N32.NE.0 ) THEN
|
|
DO 30 J = 1, N32, 32
|
|
IX = IX0
|
|
DO 20 I = I1, I2, INC
|
|
IP = IPIV( IX )
|
|
IF( IP.NE.I ) THEN
|
|
DO 10 K = J, J + 31
|
|
TEMP = A( I, K )
|
|
A( I, K ) = A( IP, K )
|
|
A( IP, K ) = TEMP
|
|
10 CONTINUE
|
|
END IF
|
|
IX = IX + INCX
|
|
20 CONTINUE
|
|
30 CONTINUE
|
|
END IF
|
|
IF( N32.NE.N ) THEN
|
|
N32 = N32 + 1
|
|
IX = IX0
|
|
DO 50 I = I1, I2, INC
|
|
IP = IPIV( IX )
|
|
IF( IP.NE.I ) THEN
|
|
DO 40 K = N32, N
|
|
TEMP = A( I, K )
|
|
A( I, K ) = A( IP, K )
|
|
A( IP, K ) = TEMP
|
|
40 CONTINUE
|
|
END IF
|
|
IX = IX + INCX
|
|
50 CONTINUE
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of DLASWP
|
|
*
|
|
END
|
|
SUBROUTINE DPOTF2( UPLO, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* February 29, 1992
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
DOUBLE PRECISION A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DPOTF2 computes the Cholesky factorization of a real symmetric
|
|
* positive definite matrix A.
|
|
*
|
|
* The factorization has the form
|
|
* A = U' * U , if UPLO = 'U', or
|
|
* A = L * L', if UPLO = 'L',
|
|
* where U is an upper triangular matrix and L is lower triangular.
|
|
*
|
|
* This is the unblocked version of the algorithm, calling Level 2 BLAS.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* Specifies whether the upper or lower triangular part of the
|
|
* symmetric matrix A is stored.
|
|
* = 'U': Upper triangular
|
|
* = 'L': Lower triangular
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the symmetric matrix A. If UPLO = 'U', the leading
|
|
* n by n upper triangular part of A contains the upper
|
|
* triangular part of the matrix A, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading n by n lower triangular part of A contains the lower
|
|
* triangular part of the matrix A, and the strictly upper
|
|
* triangular part of A is not referenced.
|
|
*
|
|
* On exit, if INFO = 0, the factor U or L from the Cholesky
|
|
* factorization A = U'*U or A = L*L'.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -k, the k-th argument had an illegal value
|
|
* > 0: if INFO = k, the leading minor of order k is not
|
|
* positive definite, and the factorization could not be
|
|
* completed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE, ZERO
|
|
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL UPPER
|
|
INTEGER J
|
|
DOUBLE PRECISION AJJ
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
DOUBLE PRECISION DDOT
|
|
EXTERNAL LSAME, DDOT
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DGEMV, DSCAL, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, SQRT
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DPOTF2', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute the Cholesky factorization A = U'*U.
|
|
*
|
|
DO 10 J = 1, N
|
|
*
|
|
* Compute U(J,J) and test for non-positive-definiteness.
|
|
*
|
|
AJJ = A( J, J ) - DDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 )
|
|
IF( AJJ.LE.ZERO ) THEN
|
|
A( J, J ) = AJJ
|
|
GO TO 30
|
|
END IF
|
|
AJJ = SQRT( AJJ )
|
|
A( J, J ) = AJJ
|
|
*
|
|
* Compute elements J+1:N of row J.
|
|
*
|
|
IF( J.LT.N ) THEN
|
|
CALL DGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ),
|
|
$ LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
|
|
CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
|
|
END IF
|
|
10 CONTINUE
|
|
ELSE
|
|
*
|
|
* Compute the Cholesky factorization A = L*L'.
|
|
*
|
|
DO 20 J = 1, N
|
|
*
|
|
* Compute L(J,J) and test for non-positive-definiteness.
|
|
*
|
|
AJJ = A( J, J ) - DDOT( J-1, A( J, 1 ), LDA, A( J, 1 ),
|
|
$ LDA )
|
|
IF( AJJ.LE.ZERO ) THEN
|
|
A( J, J ) = AJJ
|
|
GO TO 30
|
|
END IF
|
|
AJJ = SQRT( AJJ )
|
|
A( J, J ) = AJJ
|
|
*
|
|
* Compute elements J+1:N of column J.
|
|
*
|
|
IF( J.LT.N ) THEN
|
|
CALL DGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ),
|
|
$ LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
|
|
CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
GO TO 40
|
|
*
|
|
30 CONTINUE
|
|
INFO = J
|
|
*
|
|
40 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of DPOTF2
|
|
*
|
|
END
|
|
SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* March 31, 1993
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
DOUBLE PRECISION A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DPOTRF computes the Cholesky factorization of a real symmetric
|
|
* positive definite matrix A.
|
|
*
|
|
* The factorization has the form
|
|
* A = U**T * U, if UPLO = 'U', or
|
|
* A = L * L**T, if UPLO = 'L',
|
|
* where U is an upper triangular matrix and L is lower triangular.
|
|
*
|
|
* This is the block version of the algorithm, calling Level 3 BLAS.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* = 'U': Upper triangle of A is stored;
|
|
* = 'L': Lower triangle of A is stored.
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the symmetric matrix A. If UPLO = 'U', the leading
|
|
* N-by-N upper triangular part of A contains the upper
|
|
* triangular part of the matrix A, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading N-by-N lower triangular part of A contains the lower
|
|
* triangular part of the matrix A, and the strictly upper
|
|
* triangular part of A is not referenced.
|
|
*
|
|
* On exit, if INFO = 0, the factor U or L from the Cholesky
|
|
* factorization A = U**T*U or A = L*L**T.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, the leading minor of order i is not
|
|
* positive definite, and the factorization could not be
|
|
* completed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE
|
|
PARAMETER ( ONE = 1.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL UPPER
|
|
INTEGER J, JB, NB
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
INTEGER ILAENV
|
|
EXTERNAL LSAME, ILAENV
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DPOTRF', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Determine the block size for this environment.
|
|
*
|
|
NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 )
|
|
IF( NB.LE.1 .OR. NB.GE.N ) THEN
|
|
*
|
|
* Use unblocked code.
|
|
*
|
|
CALL DPOTF2( UPLO, N, A, LDA, INFO )
|
|
ELSE
|
|
*
|
|
* Use blocked code.
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute the Cholesky factorization A = U'*U.
|
|
*
|
|
DO 10 J = 1, N, NB
|
|
*
|
|
* Update and factorize the current diagonal block and test
|
|
* for non-positive-definiteness.
|
|
*
|
|
JB = MIN( NB, N-J+1 )
|
|
CALL DSYRK( 'Upper', 'Transpose', JB, J-1, -ONE,
|
|
$ A( 1, J ), LDA, ONE, A( J, J ), LDA )
|
|
CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
|
|
IF( INFO.NE.0 )
|
|
$ GO TO 30
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Compute the current block row.
|
|
*
|
|
CALL DGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1,
|
|
$ J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ),
|
|
$ LDA, ONE, A( J, J+JB ), LDA )
|
|
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit',
|
|
$ JB, N-J-JB+1, ONE, A( J, J ), LDA,
|
|
$ A( J, J+JB ), LDA )
|
|
END IF
|
|
10 CONTINUE
|
|
*
|
|
ELSE
|
|
*
|
|
* Compute the Cholesky factorization A = L*L'.
|
|
*
|
|
DO 20 J = 1, N, NB
|
|
*
|
|
* Update and factorize the current diagonal block and test
|
|
* for non-positive-definiteness.
|
|
*
|
|
JB = MIN( NB, N-J+1 )
|
|
CALL DSYRK( 'Lower', 'No transpose', JB, J-1, -ONE,
|
|
$ A( J, 1 ), LDA, ONE, A( J, J ), LDA )
|
|
CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
|
|
IF( INFO.NE.0 )
|
|
$ GO TO 30
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Compute the current block column.
|
|
*
|
|
CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
|
|
$ J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ),
|
|
$ LDA, ONE, A( J+JB, J ), LDA )
|
|
CALL DTRSM( 'Right', 'Lower', 'Transpose', 'Non-unit',
|
|
$ N-J-JB+1, JB, ONE, A( J, J ), LDA,
|
|
$ A( J+JB, J ), LDA )
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
END IF
|
|
GO TO 40
|
|
*
|
|
30 CONTINUE
|
|
INFO = INFO + J - 1
|
|
*
|
|
40 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of DPOTRF
|
|
*
|
|
END
|
|
SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* March 31, 1993
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER UPLO
|
|
INTEGER INFO, LDA, LDB, N, NRHS
|
|
* ..
|
|
* .. Array Arguments ..
|
|
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DPOTRS solves a system of linear equations A*X = B with a symmetric
|
|
* positive definite matrix A using the Cholesky factorization
|
|
* A = U**T*U or A = L*L**T computed by DPOTRF.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* = 'U': Upper triangle of A is stored;
|
|
* = 'L': Lower triangle of A is stored.
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* NRHS (input) INTEGER
|
|
* The number of right hand sides, i.e., the number of columns
|
|
* of the matrix B. NRHS >= 0.
|
|
*
|
|
* A (input) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* The triangular factor U or L from the Cholesky factorization
|
|
* A = U**T*U or A = L*L**T, as computed by DPOTRF.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
|
|
* On entry, the right hand side matrix B.
|
|
* On exit, the solution matrix X.
|
|
*
|
|
* LDB (input) INTEGER
|
|
* The leading dimension of the array B. LDB >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE
|
|
PARAMETER ( ONE = 1.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL UPPER
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
EXTERNAL LSAME
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DTRSM, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( NRHS.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
|
|
INFO = -7
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DPOTRS', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 .OR. NRHS.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Solve A*X = B where A = U'*U.
|
|
*
|
|
* Solve U'*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS,
|
|
$ ONE, A, LDA, B, LDB )
|
|
*
|
|
* Solve U*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
|
|
$ NRHS, ONE, A, LDA, B, LDB )
|
|
ELSE
|
|
*
|
|
* Solve A*X = B where A = L*L'.
|
|
*
|
|
* Solve L*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N,
|
|
$ NRHS, ONE, A, LDA, B, LDB )
|
|
*
|
|
* Solve L'*X = B, overwriting B with X.
|
|
*
|
|
CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', N, NRHS,
|
|
$ ONE, A, LDA, B, LDB )
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of DPOTRS
|
|
*
|
|
END
|
|
|
|
SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* February 29, 1992
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER DIAG, UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
DOUBLE PRECISION A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DTRTI2 computes the inverse of a real upper or lower triangular
|
|
* matrix.
|
|
*
|
|
* This is the Level 2 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* Specifies whether the matrix A is upper or lower triangular.
|
|
* = 'U': Upper triangular
|
|
* = 'L': Lower triangular
|
|
*
|
|
* DIAG (input) CHARACTER*1
|
|
* Specifies whether or not the matrix A is unit triangular.
|
|
* = 'N': Non-unit triangular
|
|
* = 'U': Unit triangular
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the triangular matrix A. If UPLO = 'U', the
|
|
* leading n by n upper triangular part of the array A contains
|
|
* the upper triangular matrix, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading n by n lower triangular part of the array A contains
|
|
* the lower triangular matrix, and the strictly upper
|
|
* triangular part of A is not referenced. If DIAG = 'U', the
|
|
* diagonal elements of A are also not referenced and are
|
|
* assumed to be 1.
|
|
*
|
|
* On exit, the (triangular) inverse of the original matrix, in
|
|
* the same storage format.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -k, the k-th argument had an illegal value
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE
|
|
PARAMETER ( ONE = 1.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL NOUNIT, UPPER
|
|
INTEGER J
|
|
DOUBLE PRECISION AJJ
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
EXTERNAL LSAME
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DSCAL, DTRMV, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
NOUNIT = LSAME( DIAG, 'N' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
|
|
INFO = -2
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DTRTI2', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute inverse of upper triangular matrix.
|
|
*
|
|
DO 10 J = 1, N
|
|
IF( NOUNIT ) THEN
|
|
A( J, J ) = ONE / A( J, J )
|
|
AJJ = -A( J, J )
|
|
ELSE
|
|
AJJ = -ONE
|
|
END IF
|
|
*
|
|
* Compute elements 1:j-1 of j-th column.
|
|
*
|
|
CALL DTRMV( 'Upper', 'No transpose', DIAG, J-1, A, LDA,
|
|
$ A( 1, J ), 1 )
|
|
CALL DSCAL( J-1, AJJ, A( 1, J ), 1 )
|
|
10 CONTINUE
|
|
ELSE
|
|
*
|
|
* Compute inverse of lower triangular matrix.
|
|
*
|
|
DO 20 J = N, 1, -1
|
|
IF( NOUNIT ) THEN
|
|
A( J, J ) = ONE / A( J, J )
|
|
AJJ = -A( J, J )
|
|
ELSE
|
|
AJJ = -ONE
|
|
END IF
|
|
IF( J.LT.N ) THEN
|
|
*
|
|
* Compute elements j+1:n of j-th column.
|
|
*
|
|
CALL DTRMV( 'Lower', 'No transpose', DIAG, N-J,
|
|
$ A( J+1, J+1 ), LDA, A( J+1, J ), 1 )
|
|
CALL DSCAL( N-J, AJJ, A( J+1, J ), 1 )
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of DTRTI2
|
|
*
|
|
END
|
|
SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* March 31, 1993
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER DIAG, UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
DOUBLE PRECISION A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* DTRTRI computes the inverse of a real upper or lower triangular
|
|
* matrix A.
|
|
*
|
|
* This is the Level 3 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* = 'U': A is upper triangular;
|
|
* = 'L': A is lower triangular.
|
|
*
|
|
* DIAG (input) CHARACTER*1
|
|
* = 'N': A is non-unit triangular;
|
|
* = 'U': A is unit triangular.
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
|
* On entry, the triangular matrix A. If UPLO = 'U', the
|
|
* leading N-by-N upper triangular part of the array A contains
|
|
* the upper triangular matrix, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading N-by-N lower triangular part of the array A contains
|
|
* the lower triangular matrix, and the strictly upper
|
|
* triangular part of A is not referenced. If DIAG = 'U', the
|
|
* diagonal elements of A are also not referenced and are
|
|
* assumed to be 1.
|
|
* On exit, the (triangular) inverse of the original matrix, in
|
|
* the same storage format.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, A(i,i) is exactly zero. The triangular
|
|
* matrix is singular and its inverse can not be computed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE, ZERO
|
|
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL NOUNIT, UPPER
|
|
INTEGER J, JB, NB, NN
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
INTEGER ILAENV
|
|
EXTERNAL LSAME, ILAENV
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL DTRMM, DTRSM, DTRTI2, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
NOUNIT = LSAME( DIAG, 'N' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
|
|
INFO = -2
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DTRTRI', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Check for singularity if non-unit.
|
|
*
|
|
IF( NOUNIT ) THEN
|
|
DO 10 INFO = 1, N
|
|
IF( A( INFO, INFO ).EQ.ZERO )
|
|
$ RETURN
|
|
10 CONTINUE
|
|
INFO = 0
|
|
END IF
|
|
*
|
|
* Determine the block size for this environment.
|
|
*
|
|
NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 )
|
|
IF( NB.LE.1 .OR. NB.GE.N ) THEN
|
|
*
|
|
* Use unblocked code
|
|
*
|
|
CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
|
|
ELSE
|
|
*
|
|
* Use blocked code
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute inverse of upper triangular matrix
|
|
*
|
|
DO 20 J = 1, N, NB
|
|
JB = MIN( NB, N-J+1 )
|
|
*
|
|
* Compute rows 1:j-1 of current block column
|
|
*
|
|
CALL DTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1,
|
|
$ JB, ONE, A, LDA, A( 1, J ), LDA )
|
|
CALL DTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1,
|
|
$ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA )
|
|
*
|
|
* Compute inverse of current diagonal block
|
|
*
|
|
CALL DTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO )
|
|
20 CONTINUE
|
|
ELSE
|
|
*
|
|
* Compute inverse of lower triangular matrix
|
|
*
|
|
NN = ( ( N-1 ) / NB )*NB + 1
|
|
DO 30 J = NN, 1, -NB
|
|
JB = MIN( NB, N-J+1 )
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Compute rows j+jb:n of current block column
|
|
*
|
|
CALL DTRMM( 'Left', 'Lower', 'No transpose', DIAG,
|
|
$ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA,
|
|
$ A( J+JB, J ), LDA )
|
|
CALL DTRSM( 'Right', 'Lower', 'No transpose', DIAG,
|
|
$ N-J-JB+1, JB, -ONE, A( J, J ), LDA,
|
|
$ A( J+JB, J ), LDA )
|
|
END IF
|
|
*
|
|
* Compute inverse of current diagonal block
|
|
*
|
|
CALL DTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO )
|
|
30 CONTINUE
|
|
END IF
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of DTRTRI
|
|
*
|
|
END
|
|
LOGICAL FUNCTION LSAME( CA, CB )
|
|
*
|
|
* -- LAPACK auxiliary routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER CA, CB
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* LSAME returns .TRUE. if CA is the same letter as CB regardless of
|
|
* case.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* CA (input) CHARACTER*1
|
|
* CB (input) CHARACTER*1
|
|
* CA and CB specify the single characters to be compared.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC ICHAR
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER INTA, INTB, ZCODE
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test if the characters are equal
|
|
*
|
|
LSAME = CA.EQ.CB
|
|
IF( LSAME )
|
|
$ RETURN
|
|
*
|
|
* Now test for equivalence if both characters are alphabetic.
|
|
*
|
|
ZCODE = ICHAR( 'Z' )
|
|
*
|
|
* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
|
|
* machines, on which ICHAR returns a value with bit 8 set.
|
|
* ICHAR('A') on Prime machines returns 193 which is the same as
|
|
* ICHAR('A') on an EBCDIC machine.
|
|
*
|
|
INTA = ICHAR( CA )
|
|
INTB = ICHAR( CB )
|
|
*
|
|
IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN
|
|
*
|
|
* ASCII is assumed - ZCODE is the ASCII code of either lower or
|
|
* upper case 'Z'.
|
|
*
|
|
IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32
|
|
IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32
|
|
*
|
|
ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN
|
|
*
|
|
* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
|
|
* upper case 'Z'.
|
|
*
|
|
IF( INTA.GE.129 .AND. INTA.LE.137 .OR.
|
|
$ INTA.GE.145 .AND. INTA.LE.153 .OR.
|
|
$ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64
|
|
IF( INTB.GE.129 .AND. INTB.LE.137 .OR.
|
|
$ INTB.GE.145 .AND. INTB.LE.153 .OR.
|
|
$ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64
|
|
*
|
|
ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN
|
|
*
|
|
* ASCII is assumed, on Prime machines - ZCODE is the ASCII code
|
|
* plus 128 of either lower or upper case 'Z'.
|
|
*
|
|
IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32
|
|
IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32
|
|
END IF
|
|
LSAME = INTA.EQ.INTB
|
|
*
|
|
* RETURN
|
|
*
|
|
* End of LSAME
|
|
*
|
|
END
|
|
SUBROUTINE XERBLA( SRNAME, INFO )
|
|
*
|
|
* -- LAPACK auxiliary routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER*6 SRNAME
|
|
INTEGER INFO
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* XERBLA is an error handler for the LAPACK routines.
|
|
* It is called by an LAPACK routine if an input parameter has an
|
|
* invalid value. A message is printed and execution stops.
|
|
*
|
|
* Installers may consider modifying the STOP statement in order to
|
|
* call system-specific exception-handling facilities.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* SRNAME (input) CHARACTER*6
|
|
* The name of the routine which called XERBLA.
|
|
*
|
|
* INFO (input) INTEGER
|
|
* The position of the invalid parameter in the parameter list
|
|
* of the calling routine.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Executable Statements ..
|
|
*
|
|
WRITE( *, FMT = 9999 )SRNAME, INFO
|
|
*
|
|
STOP
|
|
*
|
|
9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ',
|
|
$ 'an illegal value' )
|
|
*
|
|
* End of XERBLA
|
|
*
|
|
END
|
|
SUBROUTINE ZDROT( N, CX, INCX, CY, INCY, C, S )
|
|
*
|
|
* applies a plane rotation, where the cos and sin (c and s) are real
|
|
* and the vectors cx and cy are complex.
|
|
* jack dongarra, linpack, 3/11/78.
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INCX, INCY, N
|
|
DOUBLE PRECISION C, S
|
|
* ..
|
|
* .. Array Arguments ..
|
|
COMPLEX*16 CX( * ), CY( * )
|
|
*
|
|
* =====================================================================
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER I, IX, IY
|
|
COMPLEX*16 CTEMP
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
IF( N.LE.0 )
|
|
$ RETURN
|
|
IF( INCX.EQ.1 .AND. INCY.EQ.1 )
|
|
$ GO TO 20
|
|
*
|
|
* code for unequal increments or equal increments not equal
|
|
* to 1
|
|
*
|
|
IX = 1
|
|
IY = 1
|
|
IF( INCX.LT.0 )
|
|
$ IX = ( -N+1 )*INCX + 1
|
|
IF( INCY.LT.0 )
|
|
$ IY = ( -N+1 )*INCY + 1
|
|
DO 10 I = 1, N
|
|
CTEMP = C*CX( IX ) + S*CY( IY )
|
|
CY( IY ) = C*CY( IY ) - S*CX( IX )
|
|
CX( IX ) = CTEMP
|
|
IX = IX + INCX
|
|
IY = IY + INCY
|
|
10 CONTINUE
|
|
RETURN
|
|
*
|
|
* code for both increments equal to 1
|
|
*
|
|
20 CONTINUE
|
|
DO 30 I = 1, N
|
|
CTEMP = C*CX( I ) + S*CY( I )
|
|
CY( I ) = C*CY( I ) - S*CX( I )
|
|
CX( I ) = CTEMP
|
|
30 CONTINUE
|
|
RETURN
|
|
END
|
|
SUBROUTINE ZGETF2( M, N, A, LDA, IPIV, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, M, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
COMPLEX*16 A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZGETF2 computes an LU factorization of a general m-by-n matrix A
|
|
* using partial pivoting with row interchanges.
|
|
*
|
|
* The factorization has the form
|
|
* A = P * L * U
|
|
* where P is a permutation matrix, L is lower triangular with unit
|
|
* diagonal elements (lower trapezoidal if m > n), and U is upper
|
|
* triangular (upper trapezoidal if m < n).
|
|
*
|
|
* This is the right-looking Level 2 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* M (input) INTEGER
|
|
* The number of rows of the matrix A. M >= 0.
|
|
*
|
|
* N (input) INTEGER
|
|
* The number of columns of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
|
|
* On entry, the m by n matrix to be factored.
|
|
* On exit, the factors L and U from the factorization
|
|
* A = P*L*U; the unit diagonal elements of L are not stored.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,M).
|
|
*
|
|
* IPIV (output) INTEGER array, dimension (min(M,N))
|
|
* The pivot indices; for 1 <= i <= min(M,N), row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -k, the k-th argument had an illegal value
|
|
* > 0: if INFO = k, U(k,k) is exactly zero. The factorization
|
|
* has been completed, but the factor U is exactly
|
|
* singular, and division by zero will occur if it is used
|
|
* to solve a system of equations.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
COMPLEX*16 ONE, ZERO
|
|
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
|
|
$ ZERO = ( 0.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER J, JP
|
|
* ..
|
|
* .. External Functions ..
|
|
INTEGER IZAMAX
|
|
EXTERNAL IZAMAX
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZGERU, ZSCAL, ZSWAP
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
IF( M.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZGETF2', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( M.EQ.0 .OR. N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
DO 10 J = 1, MIN( M, N )
|
|
*
|
|
* Find pivot and test for singularity.
|
|
*
|
|
JP = J - 1 + IZAMAX( M-J+1, A( J, J ), 1 )
|
|
IPIV( J ) = JP
|
|
IF( A( JP, J ).NE.ZERO ) THEN
|
|
*
|
|
* Apply the interchange to columns 1:N.
|
|
*
|
|
IF( JP.NE.J )
|
|
$ CALL ZSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
|
|
*
|
|
* Compute elements J+1:M of J-th column.
|
|
*
|
|
IF( J.LT.M )
|
|
$ CALL ZSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 )
|
|
*
|
|
ELSE IF( INFO.EQ.0 ) THEN
|
|
*
|
|
INFO = J
|
|
END IF
|
|
*
|
|
IF( J.LT.MIN( M, N ) ) THEN
|
|
*
|
|
* Update trailing submatrix.
|
|
*
|
|
CALL ZGERU( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ),
|
|
$ LDA, A( J+1, J+1 ), LDA )
|
|
END IF
|
|
10 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of ZGETF2
|
|
*
|
|
END
|
|
SUBROUTINE ZGETRF( M, N, A, LDA, IPIV, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, M, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
COMPLEX*16 A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZGETRF computes an LU factorization of a general M-by-N matrix A
|
|
* using partial pivoting with row interchanges.
|
|
*
|
|
* The factorization has the form
|
|
* A = P * L * U
|
|
* where P is a permutation matrix, L is lower triangular with unit
|
|
* diagonal elements (lower trapezoidal if m > n), and U is upper
|
|
* triangular (upper trapezoidal if m < n).
|
|
*
|
|
* This is the right-looking Level 3 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* M (input) INTEGER
|
|
* The number of rows of the matrix A. M >= 0.
|
|
*
|
|
* N (input) INTEGER
|
|
* The number of columns of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
|
|
* On entry, the M-by-N matrix to be factored.
|
|
* On exit, the factors L and U from the factorization
|
|
* A = P*L*U; the unit diagonal elements of L are not stored.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,M).
|
|
*
|
|
* IPIV (output) INTEGER array, dimension (min(M,N))
|
|
* The pivot indices; for 1 <= i <= min(M,N), row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, U(i,i) is exactly zero. The factorization
|
|
* has been completed, but the factor U is exactly
|
|
* singular, and division by zero will occur if it is used
|
|
* to solve a system of equations.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
COMPLEX*16 ONE
|
|
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
INTEGER I, IINFO, J, JB, NB
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZGEMM, ZGETF2, ZLASWP, ZTRSM
|
|
* ..
|
|
* .. External Functions ..
|
|
INTEGER ILAENV
|
|
EXTERNAL ILAENV
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
IF( M.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZGETRF', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( M.EQ.0 .OR. N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Determine the block size for this environment.
|
|
*
|
|
NB = ILAENV( 1, 'ZGETRF', ' ', M, N, -1, -1 )
|
|
IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
|
|
*
|
|
* Use unblocked code.
|
|
*
|
|
CALL ZGETF2( M, N, A, LDA, IPIV, INFO )
|
|
ELSE
|
|
*
|
|
* Use blocked code.
|
|
*
|
|
DO 20 J = 1, MIN( M, N ), NB
|
|
JB = MIN( MIN( M, N )-J+1, NB )
|
|
*
|
|
* Factor diagonal and subdiagonal blocks and test for exact
|
|
* singularity.
|
|
*
|
|
CALL ZGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO )
|
|
*
|
|
* Adjust INFO and the pivot indices.
|
|
*
|
|
IF( INFO.EQ.0 .AND. IINFO.GT.0 )
|
|
$ INFO = IINFO + J - 1
|
|
DO 10 I = J, MIN( M, J+JB-1 )
|
|
IPIV( I ) = J - 1 + IPIV( I )
|
|
10 CONTINUE
|
|
*
|
|
* Apply interchanges to columns 1:J-1.
|
|
*
|
|
CALL ZLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
|
|
*
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Apply interchanges to columns J+JB:N.
|
|
*
|
|
CALL ZLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1,
|
|
$ IPIV, 1 )
|
|
*
|
|
* Compute block row of U.
|
|
*
|
|
CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
|
|
$ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
|
|
$ LDA )
|
|
IF( J+JB.LE.M ) THEN
|
|
*
|
|
* Update trailing submatrix.
|
|
*
|
|
CALL ZGEMM( 'No transpose', 'No transpose', M-J-JB+1,
|
|
$ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
|
|
$ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
|
|
$ LDA )
|
|
END IF
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
RETURN
|
|
*
|
|
* End of ZGETRF
|
|
*
|
|
END
|
|
|
|
SUBROUTINE ZGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* June 30, 1999
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, LWORK, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
COMPLEX*16 A( LDA, * ), WORK( * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZGETRI computes the inverse of a matrix using the LU factorization
|
|
* computed by ZGETRF.
|
|
*
|
|
* This method inverts U and then computes inv(A) by solving the system
|
|
* inv(A)*L = inv(U) for inv(A).
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
|
|
* On entry, the factors L and U from the factorization
|
|
* A = P*L*U as computed by ZGETRF.
|
|
* On exit, if INFO = 0, the inverse of the original matrix A.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* IPIV (input) INTEGER array, dimension (N)
|
|
* The pivot indices from ZGETRF; for 1<=i<=N, row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* WORK (workspace/output) COMPLEX*16 array, dimension (LWORK)
|
|
* On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
|
|
*
|
|
* LWORK (input) INTEGER
|
|
* The dimension of the array WORK. LWORK >= max(1,N).
|
|
* For optimal performance LWORK >= N*NB, where NB is
|
|
* the optimal blocksize returned by ILAENV.
|
|
*
|
|
* If LWORK = -1, then a workspace query is assumed; the routine
|
|
* only calculates the optimal size of the WORK array, returns
|
|
* this value as the first entry of the WORK array, and no error
|
|
* message related to LWORK is issued by XERBLA.
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
|
|
* singular and its inverse could not be computed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
COMPLEX*16 ZERO, ONE
|
|
PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
|
|
$ ONE = ( 1.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL LQUERY
|
|
INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
|
|
$ NBMIN, NN
|
|
* ..
|
|
* .. External Functions ..
|
|
INTEGER ILAENV
|
|
EXTERNAL ILAENV
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZGEMM, ZGEMV, ZSWAP, ZTRSM, ZTRTRI
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
NB = ILAENV( 1, 'ZGETRI', ' ', N, -1, -1, -1 )
|
|
LWKOPT = N*NB
|
|
WORK( 1 ) = LWKOPT
|
|
LQUERY = ( LWORK.EQ.-1 )
|
|
IF( N.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -3
|
|
ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
|
|
INFO = -6
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZGETRI', -INFO )
|
|
RETURN
|
|
ELSE IF( LQUERY ) THEN
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Form inv(U). If INFO > 0 from ZTRTRI, then U is singular,
|
|
* and the inverse is not computed.
|
|
*
|
|
CALL ZTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
|
|
IF( INFO.GT.0 )
|
|
$ RETURN
|
|
*
|
|
NBMIN = 2
|
|
LDWORK = N
|
|
IF( NB.GT.1 .AND. NB.LT.N ) THEN
|
|
IWS = MAX( LDWORK*NB, 1 )
|
|
IF( LWORK.LT.IWS ) THEN
|
|
NB = LWORK / LDWORK
|
|
NBMIN = MAX( 2, ILAENV( 2, 'ZGETRI', ' ', N, -1, -1, -1 ) )
|
|
END IF
|
|
ELSE
|
|
IWS = N
|
|
END IF
|
|
*
|
|
* Solve the equation inv(A)*L = inv(U) for inv(A).
|
|
*
|
|
IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
|
|
*
|
|
* Use unblocked code.
|
|
*
|
|
DO 20 J = N, 1, -1
|
|
*
|
|
* Copy current column of L to WORK and replace with zeros.
|
|
*
|
|
DO 10 I = J + 1, N
|
|
WORK( I ) = A( I, J )
|
|
A( I, J ) = ZERO
|
|
10 CONTINUE
|
|
*
|
|
* Compute current column of inv(A).
|
|
*
|
|
IF( J.LT.N )
|
|
$ CALL ZGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
|
|
$ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
|
|
20 CONTINUE
|
|
ELSE
|
|
*
|
|
* Use blocked code.
|
|
*
|
|
NN = ( ( N-1 ) / NB )*NB + 1
|
|
DO 50 J = NN, 1, -NB
|
|
JB = MIN( NB, N-J+1 )
|
|
*
|
|
* Copy current block column of L to WORK and replace with
|
|
* zeros.
|
|
*
|
|
DO 40 JJ = J, J + JB - 1
|
|
DO 30 I = JJ + 1, N
|
|
WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
|
|
A( I, JJ ) = ZERO
|
|
30 CONTINUE
|
|
40 CONTINUE
|
|
*
|
|
* Compute current block column of inv(A).
|
|
*
|
|
IF( J+JB.LE.N )
|
|
$ CALL ZGEMM( 'No transpose', 'No transpose', N, JB,
|
|
$ N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
|
|
$ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
|
|
CALL ZTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
|
|
$ ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
|
|
50 CONTINUE
|
|
END IF
|
|
*
|
|
* Apply column interchanges.
|
|
*
|
|
DO 60 J = N - 1, 1, -1
|
|
JP = IPIV( J )
|
|
IF( JP.NE.J )
|
|
$ CALL ZSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
|
|
60 CONTINUE
|
|
*
|
|
WORK( 1 ) = IWS
|
|
RETURN
|
|
*
|
|
* End of ZGETRI
|
|
*
|
|
END
|
|
SUBROUTINE ZGETRS( TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER TRANS
|
|
INTEGER INFO, LDA, LDB, N, NRHS
|
|
* ..
|
|
* .. Array Arguments ..
|
|
INTEGER IPIV( * )
|
|
COMPLEX*16 A( LDA, * ), B( LDB, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZGETRS solves a system of linear equations
|
|
* A * X = B, A**T * X = B, or A**H * X = B
|
|
* with a general N-by-N matrix A using the LU factorization computed
|
|
* by ZGETRF.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* TRANS (input) CHARACTER*1
|
|
* Specifies the form of the system of equations:
|
|
* = 'N': A * X = B (No transpose)
|
|
* = 'T': A**T * X = B (Transpose)
|
|
* = 'C': A**H * X = B (Conjugate transpose)
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* NRHS (input) INTEGER
|
|
* The number of right hand sides, i.e., the number of columns
|
|
* of the matrix B. NRHS >= 0.
|
|
*
|
|
* A (input) COMPLEX*16 array, dimension (LDA,N)
|
|
* The factors L and U from the factorization A = P*L*U
|
|
* as computed by ZGETRF.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* IPIV (input) INTEGER array, dimension (N)
|
|
* The pivot indices from ZGETRF; for 1<=i<=N, row i of the
|
|
* matrix was interchanged with row IPIV(i).
|
|
*
|
|
* B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
|
|
* On entry, the right hand side matrix B.
|
|
* On exit, the solution matrix X.
|
|
*
|
|
* LDB (input) INTEGER
|
|
* The leading dimension of the array B. LDB >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
COMPLEX*16 ONE
|
|
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL NOTRAN
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
EXTERNAL LSAME
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZLASWP, ZTRSM
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
NOTRAN = LSAME( TRANS, 'N' )
|
|
IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
|
|
$ LSAME( TRANS, 'C' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( NRHS.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
|
|
INFO = -8
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZGETRS', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 .OR. NRHS.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
IF( NOTRAN ) THEN
|
|
*
|
|
* Solve A * X = B.
|
|
*
|
|
* Apply row interchanges to the right hand sides.
|
|
*
|
|
CALL ZLASWP( NRHS, B, LDB, 1, N, IPIV, 1 )
|
|
*
|
|
* Solve L*X = B, overwriting B with X.
|
|
*
|
|
CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', N, NRHS,
|
|
$ ONE, A, LDA, B, LDB )
|
|
*
|
|
* Solve U*X = B, overwriting B with X.
|
|
*
|
|
CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
|
|
$ NRHS, ONE, A, LDA, B, LDB )
|
|
ELSE
|
|
*
|
|
* Solve A**T * X = B or A**H * X = B.
|
|
*
|
|
* Solve U'*X = B, overwriting B with X.
|
|
*
|
|
CALL ZTRSM( 'Left', 'Upper', TRANS, 'Non-unit', N, NRHS, ONE,
|
|
$ A, LDA, B, LDB )
|
|
*
|
|
* Solve L'*X = B, overwriting B with X.
|
|
*
|
|
CALL ZTRSM( 'Left', 'Lower', TRANS, 'Unit', N, NRHS, ONE, A,
|
|
$ LDA, B, LDB )
|
|
*
|
|
* Apply row interchanges to the solution vectors.
|
|
*
|
|
CALL ZLASWP( NRHS, B, LDB, 1, N, IPIV, -1 )
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of ZGETRS
|
|
*
|
|
END
|
|
|
|
|
|
SUBROUTINE ZPOTF2( UPLO, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
COMPLEX*16 A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZPOTF2 computes the Cholesky factorization of a complex Hermitian
|
|
* positive definite matrix A.
|
|
*
|
|
* The factorization has the form
|
|
* A = U' * U , if UPLO = 'U', or
|
|
* A = L * L', if UPLO = 'L',
|
|
* where U is an upper triangular matrix and L is lower triangular.
|
|
*
|
|
* This is the unblocked version of the algorithm, calling Level 2 BLAS.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* Specifies whether the upper or lower triangular part of the
|
|
* Hermitian matrix A is stored.
|
|
* = 'U': Upper triangular
|
|
* = 'L': Lower triangular
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
|
|
* On entry, the Hermitian matrix A. If UPLO = 'U', the leading
|
|
* n by n upper triangular part of A contains the upper
|
|
* triangular part of the matrix A, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading n by n lower triangular part of A contains the lower
|
|
* triangular part of the matrix A, and the strictly upper
|
|
* triangular part of A is not referenced.
|
|
*
|
|
* On exit, if INFO = 0, the factor U or L from the Cholesky
|
|
* factorization A = U'*U or A = L*L'.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -k, the k-th argument had an illegal value
|
|
* > 0: if INFO = k, the leading minor of order k is not
|
|
* positive definite, and the factorization could not be
|
|
* completed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE, ZERO
|
|
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
|
|
COMPLEX*16 CONE
|
|
PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL UPPER
|
|
INTEGER J
|
|
DOUBLE PRECISION AJJ
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
COMPLEX*16 ZDOTC
|
|
EXTERNAL LSAME, ZDOTC
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZDSCAL, ZGEMV, ZLACGV
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC DBLE, MAX, SQRT
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZPOTF2', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute the Cholesky factorization A = U'*U.
|
|
*
|
|
DO 10 J = 1, N
|
|
*
|
|
* Compute U(J,J) and test for non-positive-definiteness.
|
|
*
|
|
AJJ = DBLE( A( J, J ) ) - ZDOTC( J-1, A( 1, J ), 1,
|
|
$ A( 1, J ), 1 )
|
|
IF( AJJ.LE.ZERO ) THEN
|
|
A( J, J ) = AJJ
|
|
GO TO 30
|
|
END IF
|
|
AJJ = SQRT( AJJ )
|
|
A( J, J ) = AJJ
|
|
*
|
|
* Compute elements J+1:N of row J.
|
|
*
|
|
IF( J.LT.N ) THEN
|
|
CALL ZLACGV( J-1, A( 1, J ), 1 )
|
|
CALL ZGEMV( 'Transpose', J-1, N-J, -CONE, A( 1, J+1 ),
|
|
$ LDA, A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
|
|
CALL ZLACGV( J-1, A( 1, J ), 1 )
|
|
CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
|
|
END IF
|
|
10 CONTINUE
|
|
ELSE
|
|
*
|
|
* Compute the Cholesky factorization A = L*L'.
|
|
*
|
|
DO 20 J = 1, N
|
|
*
|
|
* Compute L(J,J) and test for non-positive-definiteness.
|
|
*
|
|
AJJ = DBLE( A( J, J ) ) - ZDOTC( J-1, A( J, 1 ), LDA,
|
|
$ A( J, 1 ), LDA )
|
|
IF( AJJ.LE.ZERO ) THEN
|
|
A( J, J ) = AJJ
|
|
GO TO 30
|
|
END IF
|
|
AJJ = SQRT( AJJ )
|
|
A( J, J ) = AJJ
|
|
*
|
|
* Compute elements J+1:N of column J.
|
|
*
|
|
IF( J.LT.N ) THEN
|
|
CALL ZLACGV( J-1, A( J, 1 ), LDA )
|
|
CALL ZGEMV( 'No transpose', N-J, J-1, -CONE, A( J+1, 1 ),
|
|
$ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
|
|
CALL ZLACGV( J-1, A( J, 1 ), LDA )
|
|
CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
GO TO 40
|
|
*
|
|
30 CONTINUE
|
|
INFO = J
|
|
*
|
|
40 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of ZPOTF2
|
|
*
|
|
END
|
|
SUBROUTINE ZPOTRF( UPLO, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
COMPLEX*16 A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZPOTRF computes the Cholesky factorization of a complex Hermitian
|
|
* positive definite matrix A.
|
|
*
|
|
* The factorization has the form
|
|
* A = U**H * U, if UPLO = 'U', or
|
|
* A = L * L**H, if UPLO = 'L',
|
|
* where U is an upper triangular matrix and L is lower triangular.
|
|
*
|
|
* This is the block version of the algorithm, calling Level 3 BLAS.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* = 'U': Upper triangle of A is stored;
|
|
* = 'L': Lower triangle of A is stored.
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
|
|
* On entry, the Hermitian matrix A. If UPLO = 'U', the leading
|
|
* N-by-N upper triangular part of A contains the upper
|
|
* triangular part of the matrix A, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading N-by-N lower triangular part of A contains the lower
|
|
* triangular part of the matrix A, and the strictly upper
|
|
* triangular part of A is not referenced.
|
|
*
|
|
* On exit, if INFO = 0, the factor U or L from the Cholesky
|
|
* factorization A = U**H*U or A = L*L**H.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, the leading minor of order i is not
|
|
* positive definite, and the factorization could not be
|
|
* completed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ONE
|
|
COMPLEX*16 CONE
|
|
PARAMETER ( ONE = 1.0D+0, CONE = ( 1.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL UPPER
|
|
INTEGER J, JB, NB
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
INTEGER ILAENV
|
|
EXTERNAL LSAME, ILAENV
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZGEMM, ZHERK, ZPOTF2, ZTRSM
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -2
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -4
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZPOTRF', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Determine the block size for this environment.
|
|
*
|
|
NB = ILAENV( 1, 'ZPOTRF', UPLO, N, -1, -1, -1 )
|
|
IF( NB.LE.1 .OR. NB.GE.N ) THEN
|
|
*
|
|
* Use unblocked code.
|
|
*
|
|
CALL ZPOTF2( UPLO, N, A, LDA, INFO )
|
|
ELSE
|
|
*
|
|
* Use blocked code.
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute the Cholesky factorization A = U'*U.
|
|
*
|
|
DO 10 J = 1, N, NB
|
|
*
|
|
* Update and factorize the current diagonal block and test
|
|
* for non-positive-definiteness.
|
|
*
|
|
JB = MIN( NB, N-J+1 )
|
|
CALL ZHERK( 'Upper', 'Conjugate transpose', JB, J-1,
|
|
$ -ONE, A( 1, J ), LDA, ONE, A( J, J ), LDA )
|
|
CALL ZPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
|
|
IF( INFO.NE.0 )
|
|
$ GO TO 30
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Compute the current block row.
|
|
*
|
|
CALL ZGEMM( 'Conjugate transpose', 'No transpose', JB,
|
|
$ N-J-JB+1, J-1, -CONE, A( 1, J ), LDA,
|
|
$ A( 1, J+JB ), LDA, CONE, A( J, J+JB ),
|
|
$ LDA )
|
|
CALL ZTRSM( 'Left', 'Upper', 'Conjugate transpose',
|
|
$ 'Non-unit', JB, N-J-JB+1, CONE, A( J, J ),
|
|
$ LDA, A( J, J+JB ), LDA )
|
|
END IF
|
|
10 CONTINUE
|
|
*
|
|
ELSE
|
|
*
|
|
* Compute the Cholesky factorization A = L*L'.
|
|
*
|
|
DO 20 J = 1, N, NB
|
|
*
|
|
* Update and factorize the current diagonal block and test
|
|
* for non-positive-definiteness.
|
|
*
|
|
JB = MIN( NB, N-J+1 )
|
|
CALL ZHERK( 'Lower', 'No transpose', JB, J-1, -ONE,
|
|
$ A( J, 1 ), LDA, ONE, A( J, J ), LDA )
|
|
CALL ZPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
|
|
IF( INFO.NE.0 )
|
|
$ GO TO 30
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Compute the current block column.
|
|
*
|
|
CALL ZGEMM( 'No transpose', 'Conjugate transpose',
|
|
$ N-J-JB+1, JB, J-1, -CONE, A( J+JB, 1 ),
|
|
$ LDA, A( J, 1 ), LDA, CONE, A( J+JB, J ),
|
|
$ LDA )
|
|
CALL ZTRSM( 'Right', 'Lower', 'Conjugate transpose',
|
|
$ 'Non-unit', N-J-JB+1, JB, CONE, A( J, J ),
|
|
$ LDA, A( J+JB, J ), LDA )
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
END IF
|
|
GO TO 40
|
|
*
|
|
30 CONTINUE
|
|
INFO = INFO + J - 1
|
|
*
|
|
40 CONTINUE
|
|
RETURN
|
|
*
|
|
* End of ZPOTRF
|
|
*
|
|
END
|
|
|
|
SUBROUTINE ZTRTI2( UPLO, DIAG, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER DIAG, UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
COMPLEX*16 A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZTRTI2 computes the inverse of a complex upper or lower triangular
|
|
* matrix.
|
|
*
|
|
* This is the Level 2 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* Specifies whether the matrix A is upper or lower triangular.
|
|
* = 'U': Upper triangular
|
|
* = 'L': Lower triangular
|
|
*
|
|
* DIAG (input) CHARACTER*1
|
|
* Specifies whether or not the matrix A is unit triangular.
|
|
* = 'N': Non-unit triangular
|
|
* = 'U': Unit triangular
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
|
|
* On entry, the triangular matrix A. If UPLO = 'U', the
|
|
* leading n by n upper triangular part of the array A contains
|
|
* the upper triangular matrix, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading n by n lower triangular part of the array A contains
|
|
* the lower triangular matrix, and the strictly upper
|
|
* triangular part of A is not referenced. If DIAG = 'U', the
|
|
* diagonal elements of A are also not referenced and are
|
|
* assumed to be 1.
|
|
*
|
|
* On exit, the (triangular) inverse of the original matrix, in
|
|
* the same storage format.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -k, the k-th argument had an illegal value
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
COMPLEX*16 ONE
|
|
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL NOUNIT, UPPER
|
|
INTEGER J
|
|
COMPLEX*16 AJJ
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
EXTERNAL LSAME
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZSCAL, ZTRMV
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
NOUNIT = LSAME( DIAG, 'N' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
|
|
INFO = -2
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZTRTI2', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute inverse of upper triangular matrix.
|
|
*
|
|
DO 10 J = 1, N
|
|
IF( NOUNIT ) THEN
|
|
A( J, J ) = ONE / A( J, J )
|
|
AJJ = -A( J, J )
|
|
ELSE
|
|
AJJ = -ONE
|
|
END IF
|
|
*
|
|
* Compute elements 1:j-1 of j-th column.
|
|
*
|
|
CALL ZTRMV( 'Upper', 'No transpose', DIAG, J-1, A, LDA,
|
|
$ A( 1, J ), 1 )
|
|
CALL ZSCAL( J-1, AJJ, A( 1, J ), 1 )
|
|
10 CONTINUE
|
|
ELSE
|
|
*
|
|
* Compute inverse of lower triangular matrix.
|
|
*
|
|
DO 20 J = N, 1, -1
|
|
IF( NOUNIT ) THEN
|
|
A( J, J ) = ONE / A( J, J )
|
|
AJJ = -A( J, J )
|
|
ELSE
|
|
AJJ = -ONE
|
|
END IF
|
|
IF( J.LT.N ) THEN
|
|
*
|
|
* Compute elements j+1:n of j-th column.
|
|
*
|
|
CALL ZTRMV( 'Lower', 'No transpose', DIAG, N-J,
|
|
$ A( J+1, J+1 ), LDA, A( J+1, J ), 1 )
|
|
CALL ZSCAL( N-J, AJJ, A( J+1, J ), 1 )
|
|
END IF
|
|
20 CONTINUE
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of ZTRTI2
|
|
*
|
|
END
|
|
|
|
SUBROUTINE ZTRTRI( UPLO, DIAG, N, A, LDA, INFO )
|
|
*
|
|
* -- LAPACK routine (version 3.0) --
|
|
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
|
* Courant Institute, Argonne National Lab, and Rice University
|
|
* September 30, 1994
|
|
*
|
|
* .. Scalar Arguments ..
|
|
CHARACTER DIAG, UPLO
|
|
INTEGER INFO, LDA, N
|
|
* ..
|
|
* .. Array Arguments ..
|
|
COMPLEX*16 A( LDA, * )
|
|
* ..
|
|
*
|
|
* Purpose
|
|
* =======
|
|
*
|
|
* ZTRTRI computes the inverse of a complex upper or lower triangular
|
|
* matrix A.
|
|
*
|
|
* This is the Level 3 BLAS version of the algorithm.
|
|
*
|
|
* Arguments
|
|
* =========
|
|
*
|
|
* UPLO (input) CHARACTER*1
|
|
* = 'U': A is upper triangular;
|
|
* = 'L': A is lower triangular.
|
|
*
|
|
* DIAG (input) CHARACTER*1
|
|
* = 'N': A is non-unit triangular;
|
|
* = 'U': A is unit triangular.
|
|
*
|
|
* N (input) INTEGER
|
|
* The order of the matrix A. N >= 0.
|
|
*
|
|
* A (input/output) COMPLEX*16 array, dimension (LDA,N)
|
|
* On entry, the triangular matrix A. If UPLO = 'U', the
|
|
* leading N-by-N upper triangular part of the array A contains
|
|
* the upper triangular matrix, and the strictly lower
|
|
* triangular part of A is not referenced. If UPLO = 'L', the
|
|
* leading N-by-N lower triangular part of the array A contains
|
|
* the lower triangular matrix, and the strictly upper
|
|
* triangular part of A is not referenced. If DIAG = 'U', the
|
|
* diagonal elements of A are also not referenced and are
|
|
* assumed to be 1.
|
|
* On exit, the (triangular) inverse of the original matrix, in
|
|
* the same storage format.
|
|
*
|
|
* LDA (input) INTEGER
|
|
* The leading dimension of the array A. LDA >= max(1,N).
|
|
*
|
|
* INFO (output) INTEGER
|
|
* = 0: successful exit
|
|
* < 0: if INFO = -i, the i-th argument had an illegal value
|
|
* > 0: if INFO = i, A(i,i) is exactly zero. The triangular
|
|
* matrix is singular and its inverse can not be computed.
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
COMPLEX*16 ONE, ZERO
|
|
PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
|
|
$ ZERO = ( 0.0D+0, 0.0D+0 ) )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL NOUNIT, UPPER
|
|
INTEGER J, JB, NB, NN
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL LSAME
|
|
INTEGER ILAENV
|
|
EXTERNAL LSAME, ILAENV
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL XERBLA, ZTRMM, ZTRSM, ZTRTI2
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC MAX, MIN
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Test the input parameters.
|
|
*
|
|
INFO = 0
|
|
UPPER = LSAME( UPLO, 'U' )
|
|
NOUNIT = LSAME( DIAG, 'N' )
|
|
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
|
INFO = -1
|
|
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
|
|
INFO = -2
|
|
ELSE IF( N.LT.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
|
INFO = -5
|
|
END IF
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'ZTRTRI', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Quick return if possible
|
|
*
|
|
IF( N.EQ.0 )
|
|
$ RETURN
|
|
*
|
|
* Check for singularity if non-unit.
|
|
*
|
|
IF( NOUNIT ) THEN
|
|
DO 10 INFO = 1, N
|
|
IF( A( INFO, INFO ).EQ.ZERO )
|
|
$ RETURN
|
|
10 CONTINUE
|
|
INFO = 0
|
|
END IF
|
|
*
|
|
* Determine the block size for this environment.
|
|
*
|
|
NB = ILAENV( 1, 'ZTRTRI', UPLO // DIAG, N, -1, -1, -1 )
|
|
IF( NB.LE.1 .OR. NB.GE.N ) THEN
|
|
*
|
|
* Use unblocked code
|
|
*
|
|
CALL ZTRTI2( UPLO, DIAG, N, A, LDA, INFO )
|
|
ELSE
|
|
*
|
|
* Use blocked code
|
|
*
|
|
IF( UPPER ) THEN
|
|
*
|
|
* Compute inverse of upper triangular matrix
|
|
*
|
|
DO 20 J = 1, N, NB
|
|
JB = MIN( NB, N-J+1 )
|
|
*
|
|
* Compute rows 1:j-1 of current block column
|
|
*
|
|
CALL ZTRMM( 'Left', 'Upper', 'No transpose', DIAG, J-1,
|
|
$ JB, ONE, A, LDA, A( 1, J ), LDA )
|
|
CALL ZTRSM( 'Right', 'Upper', 'No transpose', DIAG, J-1,
|
|
$ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA )
|
|
*
|
|
* Compute inverse of current diagonal block
|
|
*
|
|
CALL ZTRTI2( 'Upper', DIAG, JB, A( J, J ), LDA, INFO )
|
|
20 CONTINUE
|
|
ELSE
|
|
*
|
|
* Compute inverse of lower triangular matrix
|
|
*
|
|
NN = ( ( N-1 ) / NB )*NB + 1
|
|
DO 30 J = NN, 1, -NB
|
|
JB = MIN( NB, N-J+1 )
|
|
IF( J+JB.LE.N ) THEN
|
|
*
|
|
* Compute rows j+jb:n of current block column
|
|
*
|
|
CALL ZTRMM( 'Left', 'Lower', 'No transpose', DIAG,
|
|
$ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA,
|
|
$ A( J+JB, J ), LDA )
|
|
CALL ZTRSM( 'Right', 'Lower', 'No transpose', DIAG,
|
|
$ N-J-JB+1, JB, -ONE, A( J, J ), LDA,
|
|
$ A( J+JB, J ), LDA )
|
|
END IF
|
|
*
|
|
* Compute inverse of current diagonal block
|
|
*
|
|
CALL ZTRTI2( 'Lower', DIAG, JB, A( J, J ), LDA, INFO )
|
|
30 CONTINUE
|
|
END IF
|
|
END IF
|
|
*
|
|
RETURN
|
|
*
|
|
* End of ZTRTRI
|
|
*
|
|
END
|