ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c The following LAPACK and BLAS subroutines and functions are contained
c in this file, in the order listed:
c
c dpotri
c dlauum
c dlauu2
c dtrtri
c dtrti2
c xerbla
c lsame
c ilaenv
c dpotrf
c dpotf2
c dgemm
c dsyrk
c dtrmm
c dtrsm
c dgemv
c dscal
c dtrmv
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
SUBROUTINE DPOTRI( UPLO, N, A, LDA, INFO )
*
* -- LAPACK routine (version 1.0b) --
* 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
* =======
*
* DPOTRI computes the inverse of a real symmetric positive definite
* matrix A using the Cholesky factorization A = U'*U or A = L*L'
* computed by DPOTRF.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the factor stored in A is upper or lower
* triangular.
* = '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 triangular factor U or L from the Cholesky
* factorization A = U'*U or A = L*L', as computed by DPOTRF.
* On exit, the upper or lower triangle of the (symmetric)
* inverse of A, overwriting the input factor U or 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 (k,k) element of the factor U or L is
* zero, and the inverse could not be computed.
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DLAUUM, DTRTRI, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF( .NOT.LSAME( UPLO, 'U' ) .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( 'DPOTRI', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
* Invert the triangular Cholesky factor U or L.
*
CALL DTRTRI( UPLO, 'Non-unit', N, A, LDA, INFO )
IF( INFO.GT.0 )
$ RETURN
*
* Form inv(U)*inv(U)' or inv(L)'*inv(L).
*
CALL DLAUUM( UPLO, N, A, LDA, INFO )
*
RETURN
*
* End of DPOTRI
*
END
c
c
SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO )
*
* -- LAPACK auxiliary routine (version 1.0b) --
* 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
* =======
*
* DLAUUM 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
* the array A.
*
* If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
* overwriting the factor U in A.
* If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
* overwriting the factor L in A.
*
* This is the blocked form of the algorithm, calling Level 3 BLAS.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the triangular factor stored in the array A
* is upper or lower triangular:
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* N (input) INTEGER
* The order of the triangular factor U or L. N >= 0.
*
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
* On entry, the triangular factor U or L.
* On exit, if UPLO = 'U', the upper triangle of A is
* overwritten with the upper triangle of the product U * U';
* if UPLO = 'L', the lower triangle of A is overwritten with
* the lower triangle of the product 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
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, IB, NB
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
* ..
* .. External Subroutines ..
EXTERNAL DGEMM, DLAUU2, DSYRK, DTRMM, 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( 'DLAUUM', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
* Determine the block size for this environment.
*
NB = ILAENV( 1, 'DLAUUM', UPLO, N, -1, -1, -1 )
*
IF( NB.LE.1 .OR. NB.GE.N ) THEN
*
* Use unblocked code
*
CALL DLAUU2( UPLO, N, A, LDA, INFO )
ELSE
*
* Use blocked code
*
IF( UPPER ) THEN
*
* Compute the product U * U'.
*
DO 10 I = 1, N, NB
IB = MIN( NB, N-I+1 )
CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit',
$ I-1, IB, ONE, A( I, I ), LDA, A( 1, I ),
$ LDA )
CALL DLAUU2( 'Upper', IB, A( I, I ), LDA, INFO )
IF( I+IB.LE.N ) THEN
CALL DGEMM( 'No transpose', 'Transpose', I-1, IB,
$ N-I-IB+1, ONE, A( 1, I+IB ), LDA,
$ A( I, I+IB ), LDA, ONE, A( 1, I ), LDA )
CALL DSYRK( 'Upper', 'No transpose', IB, N-I-IB+1,
$ ONE, A( I, I+IB ), LDA, ONE, A( I, I ),
$ LDA )
END IF
10 CONTINUE
ELSE
*
* Compute the product L' * L.
*
DO 20 I = 1, N, NB
IB = MIN( NB, N-I+1 )
CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', IB,
$ I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA )
CALL DLAUU2( 'Lower', IB, A( I, I ), LDA, INFO )
IF( I+IB.LE.N ) THEN
CALL DGEMM( 'Transpose', 'No transpose', IB, I-1,
$ N-I-IB+1, ONE, A( I+IB, I ), LDA,
$ A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA )
CALL DSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, ONE,
$ A( I+IB, I ), LDA, ONE, A( I, I ), LDA )
END IF
20 CONTINUE
END IF
END IF
*
RETURN
*
* End of DLAUUM
*
END
c
c
c
SUBROUTINE DLAUU2( UPLO, N, A, LDA, INFO )
*
* -- LAPACK auxiliary routine (version 1.0b) --
* 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
* =======
*
* DLAUU2 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
* the array A.
*
* If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
* overwriting the factor U in A.
* If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
* overwriting the factor L in A.
*
* This is the unblocked form of the algorithm, calling Level 2 BLAS.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* Specifies whether the triangular factor stored in the array A
* is upper or lower triangular:
* = 'U': Upper triangular
* = 'L': Lower triangular
*
* N (input) INTEGER
* The order of the triangular factor U or L. N >= 0.
*
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
* On entry, the triangular factor U or L.
* On exit, if UPLO = 'U', the upper triangle of A is
* overwritten with the upper triangle of the product U * U';
* if UPLO = 'L', the lower triangle of A is overwritten with
* the lower triangle of the product 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
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I
DOUBLE PRECISION AII
* ..
* .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DDOT
EXTERNAL LSAME, DDOT
* ..
* .. External Subroutines ..
EXTERNAL DGEMV, DSCAL, 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( LDA.LT.MAX( 1, N ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DLAUU2', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
IF( UPPER ) THEN
*
* Compute the product U * U'.
*
DO 10 I = 1, N
AII = A( I, I )
IF( I.LT.N ) THEN
A( I, I ) = DDOT( N-I+1, A( I, I ), LDA, A( I, I ), LDA )
CALL DGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
$ LDA, A( I, I+1 ), LDA, AII, A( 1, I ), 1 )
ELSE
CALL DSCAL( I, AII, A( 1, I ), 1 )
END IF
10 CONTINUE
*
ELSE
*
* Compute the product L' * L.
*
DO 20 I = 1, N
AII = A( I, I )
IF( I.LT.N ) THEN
A( I, I ) = DDOT( N-I+1, A( I, I ), 1, A( I, I ), 1 )
CALL DGEMV( 'Transpose', N-I, I-1, ONE, A( I+1, 1 ), LDA,
$ A( I+1, I ), 1, AII, A( I, 1 ), LDA )
ELSE
CALL DSCAL( I, AII, A( I, 1 ), LDA )
END IF
20 CONTINUE
END IF
*
RETURN
*
* End of DLAUU2
*
END
c
c
c
SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO )
*
* -- LAPACK routine (version 1.0b) --
* 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
* =======
*
* 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
* 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, A(k,k) is exactly zero. The triangular
* matrix is singular and its inverse can not be computed.
* < 0: if INFO = -k, the k-th argument had an illegal value
*
* =====================================================================
*
* .. 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
c
c
c
SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
*
* -- LAPACK routine (version 1.0b) --
* 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
c
c
c
SUBROUTINE XERBLA( SRNAME, INFO )
*
* -- LAPACK auxiliary routine (version 1.0b) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 29, 1992
*
* .. 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
c
c
c
LOGICAL FUNCTION LSAME( CA, CB )
*
* -- LAPACK auxiliary routine (version 1.0b) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 29, 1992
*
* .. 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
c
c
c
INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,
$ N4 )
*
* -- LAPACK auxiliary routine (preliminary version) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 20, 1992
*
* .. Scalar Arguments ..
CHARACTER*( * ) NAME, OPTS
INTEGER ISPEC, N1, N2, N3, N4
* ..
*
* Purpose
* =======
*
* ILAENV is called from the LAPACK routines to choose problem-dependent
* parameters for the local environment. See ISPEC for a description of
* the parameters.
*
* This version provides a set of parameters which should give good,
* but not optimal, performance on many of the currently available
* computers. Users are encouraged to modify this subroutine to set
* the tuning parameters for their particular machine using the option
* and problem size information in the arguments.
*
* This routine will not function correctly if it is converted to all
* lower case. Converting it to all upper case is allowed.
*
* Arguments
* =========
*
* ISPEC (input) INTEGER
* Specifies the parameter to be returned as the value of
* ILAENV.
* = 1: the optimal blocksize; if this value is 1, an unblocked
* algorithm will give the best performance.
* = 2: the minimum block size for which the block routine
* should be used; if the usable block size is less than
* this value, an unblocked routine should be used.
* = 3: the crossover point (in a block routine, for N less
* than this value, an unblocked routine should be used)
* = 4: the number of shifts, used in the nonsymmetric
* eigenvalue routines
* = 5: the minimum column dimension for blocking to be used;
* rectangular blocks must have dimension at least k by m,
* where k is given by ILAENV(2,...) and m by ILAENV(5,...)
* = 6: the crossover point for the SVD (when reducing an m by n
* matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
* this value, a QR factorization is used first to reduce
* the matrix to a triangular form.)
* = 7: the number of processors
* = 8: the crossover point for the multishift QR and QZ methods
* for nonsymmetric eigenvalue problems.
*
* NAME (input) CHARACTER*(*)
* The name of the calling subroutine, in either upper case or
* lower case.
*
* OPTS (input) CHARACTER*(*)
* The character options to the subroutine NAME, concatenated
* into a single character string. For example, UPLO = 'U',
* TRANS = 'T', and DIAG = 'N' for a triangular routine would
* be specified as OPTS = 'UTN'.
*
* N1 (input) INTEGER
* N2 (input) INTEGER
* N3 (input) INTEGER
* N4 (input) INTEGER
* Problem dimensions for the subroutine NAME; these may not all
* be required.
*
* (ILAENV) (output) INTEGER
* >= 0: the value of the parameter specified by ISPEC
* < 0: if ILAENV = -k, the k-th argument had an illegal value.
*
* Further Details
* ===============
*
* The following conventions have been used when calling ILAENV from the
* LAPACK routines:
* 1) OPTS is a concatenation of all of the character options to
* subroutine NAME, in the same order that they appear in the
* argument list for NAME, even if they are not used in determining
* the value of the parameter specified by ISPEC.
* 2) The problem dimensions N1, N2, N3, N4 are specified in the order
* that they appear in the argument list for NAME. N1 is used
* first, N2 second, and so on, and unused problem dimensions are
* passed a value of -1.
* 3) The parameter value returned by ILAENV is checked for validity in
* the calling subroutine. For example, ILAENV is used to retrieve
* the optimal blocksize for STRTRI as follows:
*
* NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
* IF( NB.LE.1 ) NB = MAX( 1, N )
*
* =====================================================================
*
* .. Local Scalars ..
LOGICAL CNAME, SNAME
CHARACTER*1 C1
CHARACTER*2 C2, C4
CHARACTER*3 C3
CHARACTER*6 SUBNAM
INTEGER I, IC, IZ, NB, NBMIN, NX
* ..
* .. Intrinsic Functions ..
INTRINSIC CHAR, ICHAR, INT, MIN, REAL
* ..
* .. Executable Statements ..
*
GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC
*
* Invalid value for ISPEC
*
ILAENV = -1
RETURN
*
100 CONTINUE
*
* Convert NAME to upper case if the first character is lower case.
*
ILAENV = 1
SUBNAM = NAME
IC = ICHAR( SUBNAM( 1:1 ) )
IZ = ICHAR( 'Z' )
IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN
*
* ASCII character set
*
IF( IC.GE.97 .AND. IC.LE.122 ) THEN
SUBNAM( 1:1 ) = CHAR( IC-32 )
DO 10 I = 2, 6
IC = ICHAR( SUBNAM( I:I ) )
IF( IC.GE.97 .AND. IC.LE.122 )
$ SUBNAM( I:I ) = CHAR( IC-32 )
10 CONTINUE
END IF
*
ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN
*
* EBCDIC character set
*
IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
$ ( IC.GE.145 .AND. IC.LE.153 ) .OR.
$ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN
SUBNAM( 1:1 ) = CHAR( IC+64 )
DO 20 I = 2, 6
IC = ICHAR( SUBNAM( I:I ) )
IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.
$ ( IC.GE.145 .AND. IC.LE.153 ) .OR.
$ ( IC.GE.162 .AND. IC.LE.169 ) )
$ SUBNAM( I:I ) = CHAR( IC+64 )
20 CONTINUE
END IF
*
ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN
*
* Prime machines: ASCII+128
*
IF( IC.GE.225 .AND. IC.LE.250 ) THEN
SUBNAM( 1:1 ) = CHAR( IC-32 )
DO 30 I = 2, 6
IC = ICHAR( SUBNAM( I:I ) )
IF( IC.GE.225 .AND. IC.LE.250 )
$ SUBNAM( I:I ) = CHAR( IC-32 )
30 CONTINUE
END IF
END IF
*
C1 = SUBNAM( 1:1 )
SNAME = C1.EQ.'S' .OR. C1.EQ.'D'
CNAME = C1.EQ.'C' .OR. C1.EQ.'Z'
IF( .NOT.( CNAME .OR. SNAME ) )
$ RETURN
C2 = SUBNAM( 2:3 )
C3 = SUBNAM( 4:6 )
C4 = C3( 2:3 )
*
GO TO ( 110, 200, 300 ) ISPEC
*
110 CONTINUE
*
* ISPEC = 1: block size
*
* In these examples, separate code is provided for setting NB for
* real and complex. We assume that NB will take the same value in
* single or double precision.
*
NB = 1
*
IF( C2.EQ.'GE' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
NB = 64
ELSE
NB = 64
END IF
ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
$ C3.EQ.'QLF' ) THEN
IF( SNAME ) THEN
NB = 32
ELSE
NB = 32
END IF
ELSE IF( C3.EQ.'HRD' ) THEN
IF( SNAME ) THEN
NB = 32
ELSE
NB = 32
END IF
ELSE IF( C3.EQ.'BRD' ) THEN
IF( SNAME ) THEN
NB = 32
ELSE
NB = 32
END IF
ELSE IF( C3.EQ.'TRI' ) THEN
IF( SNAME ) THEN
NB = 64
ELSE
NB = 64
END IF
END IF
ELSE IF( C2.EQ.'PO' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
NB = 64
ELSE
NB = 64
END IF
END IF
ELSE IF( C2.EQ.'SY' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
NB = 64
ELSE
NB = 64
END IF
ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
NB = 1
ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN
NB = 64
END IF
ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
IF( C3.EQ.'TRF' ) THEN
NB = 64
ELSE IF( C3.EQ.'TRD' ) THEN
NB = 1
ELSE IF( C3.EQ.'GST' ) THEN
NB = 64
END IF
ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
IF( C3( 1:1 ).EQ.'G' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NB = 32
END IF
ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NB = 32
END IF
END IF
ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
IF( C3( 1:1 ).EQ.'G' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NB = 32
END IF
ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NB = 32
END IF
END IF
ELSE IF( C2.EQ.'GB' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
IF( N4.LE.64 ) THEN
NB = 1
ELSE
NB = 32
END IF
ELSE
IF( N4.LE.64 ) THEN
NB = 1
ELSE
NB = 32
END IF
END IF
END IF
ELSE IF( C2.EQ.'PB' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
IF( N2.LE.64 ) THEN
NB = 1
ELSE
NB = 32
END IF
ELSE
IF( N2.LE.64 ) THEN
NB = 1
ELSE
NB = 32
END IF
END IF
END IF
ELSE IF( C2.EQ.'TR' ) THEN
IF( C3.EQ.'TRI' ) THEN
IF( SNAME ) THEN
NB = 64
ELSE
NB = 64
END IF
END IF
ELSE IF( C2.EQ.'LA' ) THEN
IF( C3.EQ.'UUM' ) THEN
IF( SNAME ) THEN
NB = 64
ELSE
NB = 64
END IF
END IF
ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN
IF( C3.EQ.'EBZ' ) THEN
NB = 1
END IF
END IF
ILAENV = NB
RETURN
*
200 CONTINUE
*
* ISPEC = 2: minimum block size
*
NBMIN = 2
IF( C2.EQ.'GE' ) THEN
IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
$ C3.EQ.'QLF' ) THEN
IF( SNAME ) THEN
NBMIN = 2
ELSE
NBMIN = 2
END IF
ELSE IF( C3.EQ.'HRD' ) THEN
IF( SNAME ) THEN
NBMIN = 2
ELSE
NBMIN = 2
END IF
ELSE IF( C3.EQ.'BRD' ) THEN
IF( SNAME ) THEN
NBMIN = 2
ELSE
NBMIN = 2
END IF
ELSE IF( C3.EQ.'TRI' ) THEN
IF( SNAME ) THEN
NBMIN = 2
ELSE
NBMIN = 2
END IF
END IF
ELSE IF( C2.EQ.'SY' ) THEN
IF( C3.EQ.'TRF' ) THEN
IF( SNAME ) THEN
NBMIN = 2
ELSE
NBMIN = 2
END IF
ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN
NBMIN = 2
END IF
ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
IF( C3.EQ.'TRD' ) THEN
NBMIN = 2
END IF
ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
IF( C3( 1:1 ).EQ.'G' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NBMIN = 2
END IF
ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NBMIN = 2
END IF
END IF
ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
IF( C3( 1:1 ).EQ.'G' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NBMIN = 2
END IF
ELSE IF( C3( 1:1 ).EQ.'M' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NBMIN = 2
END IF
END IF
END IF
ILAENV = NBMIN
RETURN
*
300 CONTINUE
*
* ISPEC = 3: crossover point
*
NX = 0
IF( C2.EQ.'GE' ) THEN
IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.
$ C3.EQ.'QLF' ) THEN
IF( SNAME ) THEN
NX = 128
ELSE
NX = 128
END IF
ELSE IF( C3.EQ.'HRD' ) THEN
IF( SNAME ) THEN
NX = 128
ELSE
NX = 128
END IF
ELSE IF( C3.EQ.'BRD' ) THEN
IF( SNAME ) THEN
NX = 128
ELSE
NX = 128
END IF
END IF
ELSE IF( C2.EQ.'SY' ) THEN
IF( SNAME .AND. C3.EQ.'TRD' ) THEN
NX = 1
END IF
ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN
IF( C3.EQ.'TRD' ) THEN
NX = 1
END IF
ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN
IF( C3( 1:1 ).EQ.'G' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NX = 128
END IF
END IF
ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN
IF( C3( 1:1 ).EQ.'G' ) THEN
IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.
$ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.
$ C4.EQ.'BR' ) THEN
NX = 128
END IF
END IF
END IF
ILAENV = NX
RETURN
*
400 CONTINUE
*
* ISPEC = 4: number of shifts (used by xHSEQR)
*
ILAENV = 6
RETURN
*
500 CONTINUE
*
* ISPEC = 5: minimum column dimension (not used)
*
ILAENV = 2
RETURN
*
600 CONTINUE
*
* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
*
ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 )
RETURN
*
700 CONTINUE
*
* ISPEC = 7: number of processors (not used)
*
ILAENV = 1
RETURN
*
800 CONTINUE
*
* ISPEC = 8: crossover point for multishift (used by xHSEQR)
*
ILAENV = 50
RETURN
*
* End of ILAENV
*
END
c
c
c
SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
*
* -- LAPACK routine (version 1.0b) --
* 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
* =======
*
* DPOTRF 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 block version of the algorithm, calling Level 3 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
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
c
c
c
SUBROUTINE DPOTF2( UPLO, N, A, LDA, INFO )
*
* -- LAPACK routine (version 1.0b) --
* 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
c
c
c
*
************************************************************************
*
* File of the DOUBLE PRECISION Level-3 BLAS.
* ==========================================
*
* SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
* $ BETA, C, LDC )
*
* SUBROUTINE DSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB,
* $ BETA, C, LDC )
*
* SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
* $ BETA, C, LDC )
*
* SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB,
* $ BETA, C, LDC )
*
* SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
* $ B, LDB )
*
* SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
* $ B, LDB )
*
* See:
*
* Dongarra J. J., Du Croz J. J., Duff I. and Hammarling S.
* A set of Level 3 Basic Linear Algebra Subprograms. Technical
* Memorandum No.88 (Revision 1), Mathematics and Computer Science
* Division, Argonne National Laboratory, 9700 South Cass Avenue,
* Argonne, Illinois 60439.
*
*
************************************************************************
*
SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
$ BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER*1 TRANSA, TRANSB
INTEGER M, N, K, LDA, LDB, LDC
DOUBLE PRECISION ALPHA, BETA
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* DGEMM performs one of the matrix-matrix operations
*
* C := alpha*op( A )*op( B ) + beta*C,
*
* where op( X ) is one of
*
* op( X ) = X or op( X ) = X',
*
* alpha and beta are scalars, and A, B and C are matrices, with op( A )
* an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
*
* Parameters
* ==========
*
* TRANSA - CHARACTER*1.
* On entry, TRANSA specifies the form of op( A ) to be used in
* the matrix multiplication as follows:
*
* TRANSA = 'N' or 'n', op( A ) = A.
*
* TRANSA = 'T' or 't', op( A ) = A'.
*
* TRANSA = 'C' or 'c', op( A ) = A'.
*
* Unchanged on exit.
*
* TRANSB - CHARACTER*1.
* On entry, TRANSB specifies the form of op( B ) to be used in
* the matrix multiplication as follows:
*
* TRANSB = 'N' or 'n', op( B ) = B.
*
* TRANSB = 'T' or 't', op( B ) = B'.
*
* TRANSB = 'C' or 'c', op( B ) = B'.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix
* op( A ) and of the matrix C. M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix
* op( B ) and the number of columns of the matrix C. N must be
* at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry, K specifies the number of columns of the matrix
* op( A ) and the number of rows of the matrix op( B ). K must
* be at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION.
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
* k when TRANSA = 'N' or 'n', and is m otherwise.
* Before entry with TRANSA = 'N' or 'n', the leading m by k
* part of the array A must contain the matrix A, otherwise
* the leading k by m part of the array A must contain the
* matrix A.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When TRANSA = 'N' or 'n' then
* LDA must be at least max( 1, m ), otherwise LDA must be at
* least max( 1, k ).
* Unchanged on exit.
*
* B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
* n when TRANSB = 'N' or 'n', and is k otherwise.
* Before entry with TRANSB = 'N' or 'n', the leading k by n
* part of the array B must contain the matrix B, otherwise
* the leading n by k part of the array B must contain the
* matrix B.
* Unchanged on exit.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. When TRANSB = 'N' or 'n' then
* LDB must be at least max( 1, k ), otherwise LDB must be at
* least max( 1, n ).
* Unchanged on exit.
*
* BETA - DOUBLE PRECISION.
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then C need not be set on input.
* Unchanged on exit.
*
* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
* Before entry, the leading m by n part of the array C must
* contain the matrix C, except when beta is zero, in which
* case C need not be set on entry.
* On exit, the array C is overwritten by the m by n matrix
* ( alpha*op( A )*op( B ) + beta*C ).
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Local Scalars ..
LOGICAL NOTA, NOTB
INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB
DOUBLE PRECISION TEMP
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Executable Statements ..
*
* Set NOTA and NOTB as true if A and B respectively are not
* transposed and set NROWA, NCOLA and NROWB as the number of rows
* and columns of A and the number of rows of B respectively.
*
NOTA = LSAME( TRANSA, 'N' )
NOTB = LSAME( TRANSB, 'N' )
IF( NOTA )THEN
NROWA = M
NCOLA = K
ELSE
NROWA = K
NCOLA = M
END IF
IF( NOTB )THEN
NROWB = K
ELSE
NROWB = N
END IF
*
* Test the input parameters.
*
INFO = 0
IF( ( .NOT.NOTA ).AND.
$ ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.NOTB ).AND.
$ ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
$ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN
INFO = 2
ELSE IF( M .LT.0 )THEN
INFO = 3
ELSE IF( N .LT.0 )THEN
INFO = 4
ELSE IF( K .LT.0 )THEN
INFO = 5
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 8
ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
INFO = 10
ELSE IF( LDC.LT.MAX( 1, M ) )THEN
INFO = 13
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'DGEMM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* And if alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
IF( BETA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40, J = 1, N
DO 30, I = 1, M
C( I, J ) = BETA*C( I, J )
30 CONTINUE
40 CONTINUE
END IF
RETURN
END IF
*
* Start the operations.
*
IF( NOTB )THEN
IF( NOTA )THEN
*
* Form C := alpha*A*B + beta*C.
*
DO 90, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 50, I = 1, M
C( I, J ) = ZERO
50 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 60, I = 1, M
C( I, J ) = BETA*C( I, J )
60 CONTINUE
END IF
DO 80, L = 1, K
IF( B( L, J ).NE.ZERO )THEN
TEMP = ALPHA*B( L, J )
DO 70, I = 1, M
C( I, J ) = C( I, J ) + TEMP*A( I, L )
70 CONTINUE
END IF
80 CONTINUE
90 CONTINUE
ELSE
*
* Form C := alpha*A'*B + beta*C
*
DO 120, J = 1, N
DO 110, I = 1, M
TEMP = ZERO
DO 100, L = 1, K
TEMP = TEMP + A( L, I )*B( L, J )
100 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
110 CONTINUE
120 CONTINUE
END IF
ELSE
IF( NOTA )THEN
*
* Form C := alpha*A*B' + beta*C
*
DO 170, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 130, I = 1, M
C( I, J ) = ZERO
130 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 140, I = 1, M
C( I, J ) = BETA*C( I, J )
140 CONTINUE
END IF
DO 160, L = 1, K
IF( B( J, L ).NE.ZERO )THEN
TEMP = ALPHA*B( J, L )
DO 150, I = 1, M
C( I, J ) = C( I, J ) + TEMP*A( I, L )
150 CONTINUE
END IF
160 CONTINUE
170 CONTINUE
ELSE
*
* Form C := alpha*A'*B' + beta*C
*
DO 200, J = 1, N
DO 190, I = 1, M
TEMP = ZERO
DO 180, L = 1, K
TEMP = TEMP + A( L, I )*B( J, L )
180 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
190 CONTINUE
200 CONTINUE
END IF
END IF
*
RETURN
*
* End of DGEMM .
*
END
*
************************************************************************
*
SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,
$ BETA, C, LDC )
* .. Scalar Arguments ..
CHARACTER*1 UPLO, TRANS
INTEGER N, K, LDA, LDC
DOUBLE PRECISION ALPHA, BETA
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), C( LDC, * )
* ..
*
* Purpose
* =======
*
* DSYRK performs one of the symmetric rank k operations
*
* C := alpha*A*A' + beta*C,
*
* or
*
* C := alpha*A'*A + beta*C,
*
* where alpha and beta are scalars, C is an n by n symmetric matrix
* and A is an n by k matrix in the first case and a k by n matrix
* in the second case.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the upper or lower
* triangular part of the array C is to be referenced as
* follows:
*
* UPLO = 'U' or 'u' Only the upper triangular part of C
* is to be referenced.
*
* UPLO = 'L' or 'l' Only the lower triangular part of C
* is to be referenced.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' C := alpha*A*A' + beta*C.
*
* TRANS = 'T' or 't' C := alpha*A'*A + beta*C.
*
* TRANS = 'C' or 'c' C := alpha*A'*A + beta*C.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix C. N must be
* at least zero.
* Unchanged on exit.
*
* K - INTEGER.
* On entry with TRANS = 'N' or 'n', K specifies the number
* of columns of the matrix A, and on entry with
* TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
* of rows of the matrix A. K must be at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION.
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
* k when TRANS = 'N' or 'n', and is n otherwise.
* Before entry with TRANS = 'N' or 'n', the leading n by k
* part of the array A must contain the matrix A, otherwise
* the leading k by n part of the array A must contain the
* matrix A.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When TRANS = 'N' or 'n'
* then LDA must be at least max( 1, n ), otherwise LDA must
* be at least max( 1, k ).
* Unchanged on exit.
*
* BETA - DOUBLE PRECISION.
* On entry, BETA specifies the scalar beta.
* Unchanged on exit.
*
* C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array C must contain the upper
* triangular part of the symmetric matrix and the strictly
* lower triangular part of C is not referenced. On exit, the
* upper triangular part of the array C is overwritten by the
* upper triangular part of the updated matrix.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array C must contain the lower
* triangular part of the symmetric matrix and the strictly
* upper triangular part of C is not referenced. On exit, the
* lower triangular part of the array C is overwritten by the
* lower triangular part of the updated matrix.
*
* LDC - INTEGER.
* On entry, LDC specifies the first dimension of C as declared
* in the calling (sub) program. LDC must be at least
* max( 1, n ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Local Scalars ..
LOGICAL UPPER
INTEGER I, INFO, J, L, NROWA
DOUBLE PRECISION TEMP
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
IF( LSAME( TRANS, 'N' ) )THEN
NROWA = N
ELSE
NROWA = K
END IF
UPPER = LSAME( UPLO, 'U' )
*
INFO = 0
IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.
$ ( .NOT.LSAME( TRANS, 'T' ) ).AND.
$ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN
INFO = 2
ELSE IF( N .LT.0 )THEN
INFO = 3
ELSE IF( K .LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 7
ELSE IF( LDC.LT.MAX( 1, N ) )THEN
INFO = 10
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'DSYRK ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( N.EQ.0 ).OR.
$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
IF( UPPER )THEN
IF( BETA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, J
C( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40, J = 1, N
DO 30, I = 1, J
C( I, J ) = BETA*C( I, J )
30 CONTINUE
40 CONTINUE
END IF
ELSE
IF( BETA.EQ.ZERO )THEN
DO 60, J = 1, N
DO 50, I = J, N
C( I, J ) = ZERO
50 CONTINUE
60 CONTINUE
ELSE
DO 80, J = 1, N
DO 70, I = J, N
C( I, J ) = BETA*C( I, J )
70 CONTINUE
80 CONTINUE
END IF
END IF
RETURN
END IF
*
* Start the operations.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form C := alpha*A*A' + beta*C.
*
IF( UPPER )THEN
DO 130, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 90, I = 1, J
C( I, J ) = ZERO
90 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 100, I = 1, J
C( I, J ) = BETA*C( I, J )
100 CONTINUE
END IF
DO 120, L = 1, K
IF( A( J, L ).NE.ZERO )THEN
TEMP = ALPHA*A( J, L )
DO 110, I = 1, J
C( I, J ) = C( I, J ) + TEMP*A( I, L )
110 CONTINUE
END IF
120 CONTINUE
130 CONTINUE
ELSE
DO 180, J = 1, N
IF( BETA.EQ.ZERO )THEN
DO 140, I = J, N
C( I, J ) = ZERO
140 CONTINUE
ELSE IF( BETA.NE.ONE )THEN
DO 150, I = J, N
C( I, J ) = BETA*C( I, J )
150 CONTINUE
END IF
DO 170, L = 1, K
IF( A( J, L ).NE.ZERO )THEN
TEMP = ALPHA*A( J, L )
DO 160, I = J, N
C( I, J ) = C( I, J ) + TEMP*A( I, L )
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
END IF
ELSE
*
* Form C := alpha*A'*A + beta*C.
*
IF( UPPER )THEN
DO 210, J = 1, N
DO 200, I = 1, J
TEMP = ZERO
DO 190, L = 1, K
TEMP = TEMP + A( L, I )*A( L, J )
190 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
200 CONTINUE
210 CONTINUE
ELSE
DO 240, J = 1, N
DO 230, I = J, N
TEMP = ZERO
DO 220, L = 1, K
TEMP = TEMP + A( L, I )*A( L, J )
220 CONTINUE
IF( BETA.EQ.ZERO )THEN
C( I, J ) = ALPHA*TEMP
ELSE
C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
END IF
230 CONTINUE
240 CONTINUE
END IF
END IF
*
RETURN
*
* End of DSYRK .
*
END
*
************************************************************************
*
SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
$ B, LDB )
* .. Scalar Arguments ..
CHARACTER*1 SIDE, UPLO, TRANSA, DIAG
INTEGER M, N, LDA, LDB
DOUBLE PRECISION ALPHA
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* Purpose
* =======
*
* DTRMM performs one of the matrix-matrix operations
*
* B := alpha*op( A )*B, or B := alpha*B*op( A ),
*
* where alpha is a scalar, B is an m by n matrix, A is a unit, or
* non-unit, upper or lower triangular matrix and op( A ) is one of
*
* op( A ) = A or op( A ) = A'.
*
* Parameters
* ==========
*
* SIDE - CHARACTER*1.
* On entry, SIDE specifies whether op( A ) multiplies B from
* the left or right as follows:
*
* SIDE = 'L' or 'l' B := alpha*op( A )*B.
*
* SIDE = 'R' or 'r' B := alpha*B*op( A ).
*
* Unchanged on exit.
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix A is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANSA - CHARACTER*1.
* On entry, TRANSA specifies the form of op( A ) to be used in
* the matrix multiplication as follows:
*
* TRANSA = 'N' or 'n' op( A ) = A.
*
* TRANSA = 'T' or 't' op( A ) = A'.
*
* TRANSA = 'C' or 'c' op( A ) = A'.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit triangular
* as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of B. M must be at
* least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of B. N must be
* at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION.
* On entry, ALPHA specifies the scalar alpha. When alpha is
* zero then A is not referenced and B need not be set before
* entry.
* Unchanged on exit.
*
* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
* Before entry with UPLO = 'U' or 'u', the leading k by k
* upper triangular part of the array A must contain the upper
* triangular matrix and the strictly lower triangular part of
* A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading k by k
* lower triangular part of the array A must contain the lower
* triangular matrix and the strictly upper triangular part of
* A is not referenced.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced either, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When SIDE = 'L' or 'l' then
* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
* then LDA must be at least max( 1, n ).
* Unchanged on exit.
*
* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
* Before entry, the leading m by n part of the array B must
* contain the matrix B, and on exit is overwritten by the
* transformed matrix.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. LDB must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Local Scalars ..
LOGICAL LSIDE, NOUNIT, UPPER
INTEGER I, INFO, J, K, NROWA
DOUBLE PRECISION TEMP
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
LSIDE = LSAME( SIDE , 'L' )
IF( LSIDE )THEN
NROWA = M
ELSE
NROWA = N
END IF
NOUNIT = LSAME( DIAG , 'N' )
UPPER = LSAME( UPLO , 'U' )
*
INFO = 0
IF( ( .NOT.LSIDE ).AND.
$ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
INFO = 2
ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
INFO = 3
ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
$ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
INFO = 4
ELSE IF( M .LT.0 )THEN
INFO = 5
ELSE IF( N .LT.0 )THEN
INFO = 6
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 9
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
INFO = 11
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'DTRMM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
B( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
RETURN
END IF
*
* Start the operations.
*
IF( LSIDE )THEN
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*A*B.
*
IF( UPPER )THEN
DO 50, J = 1, N
DO 40, K = 1, M
IF( B( K, J ).NE.ZERO )THEN
TEMP = ALPHA*B( K, J )
DO 30, I = 1, K - 1
B( I, J ) = B( I, J ) + TEMP*A( I, K )
30 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP*A( K, K )
B( K, J ) = TEMP
END IF
40 CONTINUE
50 CONTINUE
ELSE
DO 80, J = 1, N
DO 70 K = M, 1, -1
IF( B( K, J ).NE.ZERO )THEN
TEMP = ALPHA*B( K, J )
B( K, J ) = TEMP
IF( NOUNIT )
$ B( K, J ) = B( K, J )*A( K, K )
DO 60, I = K + 1, M
B( I, J ) = B( I, J ) + TEMP*A( I, K )
60 CONTINUE
END IF
70 CONTINUE
80 CONTINUE
END IF
ELSE
*
* Form B := alpha*B*A'.
*
IF( UPPER )THEN
DO 110, J = 1, N
DO 100, I = M, 1, -1
TEMP = B( I, J )
IF( NOUNIT )
$ TEMP = TEMP*A( I, I )
DO 90, K = 1, I - 1
TEMP = TEMP + A( K, I )*B( K, J )
90 CONTINUE
B( I, J ) = ALPHA*TEMP
100 CONTINUE
110 CONTINUE
ELSE
DO 140, J = 1, N
DO 130, I = 1, M
TEMP = B( I, J )
IF( NOUNIT )
$ TEMP = TEMP*A( I, I )
DO 120, K = I + 1, M
TEMP = TEMP + A( K, I )*B( K, J )
120 CONTINUE
B( I, J ) = ALPHA*TEMP
130 CONTINUE
140 CONTINUE
END IF
END IF
ELSE
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*B*A.
*
IF( UPPER )THEN
DO 180, J = N, 1, -1
TEMP = ALPHA
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 150, I = 1, M
B( I, J ) = TEMP*B( I, J )
150 CONTINUE
DO 170, K = 1, J - 1
IF( A( K, J ).NE.ZERO )THEN
TEMP = ALPHA*A( K, J )
DO 160, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
ELSE
DO 220, J = 1, N
TEMP = ALPHA
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 190, I = 1, M
B( I, J ) = TEMP*B( I, J )
190 CONTINUE
DO 210, K = J + 1, N
IF( A( K, J ).NE.ZERO )THEN
TEMP = ALPHA*A( K, J )
DO 200, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
200 CONTINUE
END IF
210 CONTINUE
220 CONTINUE
END IF
ELSE
*
* Form B := alpha*B*A'.
*
IF( UPPER )THEN
DO 260, K = 1, N
DO 240, J = 1, K - 1
IF( A( J, K ).NE.ZERO )THEN
TEMP = ALPHA*A( J, K )
DO 230, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
230 CONTINUE
END IF
240 CONTINUE
TEMP = ALPHA
IF( NOUNIT )
$ TEMP = TEMP*A( K, K )
IF( TEMP.NE.ONE )THEN
DO 250, I = 1, M
B( I, K ) = TEMP*B( I, K )
250 CONTINUE
END IF
260 CONTINUE
ELSE
DO 300, K = N, 1, -1
DO 280, J = K + 1, N
IF( A( J, K ).NE.ZERO )THEN
TEMP = ALPHA*A( J, K )
DO 270, I = 1, M
B( I, J ) = B( I, J ) + TEMP*B( I, K )
270 CONTINUE
END IF
280 CONTINUE
TEMP = ALPHA
IF( NOUNIT )
$ TEMP = TEMP*A( K, K )
IF( TEMP.NE.ONE )THEN
DO 290, I = 1, M
B( I, K ) = TEMP*B( I, K )
290 CONTINUE
END IF
300 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTRMM .
*
END
*
************************************************************************
*
SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
$ B, LDB )
* .. Scalar Arguments ..
CHARACTER*1 SIDE, UPLO, TRANSA, DIAG
INTEGER M, N, LDA, LDB
DOUBLE PRECISION ALPHA
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* Purpose
* =======
*
* DTRSM solves one of the matrix equations
*
* op( A )*X = alpha*B, or X*op( A ) = alpha*B,
*
* where alpha is a scalar, X and B are m by n matrices, A is a unit, or
* non-unit, upper or lower triangular matrix and op( A ) is one of
*
* op( A ) = A or op( A ) = A'.
*
* The matrix X is overwritten on B.
*
* Parameters
* ==========
*
* SIDE - CHARACTER*1.
* On entry, SIDE specifies whether op( A ) appears on the left
* or right of X as follows:
*
* SIDE = 'L' or 'l' op( A )*X = alpha*B.
*
* SIDE = 'R' or 'r' X*op( A ) = alpha*B.
*
* Unchanged on exit.
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix A is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANSA - CHARACTER*1.
* On entry, TRANSA specifies the form of op( A ) to be used in
* the matrix multiplication as follows:
*
* TRANSA = 'N' or 'n' op( A ) = A.
*
* TRANSA = 'T' or 't' op( A ) = A'.
*
* TRANSA = 'C' or 'c' op( A ) = A'.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit triangular
* as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of B. M must be at
* least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of B. N must be
* at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION.
* On entry, ALPHA specifies the scalar alpha. When alpha is
* zero then A is not referenced and B need not be set before
* entry.
* Unchanged on exit.
*
* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
* Before entry with UPLO = 'U' or 'u', the leading k by k
* upper triangular part of the array A must contain the upper
* triangular matrix and the strictly lower triangular part of
* A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading k by k
* lower triangular part of the array A must contain the lower
* triangular matrix and the strictly upper triangular part of
* A is not referenced.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced either, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. When SIDE = 'L' or 'l' then
* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
* then LDA must be at least max( 1, n ).
* Unchanged on exit.
*
* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
* Before entry, the leading m by n part of the array B must
* contain the right-hand side matrix B, and on exit is
* overwritten by the solution matrix X.
*
* LDB - INTEGER.
* On entry, LDB specifies the first dimension of B as declared
* in the calling (sub) program. LDB must be at least
* max( 1, m ).
* Unchanged on exit.
*
*
* Level 3 Blas routine.
*
*
* -- Written on 8-February-1989.
* Jack Dongarra, Argonne National Laboratory.
* Iain Duff, AERE Harwell.
* Jeremy Du Croz, Numerical Algorithms Group Ltd.
* Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Local Scalars ..
LOGICAL LSIDE, NOUNIT, UPPER
INTEGER I, INFO, J, K, NROWA
DOUBLE PRECISION TEMP
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
LSIDE = LSAME( SIDE , 'L' )
IF( LSIDE )THEN
NROWA = M
ELSE
NROWA = N
END IF
NOUNIT = LSAME( DIAG , 'N' )
UPPER = LSAME( UPLO , 'U' )
*
INFO = 0
IF( ( .NOT.LSIDE ).AND.
$ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN
INFO = 1
ELSE IF( ( .NOT.UPPER ).AND.
$ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN
INFO = 2
ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
$ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN
INFO = 3
ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.
$ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN
INFO = 4
ELSE IF( M .LT.0 )THEN
INFO = 5
ELSE IF( N .LT.0 )THEN
INFO = 6
ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
INFO = 9
ELSE IF( LDB.LT.MAX( 1, M ) )THEN
INFO = 11
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'DTRSM ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
* And when alpha.eq.zero.
*
IF( ALPHA.EQ.ZERO )THEN
DO 20, J = 1, N
DO 10, I = 1, M
B( I, J ) = ZERO
10 CONTINUE
20 CONTINUE
RETURN
END IF
*
* Start the operations.
*
IF( LSIDE )THEN
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*inv( A )*B.
*
IF( UPPER )THEN
DO 60, J = 1, N
IF( ALPHA.NE.ONE )THEN
DO 30, I = 1, M
B( I, J ) = ALPHA*B( I, J )
30 CONTINUE
END IF
DO 50, K = M, 1, -1
IF( B( K, J ).NE.ZERO )THEN
IF( NOUNIT )
$ B( K, J ) = B( K, J )/A( K, K )
DO 40, I = 1, K - 1
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
40 CONTINUE
END IF
50 CONTINUE
60 CONTINUE
ELSE
DO 100, J = 1, N
IF( ALPHA.NE.ONE )THEN
DO 70, I = 1, M
B( I, J ) = ALPHA*B( I, J )
70 CONTINUE
END IF
DO 90 K = 1, M
IF( B( K, J ).NE.ZERO )THEN
IF( NOUNIT )
$ B( K, J ) = B( K, J )/A( K, K )
DO 80, I = K + 1, M
B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
80 CONTINUE
END IF
90 CONTINUE
100 CONTINUE
END IF
ELSE
*
* Form B := alpha*inv( A' )*B.
*
IF( UPPER )THEN
DO 130, J = 1, N
DO 120, I = 1, M
TEMP = ALPHA*B( I, J )
DO 110, K = 1, I - 1
TEMP = TEMP - A( K, I )*B( K, J )
110 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( I, I )
B( I, J ) = TEMP
120 CONTINUE
130 CONTINUE
ELSE
DO 160, J = 1, N
DO 150, I = M, 1, -1
TEMP = ALPHA*B( I, J )
DO 140, K = I + 1, M
TEMP = TEMP - A( K, I )*B( K, J )
140 CONTINUE
IF( NOUNIT )
$ TEMP = TEMP/A( I, I )
B( I, J ) = TEMP
150 CONTINUE
160 CONTINUE
END IF
END IF
ELSE
IF( LSAME( TRANSA, 'N' ) )THEN
*
* Form B := alpha*B*inv( A ).
*
IF( UPPER )THEN
DO 210, J = 1, N
IF( ALPHA.NE.ONE )THEN
DO 170, I = 1, M
B( I, J ) = ALPHA*B( I, J )
170 CONTINUE
END IF
DO 190, K = 1, J - 1
IF( A( K, J ).NE.ZERO )THEN
DO 180, I = 1, M
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
180 CONTINUE
END IF
190 CONTINUE
IF( NOUNIT )THEN
TEMP = ONE/A( J, J )
DO 200, I = 1, M
B( I, J ) = TEMP*B( I, J )
200 CONTINUE
END IF
210 CONTINUE
ELSE
DO 260, J = N, 1, -1
IF( ALPHA.NE.ONE )THEN
DO 220, I = 1, M
B( I, J ) = ALPHA*B( I, J )
220 CONTINUE
END IF
DO 240, K = J + 1, N
IF( A( K, J ).NE.ZERO )THEN
DO 230, I = 1, M
B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
230 CONTINUE
END IF
240 CONTINUE
IF( NOUNIT )THEN
TEMP = ONE/A( J, J )
DO 250, I = 1, M
B( I, J ) = TEMP*B( I, J )
250 CONTINUE
END IF
260 CONTINUE
END IF
ELSE
*
* Form B := alpha*B*inv( A' ).
*
IF( UPPER )THEN
DO 310, K = N, 1, -1
IF( NOUNIT )THEN
TEMP = ONE/A( K, K )
DO 270, I = 1, M
B( I, K ) = TEMP*B( I, K )
270 CONTINUE
END IF
DO 290, J = 1, K - 1
IF( A( J, K ).NE.ZERO )THEN
TEMP = A( J, K )
DO 280, I = 1, M
B( I, J ) = B( I, J ) - TEMP*B( I, K )
280 CONTINUE
END IF
290 CONTINUE
IF( ALPHA.NE.ONE )THEN
DO 300, I = 1, M
B( I, K ) = ALPHA*B( I, K )
300 CONTINUE
END IF
310 CONTINUE
ELSE
DO 360, K = 1, N
IF( NOUNIT )THEN
TEMP = ONE/A( K, K )
DO 320, I = 1, M
B( I, K ) = TEMP*B( I, K )
320 CONTINUE
END IF
DO 340, J = K + 1, N
IF( A( J, K ).NE.ZERO )THEN
TEMP = A( J, K )
DO 330, I = 1, M
B( I, J ) = B( I, J ) - TEMP*B( I, K )
330 CONTINUE
END IF
340 CONTINUE
IF( ALPHA.NE.ONE )THEN
DO 350, I = 1, M
B( I, K ) = ALPHA*B( I, K )
350 CONTINUE
END IF
360 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTRSM .
*
END
*
************************************************************************
*
* File of the DOUBLE PRECISION Level-2 BLAS.
* ===========================================
*
* SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
* $ BETA, Y, INCY )
*
* SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
*
* SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
*
* SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
*
* SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
*
* SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
*
* SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
*
* SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
*
* SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
*
* SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA )
*
* SUBROUTINE DSPR ( UPLO, N, ALPHA, X, INCX, AP )
*
* SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
*
* SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
*
* See:
*
* Dongarra J. J., Du Croz J. J., Hammarling S. and Hanson R. J..
* An extended set of Fortran Basic Linear Algebra Subprograms.
*
* Technical Memoranda Nos. 41 (revision 3) and 81, Mathematics
* and Computer Science Division, Argonne National Laboratory,
* 9700 South Cass Avenue, Argonne, Illinois 60439, US.
*
* Or
*
* NAG Technical Reports TR3/87 and TR4/87, Numerical Algorithms
* Group Ltd., NAG Central Office, 256 Banbury Road, Oxford
* OX2 7DE, UK, and Numerical Algorithms Group Inc., 1101 31st
* Street, Suite 100, Downers Grove, Illinois 60515-1263, USA.
*
************************************************************************
*
SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
$ BETA, Y, INCY )
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA, BETA
INTEGER INCX, INCY, LDA, M, N
CHARACTER*1 TRANS
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
* ..
*
* Purpose
* =======
*
* DGEMV performs one of the matrix-vector operations
*
* y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
*
* where alpha and beta are scalars, x and y are vectors and A is an
* m by n matrix.
*
* Parameters
* ==========
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
*
* TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
*
* TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
*
* Unchanged on exit.
*
* M - INTEGER.
* On entry, M specifies the number of rows of the matrix A.
* M must be at least zero.
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the number of columns of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* ALPHA - DOUBLE PRECISION.
* On entry, ALPHA specifies the scalar alpha.
* Unchanged on exit.
*
* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
* Before entry, the leading m by n part of the array A must
* contain the matrix of coefficients.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, m ).
* Unchanged on exit.
*
* X - DOUBLE PRECISION array of DIMENSION at least
* ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
* and at least
* ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
* Before entry, the incremented array X must contain the
* vector x.
* Unchanged on exit.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
* BETA - DOUBLE PRECISION.
* On entry, BETA specifies the scalar beta. When BETA is
* supplied as zero then Y need not be set on input.
* Unchanged on exit.
*
* Y - DOUBLE PRECISION array of DIMENSION at least
* ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
* and at least
* ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
* Before entry with BETA non-zero, the incremented array Y
* must contain the vector y. On exit, Y is overwritten by the
* updated vector y.
*
* INCY - INTEGER.
* On entry, INCY specifies the increment for the elements of
* Y. INCY must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
DOUBLE PRECISION ONE , ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 1
ELSE IF( M.LT.0 )THEN
INFO = 2
ELSE IF( N.LT.0 )THEN
INFO = 3
ELSE IF( LDA.LT.MAX( 1, M ) )THEN
INFO = 6
ELSE IF( INCX.EQ.0 )THEN
INFO = 8
ELSE IF( INCY.EQ.0 )THEN
INFO = 11
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'DGEMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
$ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
$ RETURN
*
* Set LENX and LENY, the lengths of the vectors x and y, and set
* up the start points in X and Y.
*
IF( LSAME( TRANS, 'N' ) )THEN
LENX = N
LENY = M
ELSE
LENX = M
LENY = N
END IF
IF( INCX.GT.0 )THEN
KX = 1
ELSE
KX = 1 - ( LENX - 1 )*INCX
END IF
IF( INCY.GT.0 )THEN
KY = 1
ELSE
KY = 1 - ( LENY - 1 )*INCY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
* First form y := beta*y.
*
IF( BETA.NE.ONE )THEN
IF( INCY.EQ.1 )THEN
IF( BETA.EQ.ZERO )THEN
DO 10, I = 1, LENY
Y( I ) = ZERO
10 CONTINUE
ELSE
DO 20, I = 1, LENY
Y( I ) = BETA*Y( I )
20 CONTINUE
END IF
ELSE
IY = KY
IF( BETA.EQ.ZERO )THEN
DO 30, I = 1, LENY
Y( IY ) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40, I = 1, LENY
Y( IY ) = BETA*Y( IY )
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF( ALPHA.EQ.ZERO )
$ RETURN
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form y := alpha*A*x + y.
*
JX = KX
IF( INCY.EQ.1 )THEN
DO 60, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*X( JX )
DO 50, I = 1, M
Y( I ) = Y( I ) + TEMP*A( I, J )
50 CONTINUE
END IF
JX = JX + INCX
60 CONTINUE
ELSE
DO 80, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = ALPHA*X( JX )
IY = KY
DO 70, I = 1, M
Y( IY ) = Y( IY ) + TEMP*A( I, J )
IY = IY + INCY
70 CONTINUE
END IF
JX = JX + INCX
80 CONTINUE
END IF
ELSE
*
* Form y := alpha*A'*x + y.
*
JY = KY
IF( INCX.EQ.1 )THEN
DO 100, J = 1, N
TEMP = ZERO
DO 90, I = 1, M
TEMP = TEMP + A( I, J )*X( I )
90 CONTINUE
Y( JY ) = Y( JY ) + ALPHA*TEMP
JY = JY + INCY
100 CONTINUE
ELSE
DO 120, J = 1, N
TEMP = ZERO
IX = KX
DO 110, I = 1, M
TEMP = TEMP + A( I, J )*X( IX )
IX = IX + INCX
110 CONTINUE
Y( JY ) = Y( JY ) + ALPHA*TEMP
JY = JY + INCY
120 CONTINUE
END IF
END IF
*
RETURN
*
* End of DGEMV .
*
END
c
subroutine dscal(n,da,dx,incx)
c
c scales a vector by a constant.
c uses unrolled loops for increment equal to one.
c jack dongarra, linpack, 3/11/78.
c modified to correct problem with negative increment, 8/21/90.
c
double precision da,dx(1)
integer i,incx,ix,m,mp1,n
c
if(n.le.0)return
if(incx.eq.1)go to 20
c
c code for increment not equal to 1
c
ix = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
do 10 i = 1,n
dx(ix) = da*dx(ix)
ix = ix + incx
10 continue
return
c
c code for increment equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dx(i) = da*dx(i)
30 continue
if( n .lt. 5 ) return
40 mp1 = m + 1
do 50 i = mp1,n,5
dx(i) = da*dx(i)
dx(i + 1) = da*dx(i + 1)
dx(i + 2) = da*dx(i + 2)
dx(i + 3) = da*dx(i + 3)
dx(i + 4) = da*dx(i + 4)
50 continue
return
end
c
double precision function ddot(n,dx,incx,dy,incy)
c
c forms the dot product of two vectors.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
ddot = 0.0d0
dtemp = 0.0d0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
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
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
* dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
50 continue
60 ddot = dtemp
return
end
*
************************************************************************
*
SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
* .. Scalar Arguments ..
INTEGER INCX, LDA, N
CHARACTER*1 DIAG, TRANS, UPLO
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), X( * )
* ..
*
* Purpose
* =======
*
* DTRMV performs one of the matrix-vector operations
*
* x := A*x, or x := A'*x,
*
* where x is an n element vector and A is an n by n unit, or non-unit,
* upper or lower triangular matrix.
*
* Parameters
* ==========
*
* UPLO - CHARACTER*1.
* On entry, UPLO specifies whether the matrix is an upper or
* lower triangular matrix as follows:
*
* UPLO = 'U' or 'u' A is an upper triangular matrix.
*
* UPLO = 'L' or 'l' A is a lower triangular matrix.
*
* Unchanged on exit.
*
* TRANS - CHARACTER*1.
* On entry, TRANS specifies the operation to be performed as
* follows:
*
* TRANS = 'N' or 'n' x := A*x.
*
* TRANS = 'T' or 't' x := A'*x.
*
* TRANS = 'C' or 'c' x := A'*x.
*
* Unchanged on exit.
*
* DIAG - CHARACTER*1.
* On entry, DIAG specifies whether or not A is unit
* triangular as follows:
*
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
*
* DIAG = 'N' or 'n' A is not assumed to be unit
* triangular.
*
* Unchanged on exit.
*
* N - INTEGER.
* On entry, N specifies the order of the matrix A.
* N must be at least zero.
* Unchanged on exit.
*
* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
* Before entry with UPLO = 'U' or 'u', the leading n by n
* upper triangular part of the array A must contain the upper
* triangular matrix and the strictly lower triangular part of
* A is not referenced.
* Before entry with UPLO = 'L' or 'l', the leading n by n
* lower triangular part of the array A must contain the lower
* triangular matrix and the strictly upper triangular part of
* A is not referenced.
* Note that when DIAG = 'U' or 'u', the diagonal elements of
* A are not referenced either, but are assumed to be unity.
* Unchanged on exit.
*
* LDA - INTEGER.
* On entry, LDA specifies the first dimension of A as declared
* in the calling (sub) program. LDA must be at least
* max( 1, n ).
* Unchanged on exit.
*
* X - DOUBLE PRECISION array of dimension at least
* ( 1 + ( n - 1 )*abs( INCX ) ).
* Before entry, the incremented array X must contain the n
* element vector x. On exit, X is overwritten with the
* tranformed vector x.
*
* INCX - INTEGER.
* On entry, INCX specifies the increment for the elements of
* X. INCX must not be zero.
* Unchanged on exit.
*
*
* Level 2 Blas routine.
*
* -- Written on 22-October-1986.
* Jack Dongarra, Argonne National Lab.
* Jeremy Du Croz, Nag Central Office.
* Sven Hammarling, Nag Central Office.
* Richard Hanson, Sandia National Labs.
*
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D+0 )
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I, INFO, IX, J, JX, KX
LOGICAL NOUNIT
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL XERBLA
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF ( .NOT.LSAME( UPLO , 'U' ).AND.
$ .NOT.LSAME( UPLO , 'L' ) )THEN
INFO = 1
ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
$ .NOT.LSAME( TRANS, 'T' ).AND.
$ .NOT.LSAME( TRANS, 'C' ) )THEN
INFO = 2
ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
$ .NOT.LSAME( DIAG , 'N' ) )THEN
INFO = 3
ELSE IF( N.LT.0 )THEN
INFO = 4
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
INFO = 6
ELSE IF( INCX.EQ.0 )THEN
INFO = 8
END IF
IF( INFO.NE.0 )THEN
CALL XERBLA( 'DTRMV ', INFO )
RETURN
END IF
*
* Quick return if possible.
*
IF( N.EQ.0 )
$ RETURN
*
NOUNIT = LSAME( DIAG, 'N' )
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF( INCX.LE.0 )THEN
KX = 1 - ( N - 1 )*INCX
ELSE IF( INCX.NE.1 )THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF( LSAME( TRANS, 'N' ) )THEN
*
* Form x := A*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
IF( INCX.EQ.1 )THEN
DO 20, J = 1, N
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
DO 10, I = 1, J - 1
X( I ) = X( I ) + TEMP*A( I, J )
10 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*A( J, J )
END IF
20 CONTINUE
ELSE
JX = KX
DO 40, J = 1, N
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
DO 30, I = 1, J - 1
X( IX ) = X( IX ) + TEMP*A( I, J )
IX = IX + INCX
30 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*A( J, J )
END IF
JX = JX + INCX
40 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 60, J = N, 1, -1
IF( X( J ).NE.ZERO )THEN
TEMP = X( J )
DO 50, I = N, J + 1, -1
X( I ) = X( I ) + TEMP*A( I, J )
50 CONTINUE
IF( NOUNIT )
$ X( J ) = X( J )*A( J, J )
END IF
60 CONTINUE
ELSE
KX = KX + ( N - 1 )*INCX
JX = KX
DO 80, J = N, 1, -1
IF( X( JX ).NE.ZERO )THEN
TEMP = X( JX )
IX = KX
DO 70, I = N, J + 1, -1
X( IX ) = X( IX ) + TEMP*A( I, J )
IX = IX - INCX
70 CONTINUE
IF( NOUNIT )
$ X( JX ) = X( JX )*A( J, J )
END IF
JX = JX - INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := A'*x.
*
IF( LSAME( UPLO, 'U' ) )THEN
IF( INCX.EQ.1 )THEN
DO 100, J = N, 1, -1
TEMP = X( J )
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 90, I = J - 1, 1, -1
TEMP = TEMP + A( I, J )*X( I )
90 CONTINUE
X( J ) = TEMP
100 CONTINUE
ELSE
JX = KX + ( N - 1 )*INCX
DO 120, J = N, 1, -1
TEMP = X( JX )
IX = JX
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 110, I = J - 1, 1, -1
IX = IX - INCX
TEMP = TEMP + A( I, J )*X( IX )
110 CONTINUE
X( JX ) = TEMP
JX = JX - INCX
120 CONTINUE
END IF
ELSE
IF( INCX.EQ.1 )THEN
DO 140, J = 1, N
TEMP = X( J )
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 130, I = J + 1, N
TEMP = TEMP + A( I, J )*X( I )
130 CONTINUE
X( J ) = TEMP
140 CONTINUE
ELSE
JX = KX
DO 160, J = 1, N
TEMP = X( JX )
IX = JX
IF( NOUNIT )
$ TEMP = TEMP*A( J, J )
DO 150, I = J + 1, N
IX = IX + INCX
TEMP = TEMP + A( I, J )*X( IX )
150 CONTINUE
X( JX ) = TEMP
JX = JX + INCX
160 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTRMV .
*
END