Mercurial > octave-nkf
changeset 8692:54227442f7ed
add missing BLAS sources
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 06 Feb 2009 07:30:55 +0100 |
| parents | 7838271ee25c |
| children | e5ffb52c9c61 |
| files | libcruft/ChangeLog libcruft/blas/Makefile.in libcruft/blas/chemm.f libcruft/blas/dsymm.f libcruft/blas/ssymm.f libcruft/blas/zhemm.f |
| diffstat | 6 files changed, 1193 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/ChangeLog Thu Feb 05 19:58:04 2009 -0500 +++ b/libcruft/ChangeLog Fri Feb 06 07:30:55 2009 +0100 @@ -1,3 +1,8 @@ +2009-02-06 Jaroslav Hajek <highegg@gmail.com> + + * blas/ssymm.f, blas/dsymm.f, blas/chemm.f, blas/zhemm.f: New sources. + * blas/Makefile.in: Include them. + 2009-01-28 John W. Eaton <jwe@octave.org> * Makefile.in (LIBRARIES, install, uninstall): use SHLLIBPRE and
--- a/libcruft/blas/Makefile.in Thu Feb 05 19:58:04 2009 -0500 +++ b/libcruft/blas/Makefile.in Fri Feb 06 07:30:55 2009 +0100 @@ -27,16 +27,16 @@ EXTERNAL_DISTFILES = $(DISTFILES) FSRC = dasum.f daxpy.f dcabs1.f dcopy.f ddot.f dgemm.f dgemv.f \ - dger.f dmach.f dnrm2.f drot.f dscal.f dswap.f dsymv.f dsyr.f \ + dger.f dmach.f dnrm2.f drot.f dscal.f dswap.f dsymv.f dsyr.f dsymm.f \ dsyr2.f dsyr2k.f dsyrk.f dtbsv.f dtrmm.f dtrmv.f dtrsm.f dtrsv.f \ dzasum.f dznrm2.f icamax.f idamax.f isamax.f izamax.f lsame.f sdot.f \ sgemm.f sgemv.f sscal.f ssyrk.f strsm.f zaxpy.f zcopy.f zdotc.f \ - zdotu.f zdrot.f zdscal.f zgemm.f zgemv.f zgerc.f zgeru.f zhemv.f \ + zdotu.f zdrot.f zdscal.f zgemm.f zgemv.f zgerc.f zgeru.f zhemv.f zhemm.f \ zher.f zher2.f zher2k.f zherk.f zscal.f zswap.f zsyrk.f ztbsv.f ztrmm.f \ ztrmv.f ztrsm.f ztrsv.f sasum.f saxpy.f scabs1.f scopy.f \ - sger.f smach.f snrm2.f srot.f sswap.f ssymv.f ssyr.f \ + sger.f smach.f snrm2.f srot.f sswap.f ssymv.f ssyr.f ssymm.f \ ssyr2.f ssyr2k.f stbsv.f strmm.f strmv.f strsv.f \ - scasum.f scnrm2.f caxpy.f ccopy.f cdotc.f cdotu.f \ + scasum.f scnrm2.f caxpy.f ccopy.f cdotc.f cdotu.f chemm.f \ csrot.f csscal.f cgemm.f cgemv.f cgerc.f cgeru.f chemv.f cher.f \ cher2.f cher2k.f cherk.f cscal.f cswap.f csyrk.f ctbsv.f ctrmm.f ctrmv.f \ ctrsm.f ctrsv.f
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/chemm.f Fri Feb 06 07:30:55 2009 +0100 @@ -0,0 +1,298 @@ + SUBROUTINE CHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* .. Scalar Arguments .. + COMPLEX ALPHA,BETA + INTEGER LDA,LDB,LDC,M,N + CHARACTER SIDE,UPLO +* .. +* .. Array Arguments .. + COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* Purpose +* ======= +* +* CHEMM performs one of the matrix-matrix operations +* +* C := alpha*A*B + beta*C, +* +* or +* +* C := alpha*B*A + beta*C, +* +* where alpha and beta are scalars, A is an hermitian matrix and B and +* C are m by n matrices. +* +* Arguments +* ========== +* +* SIDE - CHARACTER*1. +* On entry, SIDE specifies whether the hermitian matrix A +* appears on the left or right in the operation as follows: +* +* SIDE = 'L' or 'l' C := alpha*A*B + beta*C, +* +* SIDE = 'R' or 'r' C := alpha*B*A + beta*C, +* +* Unchanged on exit. +* +* UPLO - CHARACTER*1. +* On entry, UPLO specifies whether the upper or lower +* triangular part of the hermitian matrix A is to be +* referenced as follows: +* +* UPLO = 'U' or 'u' Only the upper triangular part of the +* hermitian matrix is to be referenced. +* +* UPLO = 'L' or 'l' Only the lower triangular part of the +* hermitian matrix is to be referenced. +* +* Unchanged on exit. +* +* M - INTEGER. +* On entry, M specifies the number of rows 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 C. +* N must be at least zero. +* Unchanged on exit. +* +* ALPHA - COMPLEX . +* On entry, ALPHA specifies the scalar alpha. +* Unchanged on exit. +* +* A - COMPLEX array of DIMENSION ( LDA, ka ), where ka is +* m when SIDE = 'L' or 'l' and is n otherwise. +* Before entry with SIDE = 'L' or 'l', the m by m part of +* the array A must contain the hermitian matrix, such that +* when UPLO = 'U' or 'u', the leading m by m upper triangular +* part of the array A must contain the upper triangular part +* of the hermitian matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading m by m lower triangular part of the array A +* must contain the lower triangular part of the hermitian +* matrix and the strictly upper triangular part of A is not +* referenced. +* Before entry with SIDE = 'R' or 'r', the n by n part of +* the array A must contain the hermitian matrix, such that +* when UPLO = 'U' or 'u', the leading n by n upper triangular +* part of the array A must contain the upper triangular part +* of the hermitian matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading n by n lower triangular part of the array A +* must contain the lower triangular part of the hermitian +* matrix and the strictly upper triangular part of A is not +* referenced. +* Note that the imaginary parts of the diagonal elements need +* not be set, they are assumed to be zero. +* 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 ), otherwise LDA must be at +* least max( 1, n ). +* Unchanged on exit. +* +* B - COMPLEX array of DIMENSION ( LDB, n ). +* Before entry, the leading m by n 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. LDB must be at least +* max( 1, m ). +* Unchanged on exit. +* +* BETA - COMPLEX . +* 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 - COMPLEX 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 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, 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 CONJG,MAX,REAL +* .. +* .. Local Scalars .. + COMPLEX TEMP1,TEMP2 + INTEGER I,INFO,J,K,NROWA + LOGICAL UPPER +* .. +* .. Parameters .. + COMPLEX ONE + PARAMETER (ONE= (1.0E+0,0.0E+0)) + COMPLEX ZERO + PARAMETER (ZERO= (0.0E+0,0.0E+0)) +* .. +* +* Set NROWA as the number of rows of A. +* + IF (LSAME(SIDE,'L')) THEN + NROWA = M + ELSE + NROWA = N + END IF + UPPER = LSAME(UPLO,'U') +* +* Test the input parameters. +* + INFO = 0 + IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN + INFO = 1 + ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN + INFO = 2 + ELSE IF (M.LT.0) THEN + INFO = 3 + ELSE IF (N.LT.0) THEN + INFO = 4 + ELSE IF (LDA.LT.MAX(1,NROWA)) THEN + INFO = 7 + ELSE IF (LDB.LT.MAX(1,M)) THEN + INFO = 9 + ELSE IF (LDC.LT.MAX(1,M)) THEN + INFO = 12 + END IF + IF (INFO.NE.0) THEN + CALL XERBLA('CHEMM ',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 +* +* And when 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 (LSAME(SIDE,'L')) THEN +* +* Form C := alpha*A*B + beta*C. +* + IF (UPPER) THEN + DO 70 J = 1,N + DO 60 I = 1,M + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 50 K = 1,I - 1 + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*CONJG(A(K,I)) + 50 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*REAL(A(I,I)) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*REAL(A(I,I)) + + + ALPHA*TEMP2 + END IF + 60 CONTINUE + 70 CONTINUE + ELSE + DO 100 J = 1,N + DO 90 I = M,1,-1 + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 80 K = I + 1,M + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*CONJG(A(K,I)) + 80 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*REAL(A(I,I)) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*REAL(A(I,I)) + + + ALPHA*TEMP2 + END IF + 90 CONTINUE + 100 CONTINUE + END IF + ELSE +* +* Form C := alpha*B*A + beta*C. +* + DO 170 J = 1,N + TEMP1 = ALPHA*REAL(A(J,J)) + IF (BETA.EQ.ZERO) THEN + DO 110 I = 1,M + C(I,J) = TEMP1*B(I,J) + 110 CONTINUE + ELSE + DO 120 I = 1,M + C(I,J) = BETA*C(I,J) + TEMP1*B(I,J) + 120 CONTINUE + END IF + DO 140 K = 1,J - 1 + IF (UPPER) THEN + TEMP1 = ALPHA*A(K,J) + ELSE + TEMP1 = ALPHA*CONJG(A(J,K)) + END IF + DO 130 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 130 CONTINUE + 140 CONTINUE + DO 160 K = J + 1,N + IF (UPPER) THEN + TEMP1 = ALPHA*CONJG(A(J,K)) + ELSE + TEMP1 = ALPHA*A(K,J) + END IF + DO 150 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + END IF +* + RETURN +* +* End of CHEMM . +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/dsymm.f Fri Feb 06 07:30:55 2009 +0100 @@ -0,0 +1,294 @@ + SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* .. Scalar Arguments .. + DOUBLE PRECISION ALPHA,BETA + INTEGER LDA,LDB,LDC,M,N + CHARACTER SIDE,UPLO +* .. +* .. Array Arguments .. + DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* Purpose +* ======= +* +* DSYMM performs one of the matrix-matrix operations +* +* C := alpha*A*B + beta*C, +* +* or +* +* C := alpha*B*A + beta*C, +* +* where alpha and beta are scalars, A is a symmetric matrix and B and +* C are m by n matrices. +* +* Arguments +* ========== +* +* SIDE - CHARACTER*1. +* On entry, SIDE specifies whether the symmetric matrix A +* appears on the left or right in the operation as follows: +* +* SIDE = 'L' or 'l' C := alpha*A*B + beta*C, +* +* SIDE = 'R' or 'r' C := alpha*B*A + beta*C, +* +* Unchanged on exit. +* +* UPLO - CHARACTER*1. +* On entry, UPLO specifies whether the upper or lower +* triangular part of the symmetric matrix A is to be +* referenced as follows: +* +* UPLO = 'U' or 'u' Only the upper triangular part of the +* symmetric matrix is to be referenced. +* +* UPLO = 'L' or 'l' Only the lower triangular part of the +* symmetric matrix is to be referenced. +* +* Unchanged on exit. +* +* M - INTEGER. +* On entry, M specifies the number of rows 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 C. +* 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, ka ), where ka is +* m when SIDE = 'L' or 'l' and is n otherwise. +* Before entry with SIDE = 'L' or 'l', the m by m part of +* the array A must contain the symmetric matrix, such that +* when UPLO = 'U' or 'u', the leading m by m upper triangular +* part of the array A must contain the upper triangular part +* of the symmetric matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading m by m lower triangular part of the array A +* must contain the lower triangular part of the symmetric +* matrix and the strictly upper triangular part of A is not +* referenced. +* Before entry with SIDE = 'R' or 'r', the n by n part of +* the array A must contain the symmetric matrix, such that +* when UPLO = 'U' or 'u', the leading n by n upper triangular +* part of the array A must contain the upper triangular part +* of the symmetric matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading n by n lower triangular part of the array A +* must contain the lower triangular part of the symmetric +* matrix and the strictly upper triangular part of A is not +* referenced. +* 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 ), otherwise 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. +* Unchanged on exit. +* +* 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. +* +* 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 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, 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 .. + DOUBLE PRECISION TEMP1,TEMP2 + INTEGER I,INFO,J,K,NROWA + LOGICAL UPPER +* .. +* .. Parameters .. + DOUBLE PRECISION ONE,ZERO + PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) +* .. +* +* Set NROWA as the number of rows of A. +* + IF (LSAME(SIDE,'L')) THEN + NROWA = M + ELSE + NROWA = N + END IF + UPPER = LSAME(UPLO,'U') +* +* Test the input parameters. +* + INFO = 0 + IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN + INFO = 1 + ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN + INFO = 2 + ELSE IF (M.LT.0) THEN + INFO = 3 + ELSE IF (N.LT.0) THEN + INFO = 4 + ELSE IF (LDA.LT.MAX(1,NROWA)) THEN + INFO = 7 + ELSE IF (LDB.LT.MAX(1,M)) THEN + INFO = 9 + ELSE IF (LDC.LT.MAX(1,M)) THEN + INFO = 12 + END IF + IF (INFO.NE.0) THEN + CALL XERBLA('DSYMM ',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 +* +* And when 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 (LSAME(SIDE,'L')) THEN +* +* Form C := alpha*A*B + beta*C. +* + IF (UPPER) THEN + DO 70 J = 1,N + DO 60 I = 1,M + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 50 K = 1,I - 1 + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*A(K,I) + 50 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + + + ALPHA*TEMP2 + END IF + 60 CONTINUE + 70 CONTINUE + ELSE + DO 100 J = 1,N + DO 90 I = M,1,-1 + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 80 K = I + 1,M + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*A(K,I) + 80 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + + + ALPHA*TEMP2 + END IF + 90 CONTINUE + 100 CONTINUE + END IF + ELSE +* +* Form C := alpha*B*A + beta*C. +* + DO 170 J = 1,N + TEMP1 = ALPHA*A(J,J) + IF (BETA.EQ.ZERO) THEN + DO 110 I = 1,M + C(I,J) = TEMP1*B(I,J) + 110 CONTINUE + ELSE + DO 120 I = 1,M + C(I,J) = BETA*C(I,J) + TEMP1*B(I,J) + 120 CONTINUE + END IF + DO 140 K = 1,J - 1 + IF (UPPER) THEN + TEMP1 = ALPHA*A(K,J) + ELSE + TEMP1 = ALPHA*A(J,K) + END IF + DO 130 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 130 CONTINUE + 140 CONTINUE + DO 160 K = J + 1,N + IF (UPPER) THEN + TEMP1 = ALPHA*A(J,K) + ELSE + TEMP1 = ALPHA*A(K,J) + END IF + DO 150 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + END IF +* + RETURN +* +* End of DSYMM . +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/ssymm.f Fri Feb 06 07:30:55 2009 +0100 @@ -0,0 +1,294 @@ + SUBROUTINE SSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* .. Scalar Arguments .. + REAL ALPHA,BETA + INTEGER LDA,LDB,LDC,M,N + CHARACTER SIDE,UPLO +* .. +* .. Array Arguments .. + REAL A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* Purpose +* ======= +* +* SSYMM performs one of the matrix-matrix operations +* +* C := alpha*A*B + beta*C, +* +* or +* +* C := alpha*B*A + beta*C, +* +* where alpha and beta are scalars, A is a symmetric matrix and B and +* C are m by n matrices. +* +* Arguments +* ========== +* +* SIDE - CHARACTER*1. +* On entry, SIDE specifies whether the symmetric matrix A +* appears on the left or right in the operation as follows: +* +* SIDE = 'L' or 'l' C := alpha*A*B + beta*C, +* +* SIDE = 'R' or 'r' C := alpha*B*A + beta*C, +* +* Unchanged on exit. +* +* UPLO - CHARACTER*1. +* On entry, UPLO specifies whether the upper or lower +* triangular part of the symmetric matrix A is to be +* referenced as follows: +* +* UPLO = 'U' or 'u' Only the upper triangular part of the +* symmetric matrix is to be referenced. +* +* UPLO = 'L' or 'l' Only the lower triangular part of the +* symmetric matrix is to be referenced. +* +* Unchanged on exit. +* +* M - INTEGER. +* On entry, M specifies the number of rows 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 C. +* N must be at least zero. +* Unchanged on exit. +* +* ALPHA - REAL . +* On entry, ALPHA specifies the scalar alpha. +* Unchanged on exit. +* +* A - REAL array of DIMENSION ( LDA, ka ), where ka is +* m when SIDE = 'L' or 'l' and is n otherwise. +* Before entry with SIDE = 'L' or 'l', the m by m part of +* the array A must contain the symmetric matrix, such that +* when UPLO = 'U' or 'u', the leading m by m upper triangular +* part of the array A must contain the upper triangular part +* of the symmetric matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading m by m lower triangular part of the array A +* must contain the lower triangular part of the symmetric +* matrix and the strictly upper triangular part of A is not +* referenced. +* Before entry with SIDE = 'R' or 'r', the n by n part of +* the array A must contain the symmetric matrix, such that +* when UPLO = 'U' or 'u', the leading n by n upper triangular +* part of the array A must contain the upper triangular part +* of the symmetric matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading n by n lower triangular part of the array A +* must contain the lower triangular part of the symmetric +* matrix and the strictly upper triangular part of A is not +* referenced. +* 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 ), otherwise LDA must be at +* least max( 1, n ). +* Unchanged on exit. +* +* B - REAL array of DIMENSION ( LDB, n ). +* Before entry, the leading m by n 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. LDB must be at least +* max( 1, m ). +* Unchanged on exit. +* +* BETA - REAL . +* 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 - REAL 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 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, 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 .. + REAL TEMP1,TEMP2 + INTEGER I,INFO,J,K,NROWA + LOGICAL UPPER +* .. +* .. Parameters .. + REAL ONE,ZERO + PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) +* .. +* +* Set NROWA as the number of rows of A. +* + IF (LSAME(SIDE,'L')) THEN + NROWA = M + ELSE + NROWA = N + END IF + UPPER = LSAME(UPLO,'U') +* +* Test the input parameters. +* + INFO = 0 + IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN + INFO = 1 + ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN + INFO = 2 + ELSE IF (M.LT.0) THEN + INFO = 3 + ELSE IF (N.LT.0) THEN + INFO = 4 + ELSE IF (LDA.LT.MAX(1,NROWA)) THEN + INFO = 7 + ELSE IF (LDB.LT.MAX(1,M)) THEN + INFO = 9 + ELSE IF (LDC.LT.MAX(1,M)) THEN + INFO = 12 + END IF + IF (INFO.NE.0) THEN + CALL XERBLA('SSYMM ',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 +* +* And when 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 (LSAME(SIDE,'L')) THEN +* +* Form C := alpha*A*B + beta*C. +* + IF (UPPER) THEN + DO 70 J = 1,N + DO 60 I = 1,M + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 50 K = 1,I - 1 + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*A(K,I) + 50 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + + + ALPHA*TEMP2 + END IF + 60 CONTINUE + 70 CONTINUE + ELSE + DO 100 J = 1,N + DO 90 I = M,1,-1 + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 80 K = I + 1,M + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*A(K,I) + 80 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) + + + ALPHA*TEMP2 + END IF + 90 CONTINUE + 100 CONTINUE + END IF + ELSE +* +* Form C := alpha*B*A + beta*C. +* + DO 170 J = 1,N + TEMP1 = ALPHA*A(J,J) + IF (BETA.EQ.ZERO) THEN + DO 110 I = 1,M + C(I,J) = TEMP1*B(I,J) + 110 CONTINUE + ELSE + DO 120 I = 1,M + C(I,J) = BETA*C(I,J) + TEMP1*B(I,J) + 120 CONTINUE + END IF + DO 140 K = 1,J - 1 + IF (UPPER) THEN + TEMP1 = ALPHA*A(K,J) + ELSE + TEMP1 = ALPHA*A(J,K) + END IF + DO 130 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 130 CONTINUE + 140 CONTINUE + DO 160 K = J + 1,N + IF (UPPER) THEN + TEMP1 = ALPHA*A(J,K) + ELSE + TEMP1 = ALPHA*A(K,J) + END IF + DO 150 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + END IF +* + RETURN +* +* End of SSYMM . +* + END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/blas/zhemm.f Fri Feb 06 07:30:55 2009 +0100 @@ -0,0 +1,298 @@ + SUBROUTINE ZHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* .. Scalar Arguments .. + DOUBLE COMPLEX ALPHA,BETA + INTEGER LDA,LDB,LDC,M,N + CHARACTER SIDE,UPLO +* .. +* .. Array Arguments .. + DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* Purpose +* ======= +* +* ZHEMM performs one of the matrix-matrix operations +* +* C := alpha*A*B + beta*C, +* +* or +* +* C := alpha*B*A + beta*C, +* +* where alpha and beta are scalars, A is an hermitian matrix and B and +* C are m by n matrices. +* +* Arguments +* ========== +* +* SIDE - CHARACTER*1. +* On entry, SIDE specifies whether the hermitian matrix A +* appears on the left or right in the operation as follows: +* +* SIDE = 'L' or 'l' C := alpha*A*B + beta*C, +* +* SIDE = 'R' or 'r' C := alpha*B*A + beta*C, +* +* Unchanged on exit. +* +* UPLO - CHARACTER*1. +* On entry, UPLO specifies whether the upper or lower +* triangular part of the hermitian matrix A is to be +* referenced as follows: +* +* UPLO = 'U' or 'u' Only the upper triangular part of the +* hermitian matrix is to be referenced. +* +* UPLO = 'L' or 'l' Only the lower triangular part of the +* hermitian matrix is to be referenced. +* +* Unchanged on exit. +* +* M - INTEGER. +* On entry, M specifies the number of rows 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 C. +* N must be at least zero. +* Unchanged on exit. +* +* ALPHA - COMPLEX*16 . +* On entry, ALPHA specifies the scalar alpha. +* Unchanged on exit. +* +* A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is +* m when SIDE = 'L' or 'l' and is n otherwise. +* Before entry with SIDE = 'L' or 'l', the m by m part of +* the array A must contain the hermitian matrix, such that +* when UPLO = 'U' or 'u', the leading m by m upper triangular +* part of the array A must contain the upper triangular part +* of the hermitian matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading m by m lower triangular part of the array A +* must contain the lower triangular part of the hermitian +* matrix and the strictly upper triangular part of A is not +* referenced. +* Before entry with SIDE = 'R' or 'r', the n by n part of +* the array A must contain the hermitian matrix, such that +* when UPLO = 'U' or 'u', the leading n by n upper triangular +* part of the array A must contain the upper triangular part +* of the hermitian matrix and the strictly lower triangular +* part of A is not referenced, and when UPLO = 'L' or 'l', +* the leading n by n lower triangular part of the array A +* must contain the lower triangular part of the hermitian +* matrix and the strictly upper triangular part of A is not +* referenced. +* Note that the imaginary parts of the diagonal elements need +* not be set, they are assumed to be zero. +* 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 ), otherwise LDA must be at +* least max( 1, n ). +* Unchanged on exit. +* +* B - COMPLEX*16 array of DIMENSION ( LDB, n ). +* Before entry, the leading m by n 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. LDB must be at least +* max( 1, m ). +* Unchanged on exit. +* +* BETA - COMPLEX*16 . +* 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 - COMPLEX*16 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 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, 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 DBLE,DCONJG,MAX +* .. +* .. Local Scalars .. + DOUBLE COMPLEX TEMP1,TEMP2 + INTEGER I,INFO,J,K,NROWA + LOGICAL UPPER +* .. +* .. Parameters .. + DOUBLE COMPLEX ONE + PARAMETER (ONE= (1.0D+0,0.0D+0)) + DOUBLE COMPLEX ZERO + PARAMETER (ZERO= (0.0D+0,0.0D+0)) +* .. +* +* Set NROWA as the number of rows of A. +* + IF (LSAME(SIDE,'L')) THEN + NROWA = M + ELSE + NROWA = N + END IF + UPPER = LSAME(UPLO,'U') +* +* Test the input parameters. +* + INFO = 0 + IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN + INFO = 1 + ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN + INFO = 2 + ELSE IF (M.LT.0) THEN + INFO = 3 + ELSE IF (N.LT.0) THEN + INFO = 4 + ELSE IF (LDA.LT.MAX(1,NROWA)) THEN + INFO = 7 + ELSE IF (LDB.LT.MAX(1,M)) THEN + INFO = 9 + ELSE IF (LDC.LT.MAX(1,M)) THEN + INFO = 12 + END IF + IF (INFO.NE.0) THEN + CALL XERBLA('ZHEMM ',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 +* +* And when 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 (LSAME(SIDE,'L')) THEN +* +* Form C := alpha*A*B + beta*C. +* + IF (UPPER) THEN + DO 70 J = 1,N + DO 60 I = 1,M + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 50 K = 1,I - 1 + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I)) + 50 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) + + + ALPHA*TEMP2 + END IF + 60 CONTINUE + 70 CONTINUE + ELSE + DO 100 J = 1,N + DO 90 I = M,1,-1 + TEMP1 = ALPHA*B(I,J) + TEMP2 = ZERO + DO 80 K = I + 1,M + C(K,J) = C(K,J) + TEMP1*A(K,I) + TEMP2 = TEMP2 + B(K,J)*DCONJG(A(K,I)) + 80 CONTINUE + IF (BETA.EQ.ZERO) THEN + C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2 + ELSE + C(I,J) = BETA*C(I,J) + TEMP1*DBLE(A(I,I)) + + + ALPHA*TEMP2 + END IF + 90 CONTINUE + 100 CONTINUE + END IF + ELSE +* +* Form C := alpha*B*A + beta*C. +* + DO 170 J = 1,N + TEMP1 = ALPHA*DBLE(A(J,J)) + IF (BETA.EQ.ZERO) THEN + DO 110 I = 1,M + C(I,J) = TEMP1*B(I,J) + 110 CONTINUE + ELSE + DO 120 I = 1,M + C(I,J) = BETA*C(I,J) + TEMP1*B(I,J) + 120 CONTINUE + END IF + DO 140 K = 1,J - 1 + IF (UPPER) THEN + TEMP1 = ALPHA*A(K,J) + ELSE + TEMP1 = ALPHA*DCONJG(A(J,K)) + END IF + DO 130 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 130 CONTINUE + 140 CONTINUE + DO 160 K = J + 1,N + IF (UPPER) THEN + TEMP1 = ALPHA*DCONJG(A(J,K)) + ELSE + TEMP1 = ALPHA*A(K,J) + END IF + DO 150 I = 1,M + C(I,J) = C(I,J) + TEMP1*B(I,K) + 150 CONTINUE + 160 CONTINUE + 170 CONTINUE + END IF +* + RETURN +* +* End of ZHEMM . +* + END