The R Project SVN
R-packages
(root)/trunk/Matrix/src/zherk.f @ 1908 – Rev
Blame | Last modification | View Log | Download | RSS feed
C-- straight copy of http://www.netlib.org/blas/zherk.fCSUBROUTINE ZHERK( UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC )* .. Scalar Arguments ..CHARACTER TRANS, UPLOINTEGER K, LDA, LDC, NDOUBLE PRECISION ALPHA, BETA* ..* .. Array Arguments ..COMPLEX*16 A( LDA, * ), C( LDC, * )* ..** Purpose* =======** ZHERK performs one of the hermitian rank k operations** C := alpha*A*conjg( A' ) + beta*C,** or** C := alpha*conjg( A' )*A + beta*C,** where alpha and beta are real scalars, C is an n by n hermitian* 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*conjg( A' ) + beta*C.** TRANS = 'C' or 'c' C := alpha*conjg( 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 = '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 - COMPLEX*16 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 - COMPLEX*16 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 hermitian 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 hermitian 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.* Note that the imaginary parts of the diagonal elements need* not be set, they are assumed to be zero, and on exit they* are set to zero.** 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.** -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1.* Ed Anderson, Cray Research Inc.*** .. External Functions ..LOGICAL LSAMEEXTERNAL LSAME* ..* .. External Subroutines ..EXTERNAL XERBLA* ..* .. Intrinsic Functions ..INTRINSIC DBLE, DCMPLX, DCONJG, MAX* ..* .. Local Scalars ..LOGICAL UPPERINTEGER I, INFO, J, L, NROWADOUBLE PRECISION RTEMPCOMPLEX*16 TEMP* ..* .. Parameters ..DOUBLE PRECISION ONE, ZEROPARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )* ..* .. Executable Statements ..** Test the input parameters.*IF( LSAME( TRANS, 'N' ) ) THENNROWA = NELSENROWA = KEND IFUPPER = LSAME( UPLO, 'U' )*INFO = 0IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THENINFO = 1ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND.$ ( .NOT.LSAME( TRANS, 'C' ) ) ) THENINFO = 2ELSE IF( N.LT.0 ) THENINFO = 3ELSE IF( K.LT.0 ) THENINFO = 4ELSE IF( LDA.LT.MAX( 1, NROWA ) ) THENINFO = 7ELSE IF( LDC.LT.MAX( 1, N ) ) THENINFO = 10END IFIF( INFO.NE.0 ) THENCALL XERBLA( 'ZHERK ', INFO )RETURNEND 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 ) THENIF( UPPER ) THENIF( BETA.EQ.ZERO ) THENDO 20 J = 1, NDO 10 I = 1, JC( I, J ) = ZERO10 CONTINUE20 CONTINUEELSEDO 40 J = 1, NDO 30 I = 1, J - 1C( I, J ) = BETA*C( I, J )30 CONTINUEC( J, J ) = BETA*DBLE( C( J, J ) )40 CONTINUEEND IFELSEIF( BETA.EQ.ZERO ) THENDO 60 J = 1, NDO 50 I = J, NC( I, J ) = ZERO50 CONTINUE60 CONTINUEELSEDO 80 J = 1, NC( J, J ) = BETA*DBLE( C( J, J ) )DO 70 I = J + 1, NC( I, J ) = BETA*C( I, J )70 CONTINUE80 CONTINUEEND IFEND IFRETURNEND IF** Start the operations.*IF( LSAME( TRANS, 'N' ) ) THEN** Form C := alpha*A*conjg( A' ) + beta*C.*IF( UPPER ) THENDO 130 J = 1, NIF( BETA.EQ.ZERO ) THENDO 90 I = 1, JC( I, J ) = ZERO90 CONTINUEELSE IF( BETA.NE.ONE ) THENDO 100 I = 1, J - 1C( I, J ) = BETA*C( I, J )100 CONTINUEC( J, J ) = BETA*DBLE( C( J, J ) )ELSEC( J, J ) = DBLE( C( J, J ) )END IFDO 120 L = 1, KIF( A( J, L ).NE.DCMPLX( ZERO ) ) THENTEMP = ALPHA*DCONJG( A( J, L ) )DO 110 I = 1, J - 1C( I, J ) = C( I, J ) + TEMP*A( I, L )110 CONTINUEC( J, J ) = DBLE( C( J, J ) ) +$ DBLE( TEMP*A( I, L ) )END IF120 CONTINUE130 CONTINUEELSEDO 180 J = 1, NIF( BETA.EQ.ZERO ) THENDO 140 I = J, NC( I, J ) = ZERO140 CONTINUEELSE IF( BETA.NE.ONE ) THENC( J, J ) = BETA*DBLE( C( J, J ) )DO 150 I = J + 1, NC( I, J ) = BETA*C( I, J )150 CONTINUEELSEC( J, J ) = DBLE( C( J, J ) )END IFDO 170 L = 1, KIF( A( J, L ).NE.DCMPLX( ZERO ) ) THENTEMP = ALPHA*DCONJG( A( J, L ) )C( J, J ) = DBLE( C( J, J ) ) +$ DBLE( TEMP*A( J, L ) )DO 160 I = J + 1, NC( I, J ) = C( I, J ) + TEMP*A( I, L )160 CONTINUEEND IF170 CONTINUE180 CONTINUEEND IFELSE** Form C := alpha*conjg( A' )*A + beta*C.*IF( UPPER ) THENDO 220 J = 1, NDO 200 I = 1, J - 1TEMP = ZERODO 190 L = 1, KTEMP = TEMP + DCONJG( A( L, I ) )*A( L, J )190 CONTINUEIF( BETA.EQ.ZERO ) THENC( I, J ) = ALPHA*TEMPELSEC( I, J ) = ALPHA*TEMP + BETA*C( I, J )END IF200 CONTINUERTEMP = ZERODO 210 L = 1, KRTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J )210 CONTINUEIF( BETA.EQ.ZERO ) THENC( J, J ) = ALPHA*RTEMPELSEC( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) )END IF220 CONTINUEELSEDO 260 J = 1, NRTEMP = ZERODO 230 L = 1, KRTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J )230 CONTINUEIF( BETA.EQ.ZERO ) THENC( J, J ) = ALPHA*RTEMPELSEC( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) )END IFDO 250 I = J + 1, NTEMP = ZERODO 240 L = 1, KTEMP = TEMP + DCONJG( A( L, I ) )*A( L, J )240 CONTINUEIF( BETA.EQ.ZERO ) THENC( I, J ) = ALPHA*TEMPELSEC( I, J ) = ALPHA*TEMP + BETA*C( I, J )END IF250 CONTINUE260 CONTINUEEND IFEND IF*RETURN** End of ZHERK .*END