C ALGORITHM 784, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 24,NO. 3, September, 1998, P. 303--316. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Package/ # Package/GBBEN/ # Package/GBBEN/CBENCH/ # Package/GBBEN/CBENCH/Makefile # Package/GBBEN/CBENCH/cgb02.f # Package/GBBEN/CBENCH/cgb03.f # Package/GBBEN/CBENCH/cgb04.f # Package/GBBEN/CBENCH/cgb05.f # Package/GBBEN/CBENCH/cgb06.f # Package/GBBEN/CBENCH/cgb07.f # Package/GBBEN/CBENCH/cgb08.f # Package/GBBEN/CBENCH/cgb09.f # Package/GBBEN/CBENCH/cgb90.f # Package/GBBEN/CBENCH/cgb91.f # Package/GBBEN/CBENCH/cgbt01.f # Package/GBBEN/CBENCH/cgbt02.f # Package/GBBEN/CBENCH/cgbtim.f # Package/GBBEN/CBENCH/cgbtp1.f # Package/GBBEN/CBENCH/cgbtp2.f # Package/GBBEN/CBENCH/cmark01.in # Package/GBBEN/CBENCH/cmark02.in # Package/GBBEN/CBENCH/csbpm.f # Package/GBBEN/CBENCH/eoln.f # Package/GBBEN/CBENCH/example.in # Package/GBBEN/CBENCH/getwrd.f # Package/GBBEN/CBENCH/lsame.f # Package/GBBEN/CBENCH/newcgpm.in # Package/GBBEN/CBENCH/xerbla.f # Package/GBBEN/DBENCH/ # Package/GBBEN/DBENCH/Makefile # Package/GBBEN/DBENCH/dgb02.f # Package/GBBEN/DBENCH/dgb04.f # Package/GBBEN/DBENCH/dgb06.f # Package/GBBEN/DBENCH/dgb08.f # Package/GBBEN/DBENCH/dgb09.f # Package/GBBEN/DBENCH/dgb90.f # Package/GBBEN/DBENCH/dgb91.f # Package/GBBEN/DBENCH/dgbt01.f # Package/GBBEN/DBENCH/dgbt02.f # Package/GBBEN/DBENCH/dgbtim.f # Package/GBBEN/DBENCH/dgbtp1.f # Package/GBBEN/DBENCH/dgbtp2.f # Package/GBBEN/DBENCH/dmark01.in # Package/GBBEN/DBENCH/dmark02.in # Package/GBBEN/DBENCH/dsbpm.f # Package/GBBEN/DBENCH/eoln.f # Package/GBBEN/DBENCH/example.in # Package/GBBEN/DBENCH/getwrd.f # Package/GBBEN/DBENCH/lsame.f # Package/GBBEN/DBENCH/newdgpm.in # Package/GBBEN/DBENCH/xerbla.f # Package/GBBEN/Makefile # Package/GBBEN/SBENCH/ # Package/GBBEN/SBENCH/Makefile # Package/GBBEN/SBENCH/eoln.f # Package/GBBEN/SBENCH/example.in # Package/GBBEN/SBENCH/getwrd.f # Package/GBBEN/SBENCH/lsame.f # Package/GBBEN/SBENCH/newsgpm.in # Package/GBBEN/SBENCH/sgb02.f # Package/GBBEN/SBENCH/sgb04.f # Package/GBBEN/SBENCH/sgb06.f # Package/GBBEN/SBENCH/sgb08.f # Package/GBBEN/SBENCH/sgb09.f # Package/GBBEN/SBENCH/sgb90.f # Package/GBBEN/SBENCH/sgb91.f # Package/GBBEN/SBENCH/sgbt01.f # Package/GBBEN/SBENCH/sgbt02.f # Package/GBBEN/SBENCH/sgbtim.f # Package/GBBEN/SBENCH/sgbtp1.f # Package/GBBEN/SBENCH/sgbtp2.f # Package/GBBEN/SBENCH/smark01.in # Package/GBBEN/SBENCH/smark02.in # Package/GBBEN/SBENCH/ssbpm.f # Package/GBBEN/SBENCH/xerbla.f # Package/GBBEN/TMGLIB/ # Package/GBBEN/TMGLIB/Makefile # Package/GBBEN/TMGLIB/dsecnd.f # Package/GBBEN/TMGLIB/second.f # Package/GBBEN/ZBENCH/ # Package/GBBEN/ZBENCH/Makefile # Package/GBBEN/ZBENCH/eoln.f # Package/GBBEN/ZBENCH/example.in # Package/GBBEN/ZBENCH/getwrd.f # Package/GBBEN/ZBENCH/lsame.f # Package/GBBEN/ZBENCH/newzgpm.in # Package/GBBEN/ZBENCH/xerbla.f # Package/GBBEN/ZBENCH/zgb02.f # Package/GBBEN/ZBENCH/zgb03.f # Package/GBBEN/ZBENCH/zgb04.f # Package/GBBEN/ZBENCH/zgb05.f # Package/GBBEN/ZBENCH/zgb06.f # Package/GBBEN/ZBENCH/zgb07.f # Package/GBBEN/ZBENCH/zgb08.f # Package/GBBEN/ZBENCH/zgb09.f # Package/GBBEN/ZBENCH/zgb90.f # Package/GBBEN/ZBENCH/zgb91.f # Package/GBBEN/ZBENCH/zgbt01.f # Package/GBBEN/ZBENCH/zgbt02.f # Package/GBBEN/ZBENCH/zgbtim.f # Package/GBBEN/ZBENCH/zgbtp1.f # Package/GBBEN/ZBENCH/zgbtp2.f # Package/GBBEN/ZBENCH/zmark01.in # Package/GBBEN/ZBENCH/zmark02.in # Package/GBBEN/ZBENCH/zsbpm.f # Package/GBBEN/make.gbinc # Package/GBL3B/ # Package/GBL3B/CGBL3B/ # Package/GBL3B/CGBL3B/Makefile # Package/GBL3B/CGBL3B/cbigp.f # Package/GBL3B/CGBL3B/ccld.f # Package/GBL3B/CGBL3B/cgpm.in # Package/GBL3B/CGBL3B/chemm.f # Package/GBL3B/CGBL3B/cher2k.f # Package/GBL3B/CGBL3B/cherk.f # Package/GBL3B/CGBL3B/csgpm.f # Package/GBL3B/CGBL3B/csymm.f # Package/GBL3B/CGBL3B/csyr2k.f # Package/GBL3B/CGBL3B/csyrk.f # Package/GBL3B/CGBL3B/ctrmm.f # Package/GBL3B/CGBL3B/ctrsm.f # Package/GBL3B/CGBL3B/eoln.f # Package/GBL3B/CGBL3B/getwrd.f # Package/GBL3B/CGBL3B/lsame.f # Package/GBL3B/CGBL3B/xerbla.f # Package/GBL3B/DGBL3B/ # Package/GBL3B/DGBL3B/Makefile # Package/GBL3B/DGBL3B/dbigp.f # Package/GBL3B/DGBL3B/dcld.f # Package/GBL3B/DGBL3B/dgpm.in # Package/GBL3B/DGBL3B/dsgpm.f # Package/GBL3B/DGBL3B/dsymm.f # Package/GBL3B/DGBL3B/dsyr2k.f # Package/GBL3B/DGBL3B/dsyrk.f # Package/GBL3B/DGBL3B/dtrmm.f # Package/GBL3B/DGBL3B/dtrsm.f # Package/GBL3B/DGBL3B/eoln.f # Package/GBL3B/DGBL3B/getwrd.f # Package/GBL3B/DGBL3B/lsame.f # Package/GBL3B/DGBL3B/xerbla.f # Package/GBL3B/Makefile # Package/GBL3B/SGBL3B/ # Package/GBL3B/SGBL3B/Makefile # Package/GBL3B/SGBL3B/eoln.f # Package/GBL3B/SGBL3B/getwrd.f # Package/GBL3B/SGBL3B/lsame.f # Package/GBL3B/SGBL3B/sbigp.f # Package/GBL3B/SGBL3B/scld.f # Package/GBL3B/SGBL3B/sgpm.in # Package/GBL3B/SGBL3B/ssgpm.f # Package/GBL3B/SGBL3B/ssymm.f # Package/GBL3B/SGBL3B/ssyr2k.f # Package/GBL3B/SGBL3B/ssyrk.f # Package/GBL3B/SGBL3B/strmm.f # Package/GBL3B/SGBL3B/strsm.f # Package/GBL3B/SGBL3B/xerbla.f # Package/GBL3B/ZGBL3B/ # Package/GBL3B/ZGBL3B/Makefile # Package/GBL3B/ZGBL3B/eoln.f # Package/GBL3B/ZGBL3B/getwrd.f # Package/GBL3B/ZGBL3B/lsame.f # Package/GBL3B/ZGBL3B/xerbla.f # Package/GBL3B/ZGBL3B/zbigp.f # Package/GBL3B/ZGBL3B/zcld.f # Package/GBL3B/ZGBL3B/zgpm.in # Package/GBL3B/ZGBL3B/zhemm.f # Package/GBL3B/ZGBL3B/zher2k.f # Package/GBL3B/ZGBL3B/zherk.f # Package/GBL3B/ZGBL3B/zsgpm.f # Package/GBL3B/ZGBL3B/zsymm.f # Package/GBL3B/ZGBL3B/zsyr2k.f # Package/GBL3B/ZGBL3B/zsyrk.f # Package/GBL3B/ZGBL3B/ztrmm.f # Package/GBL3B/ZGBL3B/ztrsm.f # Package/GBL3B/make.gbinc # Package/INSTALL/ # Package/INSTALL/Makefile # Package/INSTALL/dsecnd.f # Package/INSTALL/dsecndtst.f # Package/INSTALL/lsame.f # Package/INSTALL/lsametst.f # Package/INSTALL/second.f # Package/INSTALL/secondtst.f # Package/Makefile # Package/README # Package/TESTING/ # Package/TESTING/Makefile # Package/TESTING/SRC/ # Package/TESTING/SRC/Makefile # Package/TESTING/SRC/cblat3.f # Package/TESTING/SRC/dblat3.f # Package/TESTING/SRC/sblat3.f # Package/TESTING/SRC/zblat3.f # Package/TESTING/cblat3.in # Package/TESTING/dblat3.in # Package/TESTING/sblat3.in # Package/TESTING/zblat3.in # UNDERLIB/ # UNDERLIB/Makefile # UNDERLIB/caxpy.f # UNDERLIB/ccopy.f # UNDERLIB/cgemm.f # UNDERLIB/cgemv.f # UNDERLIB/cher.f # UNDERLIB/cscal.f # UNDERLIB/ctrmv.f # UNDERLIB/ctrsv.f # UNDERLIB/daxpy.f # UNDERLIB/dcabs1.f # UNDERLIB/dcopy.f # UNDERLIB/dgemm.f # UNDERLIB/dgemv.f # UNDERLIB/dger.f # UNDERLIB/dscal.f # UNDERLIB/dsyr.f # UNDERLIB/dtrmv.f # UNDERLIB/dtrsv.f # UNDERLIB/lsame.f # UNDERLIB/saxpy.f # UNDERLIB/scopy.f # UNDERLIB/sgemm.f # UNDERLIB/sgemv.f # UNDERLIB/sger.f # UNDERLIB/sscal.f # UNDERLIB/ssyr.f # UNDERLIB/strmv.f # UNDERLIB/strsv.f # UNDERLIB/xerbla.f # UNDERLIB/zaxpy.f # UNDERLIB/zcopy.f # UNDERLIB/zgemm.f # UNDERLIB/zgemv.f # UNDERLIB/zher.f # UNDERLIB/zscal.f # UNDERLIB/ztrmv.f # UNDERLIB/ztrsv.f # make.inc # This archive created: Tue Mar 23 08:55:20 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Package' then mkdir 'Package' fi cd 'Package' if test ! -d 'GBBEN' then mkdir 'GBBEN' fi cd 'GBBEN' if test ! -d 'CBENCH' then mkdir 'CBENCH' fi cd 'CBENCH' if test -f 'Makefile' then echo shar: will not over-write existing file "'Makefile'" else cat << SHAR_EOF > 'Makefile' include ../make.gbinc ### GEMM-Based Level 3 BLAS Benchmark #################################### # # The following libraries are specified, a user specified level 3 BLAS # library to be timed (LIB3B), a library with underlying BLAS routines # (LIB12B) where the underlying BLAS routine CGEMM may be specified # separately, and the library with the timing functions SECOND and # DSECND (CSEC). # LIB3B = $(ULIB) CGEMM = $(UULIB) LIB12B = $(UULIB) CSEC = $(UTMG) # # LIB specifies the order in which the libraries are linked with the # benchmark programs. Notice that the built-in GEMM-based routines # will be linked the first CGEMM, level 1 and 2 BLAS routines found # as underlying routines. Change the order in which the libraries are # linked as desired. # LIB = $(CSEC) $(LIB3B) $(CGEMM) $(LIB12B) # ### Compiler options ##################################################### # # The compiler options are specified as follows for the different # programs and libraries: # # CBMFLG : the GEMM-based performance benchmark programs # CPRFLG : routines that print the output results # CGBFLG : the built-in GEMM-based level 3 BLAS routines # CAXFLG : GEMM-based specific auxiliary routines # AXOPT : other auxiliary routines # CBMFLG = $(GBBOPT) CPRFLG = $(AXOPT) CGBFLG = $(GBOPT) CAXFLG = $(GBOPT) AXFLG = $(AXOPT) # ######################################################################## CTIMS = cgbtim.f CTIM = cgbtim.o CBMS = cgbt01.f cgbt02.f CBM = cgbt01.o cgbt02.o CPRS = cgbtp1.f cgbtp2.f CPR = cgbtp1.o cgbtp2.o CGBS = cgb02.f cgb03.f cgb04.f cgb05.f cgb06.f cgb07.f cgb08.f cgb09.f CGB = cgb02.o cgb03.o cgb04.o cgb05.o cgb06.o cgb07.o cgb08.o cgb09.o CAUXS = cgb90.f cgb91.f CAUX = cgb90.o cgb91.o AUXS = lsame.f xerbla.f AUX = lsame.o xerbla.o CPRMS = csbpm.f CPRM = csbpm.o AUXS2 = getwrd.f eoln.f AUX2 = getwrd.o eoln.o OBJ1 = $(CTIM) $(CBM) $(CGB) $(CAUX) $(AUX) $(AUX2) $(CPR) OBJ2 = $(CPRM) $(AUX2) ######################################################################## all: cgbtim csbpm cgbtim: $(OBJ1) $(LOADER) $(LOADOPT) -o cgbtim $(OBJ1) $(LIB) csbpm: $(OBJ2) $(LOADER) $(LOADOPT) -o csbpm $(OBJ2) $(CTIM): $(CTIMS) $(FORTRAN) -c $(CBMFLG) $(CTIMS) $(CBM): $(CBMS) $(FORTRAN) -c $(CBMFLG) $(CBMS) $(CPR): $(CPRS) $(FORTRAN) -c $(CPRFLG) $(CPRS) $(CGB): $(CGBS) $(FORTRAN) -c $(CGBFLG) $(CGBS) $(CAUX): $(CAUXS) $(FORTRAN) -c $(CAXFLG) $(CAUXS) $(AUX): $(AUXS) $(FORTRAN) -c $(AXFLG) $(AUXS) $(CPRM): $(CPRMS) $(FORTRAN) -c $(AXFLG) $(CPRMS) $(AUX2): $(AUXS2) $(FORTRAN) -c $(AXFLG) $(AUXS2) clean: rm -f *.o cgbtim csbpm SHAR_EOF fi # end of overwriting check if test -f 'cgb02.f' then echo shar: will not over-write existing file "'cgb02.f'" else cat << SHAR_EOF > 'cgb02.f' SUBROUTINE CGB02 ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO INTEGER M, N, LDA, LDB, LDC COMPLEX ALPHA, BETA * .. Array Arguments .. COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * CGB02 (CSYMM) 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. * * Parameters * ========== * * 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 - 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 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 - 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. * * -- Rewritten in May-1994. * GEMM-Based Level 3 BLAS. * Per Ling, Institute of Information Processing, * University of Umea, Sweden. * * * .. Local Scalars .. INTEGER INFO, NROWA LOGICAL LSIDE, UPPER INTEGER I, II, IX, ISEC, J, JJ, JX, JSEC * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA EXTERNAL CGEMM, CCOPY * .. Parameters .. COMPLEX ZERO, ONE PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ), $ ONE = ( 1.0E+0, 0.0E+0 ) ) * .. User specified parameters for CGB02 .. INTEGER RCB, CB PARAMETER ( RCB = 128, CB = 64 ) * .. Local Arrays .. COMPLEX T1( RCB, RCB ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE, 'L' ) UPPER = LSAME( UPLO, 'U' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF 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( 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( 'CGB02 ', 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 CALL CGEMM ( 'N', 'N', M, N, 0, ZERO, A, MAX( LDA, LDB ), $ B, MAX( LDA, LDB ), BETA, C, LDC ) RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( UPPER )THEN * * Form C := alpha*A*B + beta*C. Left, Upper. * DO 40, II = 1, M, RCB ISEC = MIN( RCB, M-II+1 ) * * T1 := A, a upper triangular diagonal block of A is copied * to the upper triangular part of T1. * DO 10, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, A( II, I ), 1, T1( 1, I-II+1 ), $ 1 ) 10 CONTINUE * * T1 := A', a strictly upper triangular diagonal block of * A is copied to the strictly lower triangular part of T1. * Notice that T1 is referenced by row and that the maximum * length of a vector referenced by CCOPY is CB. * DO 30, JJ = II, II+ISEC-1, CB JSEC = MIN( CB, II+ISEC-JJ ) DO 20, J = JJ+1, II+ISEC-1 CALL CCOPY ( MIN( JSEC, J-JJ ), A( JJ, J ), 1, $ T1( J-II+1, JJ-II+1 ), RCB ) 20 CONTINUE 30 CONTINUE * * C := alpha*T1*B + beta*C, a horizontal block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', ISEC, N, ISEC, ALPHA, T1( 1, 1 ), $ RCB, B( II, 1 ), LDB, BETA, C( II, 1 ), LDC ) * * C := alpha*A'*B + C and C := alpha*A*B + C, general * matrix multiply operations involving rectangular blocks * of A. * IF( II.GT.1 )THEN CALL CGEMM ( 'T', 'N', ISEC, N, II-1, ALPHA, $ A( 1, II ), LDA, B( 1, 1 ), LDB, $ ONE, C( II, 1 ), LDC ) END IF IF( II+ISEC.LE.M )THEN CALL CGEMM ( 'N', 'N', ISEC, N, M-II-ISEC+1, ALPHA, $ A( II, II+ISEC ), LDA, B( II+ISEC, 1 ), $ LDB, ONE, C( II, 1 ), LDC ) END IF 40 CONTINUE ELSE * * Form C := alpha*A*B + beta*C. Left, Lower. * DO 80, IX = M, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 * * T1 := A, a lower triangular diagonal block of A is copied * to the lower triangular part of T1. * DO 50, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, A( I, I ), 1, $ T1( I-II+1, I-II+1 ), 1 ) 50 CONTINUE * * T1 := A', a strictly lower triangular diagonal block of * A is copied to the strictly upper triangular part of T1. * Notice that T1 is referenced by row and that the maximum * length of a vector referenced by CCOPY is CB. * DO 70, JX = II+ISEC-1, II, -CB JJ = MAX( II, JX-CB+1 ) JSEC = JX-JJ+1 DO 60, J = II, JJ+JSEC-2 CALL CCOPY ( MIN( JSEC, JJ+JSEC-1-J ), $ A( MAX( JJ, J+1 ), J ), 1, $ T1( J-II+1, MAX( JJ-II+1, J-II+2 ) ), RCB ) 60 CONTINUE 70 CONTINUE * * C := alpha*T1*B + beta*C, a horizontal block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', ISEC, N, ISEC, ALPHA, T1( 1, 1 ), $ RCB, B( II, 1 ), LDB, BETA, C( II, 1 ), LDC ) * * C := alpha*A'*B + C and C := alpha*A*B + C, general * matrix multiply operations involving rectangular blocks * of A. * IF( II+ISEC.LE.M )THEN CALL CGEMM ( 'T', 'N', ISEC, N, M-II-ISEC+1, ALPHA, $ A( II+ISEC, II ), LDA, B( II+ISEC, 1 ), $ LDB, ONE, C( II, 1 ), LDC ) END IF IF( II.GT.1 )THEN CALL CGEMM ( 'N', 'N', ISEC, N, II-1, ALPHA, $ A( II, 1 ), LDA, B( 1, 1 ), LDB, $ ONE, C( II, 1 ), LDC ) END IF 80 CONTINUE END IF ELSE IF( UPPER )THEN * * Form C := alpha*B*A + beta*C. Right, Upper. * DO 120, JJ = 1, N, RCB JSEC = MIN( RCB, N-JJ+1 ) * * T1 := A, a upper triangular diagonal block of A is copied * to the upper triangular part of T1. * DO 90, J = JJ, JJ+JSEC-1 CALL CCOPY ( J-JJ+1, A( JJ, J ), 1, T1( 1, J-JJ+1 ), $ 1 ) 90 CONTINUE * * T1 := A', a strictly upper triangular diagonal block of * A is copied to the strictly lower triangular part of T1. * Notice that T1 is referenced by row and that the maximum * length of a vector referenced by CCOPY is CB. * DO 110, II = JJ, JJ+JSEC-1, CB ISEC = MIN( CB, JJ+JSEC-II ) DO 100, I = II+1, JJ+JSEC-1 CALL CCOPY ( MIN( ISEC, I-II ), A( II, I ), 1, $ T1( I-JJ+1, II-JJ+1 ), RCB ) 100 CONTINUE 110 CONTINUE * * C := alpha*B*T1 + beta*C, a vertical block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', M, JSEC, JSEC, ALPHA, B( 1, JJ ), $ LDB, T1( 1, 1 ), RCB, BETA, C( 1, JJ ), LDC ) * * C := alpha*B*A + C and C := alpha*B*A' + C, general * matrix multiply operations involving rectangular blocks * of A. * IF( JJ.GT.1 )THEN CALL CGEMM ( 'N', 'N', M, JSEC, JJ-1, ALPHA, $ B( 1, 1 ), LDB, A( 1, JJ ), LDA, $ ONE, C( 1, JJ ), LDC ) END IF IF( JJ+JSEC.LE.N )THEN CALL CGEMM ( 'N', 'T', M, JSEC, N-JJ-JSEC+1, ALPHA, $ B( 1, JJ+JSEC ), LDB, A( JJ, JJ+JSEC ), $ LDA, ONE, C( 1, JJ ), LDC ) END IF 120 CONTINUE ELSE * * Form C := alpha*B*A + beta*C. Right, Lower. * DO 160, JX = N, 1, -RCB JJ = MAX( 1, JX-RCB+1 ) JSEC = JX-JJ+1 * * T1 := A, a lower triangular diagonal block of A is copied * to the lower triangular part of T1. * DO 130, J = JJ, JJ+JSEC-1 CALL CCOPY ( JJ+JSEC-J, A( J, J ), 1, $ T1( J-JJ+1, J-JJ+1 ), 1 ) 130 CONTINUE * * T1 := A', a strictly lower triangular diagonal block of * A is copied to the strictly upper triangular part of T1. * Notice that T1 is referenced by row and that the maximum * length of a vector referenced by CCOPY is CB. * DO 150, IX = JJ+JSEC-1, JJ, -CB II = MAX( JJ, IX-CB+1 ) ISEC = IX-II+1 DO 140, I = JJ, II+ISEC-2 CALL CCOPY ( MIN( ISEC, II+ISEC-1-I ), $ A( MAX( II, I+1 ), I ), 1, $ T1( I-JJ+1, MAX( II-JJ+1, I-JJ+2 ) ), RCB ) 140 CONTINUE 150 CONTINUE * * C := alpha*B*T1 + beta*C, a vertical block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', M, JSEC, JSEC, ALPHA, $ B( 1, JJ ), LDB, T1( 1, 1 ), RCB, $ BETA, C( 1, JJ ), LDC ) * * C := alpha*B*A + C and C := alpha*B*A' + C, general * matrix multiply operations involving rectangular blocks * of A. * IF( JJ+JSEC.LE.N )THEN CALL CGEMM ( 'N', 'N', M, JSEC, N-JJ-JSEC+1, ALPHA, $ B( 1, JJ+JSEC ), LDB, A( JJ+JSEC, JJ ), $ LDA, ONE, C( 1, JJ ), LDC ) END IF IF( JJ.GT.1 )THEN CALL CGEMM ( 'N', 'T', M, JSEC, JJ-1, ALPHA, $ B( 1, 1 ), LDB, A( JJ, 1 ), LDA, $ ONE, C( 1, JJ ), LDC ) END IF 160 CONTINUE END IF END IF * RETURN * * End of CGB02. * END SHAR_EOF fi # end of overwriting check if test -f 'cgb03.f' then echo shar: will not over-write existing file "'cgb03.f'" else cat << SHAR_EOF > 'cgb03.f' SUBROUTINE CGB03 ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO INTEGER M, N, LDA, LDB, LDC COMPLEX ALPHA, BETA * .. Array Arguments .. COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * CGB03 (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. * * Parameters * ========== * * 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. * * -- Rewritten in May-1994. * GEMM-Based Level 3 BLAS. * Per Ling, Institute of Information Processing, * University of Umea, Sweden. * * * .. Local Scalars .. INTEGER INFO, NROWA LOGICAL LSIDE, UPPER INTEGER I, II, IX, ISEC, J, JJ, JX, JSEC * .. Intrinsic Functions .. INTRINSIC MAX, MIN, REAL, CMPLX, CONJG * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA EXTERNAL CGEMM, CCOPY * .. Parameters .. REAL ZERO COMPLEX CZERO, CONE PARAMETER ( ZERO = 0.0E+0, $ CZERO = ( 0.0E+0, 0.0E+0 ), $ CONE = ( 1.0E+0, 0.0E+0 ) ) * .. User specified parameters for CGB03 .. INTEGER RCB, CB PARAMETER ( RCB = 128, CB = 64 ) * .. Local Arrays .. COMPLEX T1( RCB, RCB ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE, 'L' ) UPPER = LSAME( UPLO, 'U' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF 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( 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( 'CGB03 ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.CZERO ).AND.( BETA.EQ.CONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.CZERO )THEN CALL CGEMM ( 'N', 'N', M, N, 0, CZERO, A, MAX( LDA, LDB ), $ B, MAX( LDA, LDB ), BETA, C, LDC ) RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( UPPER )THEN * * Form C := alpha*A*B + beta*C. Left, Upper. * DO 60, II = 1, M, RCB ISEC = MIN( RCB, M-II+1 ) * * T1 := A, a upper triangular diagonal block of A is copied * to the upper triangular part of T1. * DO 10, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, A( II, I ), 1, T1( 1, I-II+1 ), $ 1 ) 10 CONTINUE * * Set the imaginary part of diagonal elements of T1 * to zero. * DO 20, I = 1, ISEC T1( I, I ) = CMPLX( REAL( T1( I, I ) ), ZERO ) 20 CONTINUE * * T1 := conjg( A' ), the conjugated transpose of a * strictly upper triangular diagonal block of A is copied * to the strictly lower triangular part of T1. Notice that * T1 is referenced by row and that the maximum length of a * vector referenced is CB. * DO 50, JJ = II, II+ISEC-1, CB JSEC = MIN( CB, II+ISEC-JJ ) DO 40, J = JJ+1, II+ISEC-1 DO 30, I = JJ, J-1 T1( J-II+1, I-II+1 ) = CONJG( A( I, J ) ) 30 CONTINUE 40 CONTINUE 50 CONTINUE * * C := alpha*T1*B + beta*C, a horizontal block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', ISEC, N, ISEC, ALPHA, T1( 1, 1 ), $ RCB, B( II, 1 ), LDB, BETA, C( II, 1 ), LDC ) * * C := alpha*conjg( A' )*B + C and C := alpha*A*B + C, * matrix multiply operations involving rectangular blocks * of A. * IF( II.GT.1 )THEN CALL CGEMM ( 'C', 'N', ISEC, N, II-1, ALPHA, $ A( 1, II ), LDA, B( 1, 1 ), LDB, $ CONE, C( II, 1 ), LDC ) END IF IF( II+ISEC.LE.M )THEN CALL CGEMM ( 'N', 'N', ISEC, N, M-II-ISEC+1, ALPHA, $ A( II, II+ISEC ), LDA, B( II+ISEC, 1 ), $ LDB, CONE, C( II, 1 ), LDC ) END IF 60 CONTINUE ELSE * * Form C := alpha*A*B + beta*C. Left, Lower. * DO 120, IX = M, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 * * T1 := A, a lower triangular diagonal block of A is copied * to the lower triangular part of T1. * DO 70, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, A( I, I ), 1, $ T1( I-II+1, I-II+1 ), 1 ) 70 CONTINUE * * Set the imaginary part of diagonal elements of T1 * to zero. * DO 80, I = 1, ISEC T1( I, I ) = CMPLX( REAL( T1( I, I ) ), ZERO ) 80 CONTINUE * * T1 := conjg( A' ), the conjugated transpose of a * strictly lower triangular diagonal block of A is copied * to the strictly upper triangular part of T1. Notice that * T1 is referenced by row and that the maximum length of a * vector referenced is CB. * DO 110, JX = II+ISEC-1, II, -CB JJ = MAX( II, JX-CB+1 ) JSEC = JX-JJ+1 DO 100, J = II, JJ+JSEC-2 DO 90, I = J+1, II+ISEC-1 T1( J-II+1, I-II+1 ) = CONJG( A( I, J ) ) 90 CONTINUE 100 CONTINUE 110 CONTINUE * * C := alpha*T1*B + beta*C, a horizontal block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', ISEC, N, ISEC, ALPHA, T1( 1, 1 ), $ RCB, B( II, 1 ), LDB, BETA, C( II, 1 ), LDC ) * * C := alpha*conjg( A' )*B + C and C := alpha*A*B + C, * matrix multiply operations involving rectangular blocks * of A. * IF( II+ISEC.LE.M )THEN CALL CGEMM ( 'C', 'N', ISEC, N, M-II-ISEC+1, ALPHA, $ A( II+ISEC, II ), LDA, B( II+ISEC, 1 ), $ LDB, CONE, C( II, 1 ), LDC ) END IF IF( II.GT.1 )THEN CALL CGEMM ( 'N', 'N', ISEC, N, II-1, ALPHA, $ A( II, 1 ), LDA, B( 1, 1 ), LDB, $ CONE, C( II, 1 ), LDC ) END IF 120 CONTINUE END IF ELSE IF( UPPER )THEN * * Form C := alpha*B*A + beta*C. Right, Upper. * DO 180, JJ = 1, N, RCB JSEC = MIN( RCB, N-JJ+1 ) * * T1 := A, a upper triangular diagonal block of A is copied * to the upper triangular part of T1. * DO 130, J = JJ, JJ+JSEC-1 CALL CCOPY ( J-JJ+1, A( JJ, J ), 1, T1( 1, J-JJ+1 ), $ 1 ) 130 CONTINUE * * Set the imaginary part of diagonal elements of T1 * to zero. * DO 140, J = 1, JSEC T1( J, J ) = CMPLX( REAL( T1( J, J ) ), ZERO ) 140 CONTINUE * * T1 := conjg( A' ), the conjugated transpose of a * strictly upper triangular diagonal block of A is copied * to the strictly lower triangular part of T1. Notice that * T1 is referenced by row and that the maximum length of a * vector referenced is CB. * DO 170, II = JJ, JJ+JSEC-1, CB ISEC = MIN( CB, JJ+JSEC-II ) DO 160, I = II+1, JJ+JSEC-1 DO 150, J = II, I-1 T1( I-JJ+1, J-JJ+1 ) = CONJG( A( J, I ) ) 150 CONTINUE 160 CONTINUE 170 CONTINUE * * C := alpha*B*T1 + beta*C, a vertical block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', M, JSEC, JSEC, ALPHA, B( 1, JJ ), $ LDB, T1( 1, 1 ), RCB, BETA, C( 1, JJ ), LDC ) * * C := alpha*B*A + C and C := alpha*B*conjg( A' ) + C, * matrix multiply operations involving rectangular blocks * of A. * IF( JJ.GT.1 )THEN CALL CGEMM ( 'N', 'N', M, JSEC, JJ-1, ALPHA, $ B( 1, 1 ), LDB, A( 1, JJ ), LDA, $ CONE, C( 1, JJ ), LDC ) END IF IF( JJ+JSEC.LE.N )THEN CALL CGEMM ( 'N', 'C', M, JSEC, N-JJ-JSEC+1, ALPHA, $ B( 1, JJ+JSEC ), LDB, A( JJ, JJ+JSEC ), $ LDA, CONE, C( 1, JJ ), LDC ) END IF 180 CONTINUE ELSE * * Form C := alpha*B*A + beta*C. Right, Lower. * DO 240, JX = N, 1, -RCB JJ = MAX( 1, JX-RCB+1 ) JSEC = JX-JJ+1 * * T1 := A, a lower triangular diagonal block of A is copied * to the lower triangular part of T1. * DO 190, J = JJ, JJ+JSEC-1 CALL CCOPY ( JJ+JSEC-J, A( J, J ), 1, $ T1( J-JJ+1, J-JJ+1 ), 1 ) 190 CONTINUE * * Set the imaginary part of diagonal elements of T1 * to zero. * DO 200, J = 1, JSEC T1( J, J ) = CMPLX( REAL( T1( J, J ) ), ZERO ) 200 CONTINUE * * T1 := conjg( A' ), the conjugated transpose of a * strictly lower triangular diagonal block of A is copied * to the strictly upper triangular part of T1. Notice that * T1 is referenced by row and that the maximum length of a * vector referenced is CB. * DO 230, IX = JJ+JSEC-1, JJ, -CB II = MAX( JJ, IX-CB+1 ) ISEC = IX-II+1 DO 220, I = JJ, II+ISEC-2 DO 210, J = I+1, JJ+JSEC-1 T1( I-JJ+1, J-JJ+1 ) = CONJG( A( J, I ) ) 210 CONTINUE 220 CONTINUE 230 CONTINUE * * C := alpha*B*T1 + beta*C, a vertical block of C is * updated using the general matrix multiply, CGEMM. T1 * corresponds to a full diagonal block of the matrix A. * CALL CGEMM ( 'N', 'N', M, JSEC, JSEC, ALPHA, $ B( 1, JJ ), LDB, T1( 1, 1 ), RCB, $ BETA, C( 1, JJ ), LDC ) * * C := alpha*B*A + C and C := alpha*B*conjg( A' ) + C, * matrix multiply operations involving rectangular blocks * of A. * IF( JJ+JSEC.LE.N )THEN CALL CGEMM ( 'N', 'N', M, JSEC, N-JJ-JSEC+1, ALPHA, $ B( 1, JJ+JSEC ), LDB, A( JJ+JSEC, JJ ), $ LDA, CONE, C( 1, JJ ), LDC ) END IF IF( JJ.GT.1 )THEN CALL CGEMM ( 'N', 'C', M, JSEC, JJ-1, ALPHA, $ B( 1, 1 ), LDB, A( JJ, 1 ), LDA, $ CONE, C( 1, JJ ), LDC ) END IF 240 CONTINUE END IF END IF * RETURN * * End of CGB03. * END SHAR_EOF fi # end of overwriting check if test -f 'cgb04.f' then echo shar: will not over-write existing file "'cgb04.f'" else cat << SHAR_EOF > 'cgb04.f' SUBROUTINE CGB04 ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC COMPLEX ALPHA, BETA * .. Array Arguments .. COMPLEX A( LDA, * ), C( LDC, * ) * .. * * Purpose * ======= * * CGB04 (CSYRK) 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. * * 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', K specifies the number of rows of the * matrix A. K 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 * 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 - COMPLEX . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX 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. * * -- Rewritten in May-1994. * GEMM-Based Level 3 BLAS. * Per Ling, Institute of Information Processing, * University of Umea, Sweden. * * * .. Local Scalars .. INTEGER INFO, NROWA INTEGER I, II, IX, ISEC, L, LL, LSEC LOGICAL UPPER, NOTR, CLDA, SMALLN, TINYK COMPLEX DELTA * .. Intrinsic Functions .. INTRINSIC MIN, MAX * .. External Functions .. LOGICAL LSAME, CGB90, CGB91 EXTERNAL LSAME, CGB90, CGB91 * .. External Subroutines .. EXTERNAL XERBLA EXTERNAL CGEMM, CGEMV, CCOPY, CSCAL * .. Parameters .. COMPLEX ONE, ZERO INTEGER CIP41, CIP42 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), $ ZERO = ( 0.0E+0, 0.0E+0 ), $ CIP41 = 41, CIP42 = 42 ) * .. User specified parameters for CGB04 .. INTEGER RB, CB, RCB PARAMETER ( RCB = 64, RB = 64, CB = 64 ) * .. Local Arrays .. COMPLEX T1( RB, CB ), T2( RCB, RCB ), T3( RCB, RCB ) * .. * .. Executable Statements .. * * Test the input parameters. * UPPER = LSAME( UPLO, 'U' ) NOTR = LSAME( TRANS, 'N' ) IF( NOTR )THEN NROWA = N ELSE NROWA = K END IF INFO = 0 IF( ( .NOT.UPPER ).AND.( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( ( .NOT.NOTR ) ).AND.( .NOT.LSAME( TRANS, 'T' ) ) )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( 'CGB04 ', 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 or k.eq.0. * IF( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) )THEN IF( UPPER )THEN DO 10, I = 1, N CALL CSCAL ( I, BETA, C( 1, I ), 1 ) 10 CONTINUE ELSE DO 20, I = 1, N CALL CSCAL ( N-I+1, BETA, C( I, I ), 1 ) 20 CONTINUE END IF RETURN END IF * * Start the operations. * IF( UPPER )THEN IF( NOTR )THEN * * Form C := alpha*A*A' + beta*C. Upper, Notr. * SMALLN = .NOT.CGB90( CIP41 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP42 , N, K ) DO 110, II = 1, N, RCB ISEC = MIN( RCB, N-II+1 ) * * C := alpha*A*A' + beta*C, general matrix multiply on * upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'N', 'T', II-1, ISEC, K, ALPHA, $ A( 1, 1 ), LDA, A( II, 1 ), LDA, $ BETA, C( 1, II ), LDC ) END IF IF( TINYK )THEN * * C := beta*C, a upper triangular diagonal block * of C is updated with beta. * IF( BETA.NE.ONE )THEN DO 30, I = II, II+ISEC-1 CALL CSCAL ( I-II+1, BETA, C( II, I ), 1 ) 30 CONTINUE END IF * * C := alpha*A*A' + C, symmetric matrix multiply. * C is a symmetric diagonal block having upper * triangular storage format. * DO 50, I = II, II+ISEC-1 DO 40, L = 1, K CALL CAXPY ( I-II+1, ALPHA*A( I, L ), $ A( II, L ), 1, C( II, I ), 1 ) 40 CONTINUE 50 CONTINUE ELSE * * T2 := C, a upper triangular diagonal block of the * symmetric matrix C is copied to the upper * triangular part of T2. * DO 60, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, C( II, I ), 1, $ T2( 1, I-II+1 ), 1 ) 60 CONTINUE * * T2 := beta*T2, the upper triangular part of T2 is * updated with beta. * IF( BETA.NE.ONE )THEN DO 70, I = II, II+ISEC-1 CALL CSCAL ( I-II+1, BETA, $ T2( 1, I-II+1 ), 1 ) 70 CONTINUE END IF * * T2 := alpha*A*A' + T2, symmetric matrix multiply. * T2 contains a symmetric block having upper * triangular storage format. * DO 90, I = II, II+ISEC-1 DO 80, L = 1, K CALL CAXPY ( I-II+1, ALPHA*A( I, L ), $ A( II, L ), 1, T2( 1, I-II+1 ), 1 ) 80 CONTINUE 90 CONTINUE * * C := T2, the upper triangular part of T2 is copied * back to C. * DO 100, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, T2( 1, I-II+1 ), 1, $ C( II, I ), 1 ) 100 CONTINUE END IF 110 CONTINUE ELSE DO 150, II = 1, N, RB ISEC = MIN( RB, N-II+1 ) * * C := alpha*A*A' + beta*C, general matrix multiply on * upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'N', 'T', II-1, ISEC, K, ALPHA, $ A( 1, 1 ), LDA, A( II, 1 ), LDA, $ BETA, C( 1, II ), LDC ) END IF DELTA = BETA DO 140, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := A, a rectangular block of A is copied to T1. * DO 120, L = LL, LL+LSEC-1 CALL CCOPY ( ISEC, A( II, L ), 1, $ T1( 1, L-LL+1 ), 1 ) 120 CONTINUE * * C := alpha*T1*T1' + delta*C, C is symmetric having * triangular storage format. Delta is used instead * of beta to avoid updating the block of C with beta * multiple times. * DO 130, I = II, II+ISEC-1 CALL CGEMV ( 'N', I-II+1, LSEC, ALPHA, $ T1( 1, 1 ), RB, T1( I-II+1, 1 ), RB, $ DELTA, C( II, I ), 1 ) 130 CONTINUE DELTA = ONE 140 CONTINUE 150 CONTINUE END IF ELSE * * Form C := alpha*A'*A + beta*C. Upper, Trans. * SMALLN = .NOT.CGB90( CIP41 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP42 , N, K ) DO 260, II = 1, N, RCB ISEC = MIN( RCB, N-II+1 ) * * C := alpha*A'*A + beta*C, general matrix multiply on * upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'T', 'N', II-1, ISEC, K, ALPHA, $ A( 1, 1 ), LDA, A( 1, II ), LDA, $ BETA, C( 1, II ), LDC ) END IF IF( TINYK )THEN * * C := beta*C, a upper triangular diagonal block * of C is updated with beta. * IF( BETA.NE.ONE )THEN DO 160, I = II, II+ISEC-1 CALL CSCAL ( I-II+1, BETA, C( II, I ), 1 ) 160 CONTINUE END IF * * C := alpha*A*A' + C, symmetric matrix multiply. * C is a symmetric diagonal block having upper * triangular storage format. * DO 180, I = II, II+ISEC-1 DO 170, L = 1, K CALL CAXPY ( I-II+1, ALPHA*A( L, I ), $ A( L, II ), LDA, C( II, I ), 1 ) 170 CONTINUE 180 CONTINUE ELSE * * T2 := C, a upper triangular diagonal block of the * symmetric matrix C is copied to the upper * triangular part of T2. * DO 190, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, C( II, I ), 1, $ T2( 1, I-II+1 ), 1 ) 190 CONTINUE * * T2 := beta*T2, the upper triangular part of T2 is * updated with beta. * IF( BETA.NE.ONE )THEN DO 200, I = II, II+ISEC-1 CALL CSCAL ( I-II+1, BETA, $ T2( 1, I-II+1 ), 1 ) 200 CONTINUE END IF DO 240, LL = 1, K, RCB LSEC = MIN( RCB, K-LL+1 ) * * T3 := A', the transpose of a square block of A * is copied to T3. * DO 210, I = II, II+ISEC-1 CALL CCOPY ( LSEC, A( LL, I ), 1, $ T3( I-II+1, 1 ), RCB ) 210 CONTINUE * * T2 := alpha*T3*T3' + T2, symmetric matrix * multiply. T2 contains a symmetric block having * upper triangular storage format. * DO 230, I = II, II+ISEC-1 DO 220, L = LL, LL+LSEC-1 CALL CAXPY ( I-II+1, $ ALPHA*T3( I-II+1, L-LL+1 ), $ T3( 1, L-LL+1 ), 1, $ T2( 1, I-II+1 ), 1 ) 220 CONTINUE 230 CONTINUE 240 CONTINUE * * C := T2, the upper triangular part of T2 is copied * back to C. * DO 250, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, T2( 1, I-II+1 ), 1, $ C( II, I ), 1 ) 250 CONTINUE END IF 260 CONTINUE ELSE CLDA = CGB91( LDA ) DO 310, II = 1, N, RB ISEC = MIN( RB, N-II+1 ) * * C := alpha*A'*A + beta*C, general matrix multiply on * upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'T', 'N', II-1, ISEC, K, ALPHA, $ A( 1, 1 ), LDA, A( 1, II ), LDA, $ BETA, C( 1, II ), LDC ) END IF DELTA = BETA DO 300, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := A', the transpose of a rectangular block * of A is copied to T1. * IF( CLDA )THEN DO 270, I = II, II+ISEC-1 CALL CCOPY ( LSEC, A( LL, I ), 1, $ T1( I-II+1, 1 ), RB ) 270 CONTINUE ELSE DO 280, L = LL, LL+LSEC-1 CALL CCOPY ( ISEC, A( L, II ), LDA, $ T1( 1, L-LL+1 ), 1 ) 280 CONTINUE END IF * * C := alpha*T1*T1' + delta*C, C is symmetric having * triangular storage format. Delta is used instead * of beta to avoid updating the block of C with beta * multiple times. * DO 290, I = II, II+ISEC-1 CALL CGEMV ( 'N', I-II+1, LSEC, ALPHA, $ T1( 1, 1 ), RB, T1( I-II+1, 1 ), RB, $ DELTA, C( II, I ), 1 ) 290 CONTINUE DELTA = ONE 300 CONTINUE 310 CONTINUE END IF END IF ELSE IF( NOTR )THEN * * Form C := alpha*A*A' + beta*C. Lower, Notr. * SMALLN = .NOT.CGB90( CIP41 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP42 , N, K ) DO 400, IX = N, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 IF( TINYK )THEN * * C := beta*C, a lower triangular diagonal block * of C is updated with beta. * IF( BETA.NE.ONE )THEN DO 320, I = II, II+ISEC-1 CALL CSCAL ( II+ISEC-I, BETA, C( I, I ), 1 ) 320 CONTINUE END IF * * C := alpha*A*A' + C, symmetric matrix multiply. * C is a symmetric diagonal block having lower * triangular storage format. * DO 340, I = II, II+ISEC-1 DO 330, L = 1, K CALL CAXPY ( II+ISEC-I, ALPHA*A( I, L ), $ A( I, L ), 1, C( I, I ), 1 ) 330 CONTINUE 340 CONTINUE ELSE * * T2 := C, a lower triangular diagonal block of the * symmetric matrix C is copied to the lower * triangular part of T2. * DO 350, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, C( I, I ), 1, $ T2( I-II+1, I-II+1 ), 1 ) 350 CONTINUE * * T2 := beta*T2, the lower triangular part of T2 is * updated with beta. * IF( BETA.NE.ONE )THEN DO 360, I = II, II+ISEC-1 CALL CSCAL ( II+ISEC-I, BETA, $ T2( I-II+1, I-II+1 ), 1 ) 360 CONTINUE END IF * * T2 := alpha*A*A' + T2, symmetric matrix multiply. * T2 contains a symmetric block having lower * triangular storage format. * DO 380, I = II, II+ISEC-1 DO 370, L = 1, K CALL CAXPY ( II+ISEC-I, ALPHA*A( I, L ), $ A( I, L ), 1, T2( I-II+1, I-II+1 ), 1 ) 370 CONTINUE 380 CONTINUE * * C := T2, the lower triangular part of T2 is copied * back to C. * DO 390, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, T2( I-II+1, I-II+1 ), $ 1, C( I, I ), 1 ) 390 CONTINUE END IF * * C := alpha*A*A' + beta*C, general matrix multiply on * lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'N', 'T', N-II-ISEC+1, ISEC, K, $ ALPHA, A( II+ISEC, 1 ), LDA, A( II, 1 ), $ LDA, BETA, C( II+ISEC, II ), LDC ) END IF 400 CONTINUE ELSE DO 440, IX = N, 1, -RB II = MAX( 1, IX-RB+1 ) ISEC = IX-II+1 DELTA = BETA DO 430, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := A, a rectangular block of A is copied to T1. * DO 410, L = LL, LL+LSEC-1 CALL CCOPY ( ISEC, A( II, L ), 1, $ T1( 1, L-LL+1 ), 1 ) 410 CONTINUE * * C := alpha*T1*T1' + delta*C, C is symmetric having * triangular storage format. Delta is used instead * of beta to avoid updating the block of C with beta * multiple times. * DO 420, I = II, II+ISEC-1 CALL CGEMV ( 'N', II+ISEC-I, LSEC, ALPHA, $ T1( I-II+1, 1 ), RB, T1( I-II+1, 1 ), RB, $ DELTA, C( I, I ), 1 ) 420 CONTINUE DELTA = ONE 430 CONTINUE * * C := alpha*A*A' + beta*C, general matrix multiply on * lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'N', 'T', N-II-ISEC+1, ISEC, K, $ ALPHA, A( II+ISEC, 1 ), LDA, A( II, 1 ), $ LDA, BETA, C( II+ISEC, II ), LDC ) END IF 440 CONTINUE END IF ELSE * * Form C := alpha*A'*A + beta*C. Lower, Trans. * SMALLN = .NOT.CGB90( CIP41 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP42 , N, K ) DO 550, IX = N, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 IF( TINYK )THEN * * C := beta*C, a lower triangular diagonal block * of C is updated with beta. * IF( BETA.NE.ONE )THEN DO 450, I = II, II+ISEC-1 CALL CSCAL ( II+ISEC-I, BETA, C( I, I ), 1 ) 450 CONTINUE END IF * * C := alpha*A*A' + C, symmetric matrix multiply. * C is a symmetric diagonal block having lower * triangular storage format. * DO 470, I = II, II+ISEC-1 DO 460, L = 1, K CALL CAXPY ( II+ISEC-I, ALPHA*A( L, I ), $ A( L, I ), LDA, C( I, I ), 1 ) 460 CONTINUE 470 CONTINUE ELSE * * T2 := C, a lower triangular diagonal block of the * symmetric matrix C is copied to the lower * triangular part of T2. * DO 480, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, C( I, I ), 1, $ T2( I-II+1, I-II+1 ), 1 ) 480 CONTINUE * * T2 := beta*T2, the lower triangular part of T2 is * updated with beta. * IF( BETA.NE.ONE )THEN DO 490, I = II, II+ISEC-1 CALL CSCAL ( II+ISEC-I, BETA, $ T2( I-II+1, I-II+1 ), 1 ) 490 CONTINUE END IF DO 530, LL = 1, K, RCB LSEC = MIN( RCB, K-LL+1 ) * * T3 := A', the transpose of a square block of A * is copied to T3. * DO 500, I = II, II+ISEC-1 CALL CCOPY ( LSEC, A( LL, I ), 1, $ T3( I-II+1, 1 ), RCB ) 500 CONTINUE * * T2 := alpha*T3*T3' + T2, symmetric matrix * multiply. T2 contains a symmetric block having * lower triangular storage format. * DO 520, I = II, II+ISEC-1 DO 510, L = LL, LL+LSEC-1 CALL CAXPY ( II+ISEC-I, $ ALPHA*T3( I-II+1, L-LL+1 ), $ T3( I-II+1, L-LL+1 ), 1, $ T2( I-II+1, I-II+1 ), 1 ) 510 CONTINUE 520 CONTINUE 530 CONTINUE * * C := T2, the lower triangular part of T2 is copied * back to C. * DO 540, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, T2( I-II+1, I-II+1 ), $ 1, C( I, I ), 1 ) 540 CONTINUE END IF * * C := alpha*A'*A + beta*C, general matrix multiply on * lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'T', 'N', N-II-ISEC+1, ISEC, K, $ ALPHA, A( 1, II+ISEC ), LDA, A( 1, II ), $ LDA, BETA, C( II+ISEC, II ), LDC ) END IF 550 CONTINUE ELSE CLDA = CGB91( LDA ) DO 600, IX = N, 1, -RB II = MAX( 1, IX-RB+1 ) ISEC = IX-II+1 DELTA = BETA DO 590, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := A', the transpose of a rectangular block * of A is copied to T1. * IF( CLDA )THEN DO 560, I = II, II+ISEC-1 CALL CCOPY ( LSEC, A( LL, I ), 1, $ T1( I-II+1, 1 ), RB ) 560 CONTINUE ELSE DO 570, L = LL, LL+LSEC-1 CALL CCOPY ( ISEC, A( L, II ), LDA, $ T1( 1, L-LL+1 ), 1 ) 570 CONTINUE END IF * * C := alpha*T1*T1' + delta*C, C is symmetric having * triangular storage format. Delta is used instead * of beta to avoid updating the block of C with beta * multiple times. * DO 580, I = II, II+ISEC-1 CALL CGEMV ( 'N', II+ISEC-I, LSEC, ALPHA, $ T1( I-II+1, 1 ), RB, T1( I-II+1, 1 ), RB, $ DELTA, C( I, I ), 1 ) 580 CONTINUE DELTA = ONE 590 CONTINUE * * C := alpha*A'*A + beta*C, general matrix multiply on * lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'T', 'N', N-II-ISEC+1, ISEC, K, $ ALPHA, A( 1, II+ISEC ), LDA, A( 1, II ), $ LDA, BETA, C( II+ISEC, II ), LDC ) END IF 600 CONTINUE END IF END IF END IF * RETURN * * End of CGB04. * END SHAR_EOF fi # end of overwriting check if test -f 'cgb05.f' then echo shar: will not over-write existing file "'cgb05.f'" else cat << SHAR_EOF > 'cgb05.f' SUBROUTINE CGB05 ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC REAL ALPHA, BETA * .. Array Arguments .. COMPLEX A( LDA, * ), C( LDC, * ) * .. * * Purpose * ======= * * CGB05 (CHERK) 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 - REAL. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX 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 - REAL. * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX 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. * * -- Rewritten in May-1994. * GEMM-Based Level 3 BLAS. * Per Ling, Institute of Information Processing, * University of Umea, Sweden. * * * .. Local Scalars .. INTEGER INFO, NROWA INTEGER I, II, IX, ISEC, L, LL, LSEC LOGICAL UPPER, NOTR, CLDA, SMALLN, TINYK COMPLEX CALPHA, CBETA, CDELTA * .. Intrinsic Functions .. INTRINSIC MIN, MAX, REAL, CMPLX, CONJG * .. External Functions .. LOGICAL LSAME, CGB90, CGB91 EXTERNAL LSAME, CGB90, CGB91 * .. External Subroutines .. EXTERNAL XERBLA EXTERNAL CGEMM, CGEMV, CHER, CCOPY, CSCAL * .. Parameters .. REAL ONE, ZERO COMPLEX CONE INTEGER CIP51, CIP52 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0, $ CONE = ( 1.0E+0, 0.0E+0 ), $ CIP51 = 51, CIP52 = 52 ) * .. User specified parameters for CGB05 .. INTEGER RB, CB, RCB PARAMETER ( RCB = 64, RB = 64, CB = 64 ) * .. Local Arrays .. COMPLEX T1( RB, CB ), T2( RCB, RCB ), T3( RCB, RCB ), $ T4( CB ) * .. * .. Executable Statements .. * * Test the input parameters. * UPPER = LSAME( UPLO, 'U' ) NOTR = LSAME( TRANS, 'N' ) IF( NOTR )THEN NROWA = N ELSE NROWA = K END IF INFO = 0 IF( ( .NOT.UPPER ).AND.( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTR ).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( 'CGB05 ', 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 * CALPHA = CMPLX( ALPHA, ZERO ) CBETA = CMPLX( BETA, ZERO ) * * And when alpha.eq.zero or k.eq.0. * IF( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) )THEN IF( UPPER )THEN C( 1, 1 ) = CMPLX( BETA*REAL( C( 1, 1 ) ), ZERO ) DO 10, I = 2, N CALL CSCAL ( I-1, CBETA, C( 1, I ), 1 ) C( I, I ) = CMPLX( BETA*REAL( C( I, I ) ), ZERO ) 10 CONTINUE ELSE DO 20, I = 1, N-1 C( I, I ) = CMPLX( BETA*REAL( C( I, I ) ), ZERO ) CALL CSCAL ( N-I, CBETA, C( I+1, I ), 1 ) 20 CONTINUE C( N, N ) = CMPLX( BETA*REAL( C( N, N ) ), ZERO ) END IF RETURN END IF * * Start the operations. * IF( UPPER )THEN IF( NOTR )THEN * * Form C := alpha*A*conjg( A' ) + beta*C. Upper, Notr. * SMALLN = .NOT.CGB90( CIP51 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP52 , N, K ) DO 90, II = 1, N, RCB ISEC = MIN( RCB, N-II+1 ) * * C := alpha*A*conjg( A' ) + beta*C, matrix multiply * updating upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'N', 'C', II-1, ISEC, K, CALPHA, $ A( 1, 1 ), LDA, A( II, 1 ), LDA, $ CBETA, C( 1, II ), LDC ) END IF IF( TINYK )THEN * * C := beta*C, a upper triangular diagonal block * of C is updated with beta. The imaginary part of * the diagonal elements of C are set to ZERO. * IF( BETA.NE.ONE )THEN C( II, II ) = $ CMPLX( BETA*REAL( C( II, II ) ), ZERO ) DO 30, I = II+1, II+ISEC-1 CALL CSCAL ( I-II, CBETA, C( II, I ), 1 ) C( I, I ) = $ CMPLX( BETA*REAL( C( I, I ) ), ZERO ) 30 CONTINUE END IF * * C := alpha*A*conjg( A' ) + C, hermitian matrix * multiply. C is a hermitian diagonal block having * upper triangular storage format. * DO 40, L = 1, K CALL CHER ( 'U', ISEC, ALPHA, A( II, L ), $ 1, C( II, II ), LDC ) 40 CONTINUE ELSE * * T2 := C, a upper triangular diagonal block of the * hermitian matrix C is copied to the upper * triangular part of T2. * DO 50, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, C( II, I ), 1, $ T2( 1, I-II+1 ), 1 ) 50 CONTINUE * * T2 := beta*T2, the upper triangular part of T2 is * updated with beta. The imaginary part of the * diagonal elements of T2 are set to ZERO. * IF( BETA.NE.ONE )THEN T2( 1, 1 ) = $ CMPLX( BETA*REAL( T2( 1, 1 ) ), ZERO ) DO 60, I = 2, ISEC CALL CSCAL ( I-1, CBETA, T2( 1, I ), 1 ) T2( I, I ) = $ CMPLX( BETA*REAL( T2( I, I ) ), ZERO ) 60 CONTINUE END IF * * T2 := alpha*A*conjg( A' ) + T2, hermitian matrix * multiply. T2 contains a hermitian block having * upper triangular storage format. * DO 70, L = 1, K CALL CHER ( 'U', ISEC, ALPHA, A( II, L ), $ 1, T2( 1, 1 ), RCB ) 70 CONTINUE * * C := T2, the upper triangular part of T2 is copied * back to C. * DO 80, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, T2( 1, I-II+1 ), 1, $ C( II, I ), 1 ) 80 CONTINUE END IF 90 CONTINUE ELSE DO 140, II = 1, N, RB ISEC = MIN( RB, N-II+1 ) * * C := alpha*A*conjg( A' ) + beta*C, matrix multiply * updating upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'N', 'C', II-1, ISEC, K, CALPHA, $ A( 1, 1 ), LDA, A( II, 1 ), LDA, $ CBETA, C( 1, II ), LDC ) END IF CDELTA = CBETA DO 130, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := A, a rectangular block of A is copied to T1. * DO 100, L = LL, LL+LSEC-1 CALL CCOPY ( ISEC, A( II, L ), 1, $ T1( 1, L-LL+1 ), 1 ) 100 CONTINUE * * C := alpha*T1*conjg( T1' ) + delta*C, C is * hermitian having triangular storage format. Delta * is used instead of beta to avoid updating the * block of C with beta multiple times. The local * array T4 is used for the conjugated transpose * of vectors of T1. * DO 120, I = II, II+ISEC-1 DO 110, L = LL, LL+LSEC-1 T4( L-LL+1 ) = $ CONJG( T1( I-II+1, L-LL+1 ) ) 110 CONTINUE CALL CGEMV ( 'N', I-II+1, LSEC, CALPHA, $ T1( 1, 1 ), RB, T4( 1 ), 1, $ CDELTA, C( II, I ), 1 ) C( I, I ) = CMPLX( REAL( C( I, I ) ), ZERO ) 120 CONTINUE CDELTA = CONE 130 CONTINUE 140 CONTINUE END IF ELSE * * Form C := alpha*conjg( A' )*A + beta*C. Upper, Trans. * SMALLN = .NOT.CGB90( CIP51 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP52 , N, K ) DO 250, II = 1, N, RCB ISEC = MIN( RCB, N-II+1 ) * * C := alpha*conjg( A' )*A + beta*C, matrix multiply * updating upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'C', 'N', II-1, ISEC, K, CALPHA, $ A( 1, 1 ), LDA, A( 1, II ), LDA, $ CBETA, C( 1, II ), LDC ) END IF IF( TINYK )THEN * * C := beta*C, a upper triangular diagonal block * of C is updated with beta. The imaginary part of * the diagonal elements of C are set to ZERO. * IF( BETA.NE.ONE )THEN C( II, II ) = $ CMPLX( BETA*REAL( C( II, II ) ), ZERO ) DO 150, I = II+1, II+ISEC-1 CALL CSCAL ( I-II, CBETA, C( II, I ), 1 ) C( I, I ) = $ CMPLX( BETA*REAL( C( I, I ) ), ZERO ) 150 CONTINUE END IF * * C := alpha*conjg( A' )*A + C, hermitian matrix * multiply. C is a hermitian diagonal block having * upper triangular storage format. The local array * T3 is used for temporary storage of the conjugate * transposed vectors of A. * DO 170, L = 1, K DO 160, I = II, II+ISEC-1 T3( I-II+1, 1 ) = CONJG( A( L, I ) ) 160 CONTINUE CALL CHER ( 'U', ISEC, ALPHA, T3( 1, 1 ), $ 1, C( II, II ), LDC ) 170 CONTINUE ELSE * * T2 := C, a upper triangular diagonal block of the * hermitian matrix C is copied to the upper * triangular part of T2. * DO 180, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, C( II, I ), 1, $ T2( 1, I-II+1 ), 1 ) 180 CONTINUE * * T2 := beta*T2, the upper triangular part of T2 is * updated with beta. * IF( BETA.NE.ONE )THEN DO 190, I = II, II+ISEC-1 CALL CSCAL ( I-II+1, CBETA, $ T2( 1, I-II+1 ), 1 ) 190 CONTINUE END IF DO 230, LL = 1, K, RCB LSEC = MIN( RCB, K-LL+1 ) * * T3 := A', the transpose of a square block of A * is copied to T3. * DO 200, I = II, II+ISEC-1 CALL CCOPY ( LSEC, A( LL, I ), 1, $ T3( I-II+1, 1 ), RCB ) 200 CONTINUE * * T2 := alpha*conjg( T3' )*T3 + T2, hermitian * matrix multiply. T2 contains a hermitian block * having upper triangular storage format. The * local array T3 is used for temporary storage of * the conjugate transposed vectors of A. * DO 220, L = LL, LL+LSEC-1 DO 210, I = 1, ISEC T3( I, L-LL+1 ) = $ CONJG( T3( I, L-LL+1 ) ) 210 CONTINUE CALL CHER ( 'U', ISEC, ALPHA, $ T3( 1, L-LL+1 ), 1, T2( 1, 1 ), RCB ) 220 CONTINUE 230 CONTINUE * * C := T2, the upper triangular part of T2 is copied * back to C. * DO 240, I = II, II+ISEC-1 CALL CCOPY ( I-II+1, T2( 1, I-II+1 ), 1, $ C( II, I ), 1 ) 240 CONTINUE END IF 250 CONTINUE ELSE CLDA = CGB91( LDA ) DO 330, II = 1, N, RB ISEC = MIN( RB, N-II+1 ) * * C := alpha*conjg( A' )*A + beta*C, matrix multiply * updating upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'C', 'N', II-1, ISEC, K, CALPHA, $ A( 1, 1 ), LDA, A( 1, II ), LDA, $ CBETA, C( 1, II ), LDC ) END IF CDELTA = CBETA DO 320, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := conjg( A' ), the conjugated transpose of a * rectangular block of A is copied to T1. * IF( CLDA )THEN DO 270, I = II, II+ISEC-1 DO 260, L = LL, LL+LSEC-1 T1( I-II+1, L-LL+1 ) = $ CONJG( A( L, I ) ) 260 CONTINUE 270 CONTINUE ELSE DO 290, L = LL, LL+LSEC-1 DO 280, I = II, II+ISEC-1 T1( I-II+1, L-LL+1 ) = $ CONJG( A( L, I ) ) 280 CONTINUE 290 CONTINUE END IF * * C := alpha*T1*conjg( T1' ) + delta*C, C is * hermitian having triangular storage format. Delta * is used instead of beta to avoid updating the * block of C with beta multiple times. The local * array T4 is used for the conjugated transpose * of vectors of T1. * DO 310, I = II, II+ISEC-1 DO 300, L = LL, LL+LSEC-1 T4( L-LL+1 ) = $ CONJG( T1( I-II+1, L-LL+1 ) ) 300 CONTINUE CALL CGEMV ( 'N', I-II+1, LSEC, CALPHA, $ T1( 1, 1 ), RB, T4( 1 ), 1, $ CDELTA, C( II, I ), 1 ) C( I, I ) = CMPLX( REAL( C( I, I ) ), ZERO ) 310 CONTINUE CDELTA = CONE 320 CONTINUE 330 CONTINUE END IF END IF ELSE IF( NOTR )THEN * * Form C := alpha*A*conjg( A' ) + beta*C. Lower, Notr. * SMALLN = .NOT.CGB90( CIP51 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP52 , N, K ) DO 400, IX = N, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 IF( TINYK )THEN * * C := beta*C, a lower triangular diagonal block * of C is updated with beta. * IF( BETA.NE.ONE )THEN DO 340, I = II, II+ISEC-2 C( I, I ) = $ CMPLX( BETA*REAL( C( I, I ) ), ZERO ) CALL CSCAL ( II+ISEC-I-1, CBETA, $ C( I+1, I ), 1 ) 340 CONTINUE C( II+ISEC-1, II+ISEC-1 ) = $ CMPLX( BETA*REAL( C( II+ISEC-1, $ II+ISEC-1 ) ), ZERO ) END IF * * C := alpha*A*conjg( A' ) + C, hermitian matrix * multiply. C is a hermitian diagonal block having * lower triangular storage format. * DO 350, L = 1, K CALL CHER ( 'L', ISEC, ALPHA, A( II, L ), $ 1, C( II, II ), LDC ) 350 CONTINUE ELSE * * T2 := C, a lower triangular diagonal block of the * hermitian matrix C is copied to the lower * triangular part of T2. * DO 360, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, C( I, I ), 1, $ T2( I-II+1, I-II+1 ), 1 ) 360 CONTINUE * * T2 := beta*T2, the lower triangular part of T2 is * updated with beta. The imaginary part of the * diagonal elements of T2 are set to ZERO. * IF( BETA.NE.ONE )THEN DO 370, I = 1, ISEC-1 T2( I, I ) = $ CMPLX( BETA*REAL( T2( I, I ) ), ZERO ) CALL CSCAL ( ISEC-I, CBETA, $ T2( I+1, I ), 1 ) 370 CONTINUE T2( ISEC, ISEC ) = $ CMPLX( BETA*REAL( T2( ISEC, ISEC ) ), $ ZERO ) END IF * * T2 := alpha*A*conjg( A' ) + T2, symmetric matrix * multiply. T2 contains a hermitian block having * lower triangular storage format. * DO 380, L = 1, K CALL CHER ( 'L', ISEC, ALPHA, A( II, L ), $ 1, T2( 1, 1 ), RCB ) 380 CONTINUE * * C := T2, the lower triangular part of T2 is copied * back to C. * DO 390, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, T2( I-II+1, I-II+1 ), $ 1, C( I, I ), 1 ) 390 CONTINUE END IF * * C := alpha*A*conjg( A' ) + beta*C, matrix multiply * on lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'N', 'C', N-II-ISEC+1, ISEC, K, $ CALPHA, A( II+ISEC, 1 ), LDA, A( II, 1 ), $ LDA, CBETA, C( II+ISEC, II ), LDC ) END IF 400 CONTINUE ELSE DO 450, IX = N, 1, -RB II = MAX( 1, IX-RB+1 ) ISEC = IX-II+1 CDELTA = CBETA DO 440, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := A, a rectangular block of A is copied to T1. * DO 410, L = LL, LL+LSEC-1 CALL CCOPY ( ISEC, A( II, L ), 1, $ T1( 1, L-LL+1 ), 1 ) 410 CONTINUE * * C := alpha*T1*conjg( T1' ) + delta*C, C is * hermitian having triangular storage format. Delta * is used instead of beta to avoid updating the * block of C with beta multiple times. The local * array T4 is used for the conjugated transpose * of vectors of T1. * DO 430, I = II, II+ISEC-1 DO 420, L = LL, LL+LSEC-1 T4( L-LL+1 ) = $ CONJG( T1( I-II+1, L-LL+1 ) ) 420 CONTINUE CALL CGEMV ( 'N', II+ISEC-I, LSEC, CALPHA, $ T1( I-II+1, 1 ), RB, T4( 1 ), 1, $ CDELTA, C( I, I ), 1 ) C( I, I ) = CMPLX( REAL( C( I, I ) ), ZERO ) 430 CONTINUE CDELTA = CONE 440 CONTINUE * * C := alpha*A*conjg( A' ) + beta*C, matrix multiply * updating lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'N', 'C', N-II-ISEC+1, ISEC, K, $ CALPHA, A( II+ISEC, 1 ), LDA, A( II, 1 ), $ LDA, CBETA, C( II+ISEC, II ), LDC ) END IF 450 CONTINUE END IF ELSE * * Form C := alpha*conjg( A' )*A + beta*C. Lower, Trans. * SMALLN = .NOT.CGB90( CIP51 , N, K ) IF( SMALLN )THEN TINYK = .NOT.CGB90( CIP52 , N, K ) DO 560, IX = N, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 IF( TINYK )THEN * * C := beta*C, a lower triangular diagonal block * of C is updated with beta. The imaginary part of * the diagonal elements of C are set to ZERO. * IF( BETA.NE.ONE )THEN DO 460, I = II, II+ISEC-2 C( I, I ) = $ CMPLX( BETA*REAL( C( I, I ) ), ZERO ) CALL CSCAL ( II+ISEC-I-1, CBETA, $ C( I+1, I ), 1 ) 460 CONTINUE C( II+ISEC-1, II+ISEC-1 ) = $ CMPLX( BETA*REAL( C( II+ISEC-1, $ II+ISEC-1 ) ), ZERO ) END IF * * C := alpha*conjg( A' )*A + C, hermitian matrix * multiply. C is a hermitian diagonal block having * lower triangular storage format. The local array * T3 is used for temporary storage of the conjugate * transposed vectors of A. * DO 480, L = 1, K DO 470, I = II, II+ISEC-1 T3( I-II+1, 1 ) = CONJG( A( L, I ) ) 470 CONTINUE CALL CHER ( 'L', ISEC, ALPHA, T3( 1, 1 ), $ 1, C( II, II ), LDC ) 480 CONTINUE ELSE * * T2 := C, a lower triangular diagonal block of the * symmetric matrix C is copied to the lower * triangular part of T2. * DO 490, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, C( I, I ), 1, $ T2( I-II+1, I-II+1 ), 1 ) 490 CONTINUE * * T2 := beta*T2, the lower triangular part of T2 is * updated with beta. * IF( BETA.NE.ONE )THEN DO 500, I = II, II+ISEC-1 CALL CSCAL ( II+ISEC-I, CBETA, $ T2( I-II+1, I-II+1 ), 1 ) 500 CONTINUE END IF DO 540, LL = 1, K, RCB LSEC = MIN( RCB, K-LL+1 ) * * T3 := A', the transpose of a square block of A * is copied to T3. * DO 510, I = II, II+ISEC-1 CALL CCOPY ( LSEC, A( LL, I ), 1, $ T3( I-II+1, 1 ), RCB ) 510 CONTINUE * * T2 := alpha*conjg( T3' )*T3 + T2, hermitian * matrix multiply. T2 contains a hermitian block * having lower triangular storage format. The * local array T3 is used for temporary storage of * the conjugate transposed vectors of A. * DO 530, L = LL, LL+LSEC-1 DO 520, I = 1, ISEC T3( I, L-LL+1 ) = $ CONJG( T3( I, L-LL+1 ) ) 520 CONTINUE CALL CHER ( 'L', ISEC, ALPHA, $ T3( 1, L-LL+1 ), 1, T2( 1, 1 ), RCB ) 530 CONTINUE 540 CONTINUE * * C := T2, the lower triangular part of T2 is copied * back to C. * DO 550, I = II, II+ISEC-1 CALL CCOPY ( II+ISEC-I, T2( I-II+1, I-II+1 ), $ 1, C( I, I ), 1 ) 550 CONTINUE END IF * * C := alpha*conjg( A' )*A + beta*C, matrix multiply * updating lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'C', 'N', N-II-ISEC+1, ISEC, K, $ CALPHA, A( 1, II+ISEC ), LDA, A( 1, II ), $ LDA, CBETA, C( II+ISEC, II ), LDC ) END IF 560 CONTINUE ELSE CLDA = CGB91( LDA ) DO 650, IX = N, 1, -RB II = MAX( 1, IX-RB+1 ) ISEC = IX-II+1 CDELTA = CBETA DO 640, LL = 1, K, CB LSEC = MIN( CB, K-LL+1 ) * * T1 := conjg( A' ), the conjugated transpose of a * rectangular block of A is copied to T1. * IF( CLDA )THEN DO 580, I = II, II+ISEC-1 DO 570, L = LL, LL+LSEC-1 T1( I-II+1, L-LL+1 ) = $ CONJG( A( L, I ) ) 570 CONTINUE 580 CONTINUE ELSE DO 600, L = LL, LL+LSEC-1 DO 590, I = II, II+ISEC-1 T1( I-II+1, L-LL+1 ) = $ CONJG( A( L, I ) ) 590 CONTINUE 600 CONTINUE END IF * * C := alpha*T1*conjg( T1' ) + delta*C, C is * hermitian having triangular storage format. Delta * is used instead of beta to avoid updating the * block of C with beta multiple times. The local * array T4 is used for the conjugated transpose * of vectors of T1. * DO 630, I = II, II+ISEC-1 DO 620, L = LL, LL+LSEC-1 T4( L-LL+1 ) = $ CONJG( T1( I-II+1, L-LL+1 ) ) 620 CONTINUE CALL CGEMV ( 'N', II+ISEC-I, LSEC, CALPHA, $ T1( I-II+1, 1 ), RB, T4( 1 ), 1, $ CDELTA, C( I, I ), 1 ) C( I, I ) = CMPLX( REAL( C( I, I ) ), ZERO ) 630 CONTINUE CDELTA = CONE 640 CONTINUE * * C := alpha*conjg( A' )*A + beta*C, matrix multiply * updating lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'C', 'N', N-II-ISEC+1, ISEC, K, $ CALPHA, A( 1, II+ISEC ), LDA, A( 1, II ), $ LDA, CBETA, C( II+ISEC, II ), LDC ) END IF 650 CONTINUE END IF END IF END IF * RETURN * * End of CGB05. * END SHAR_EOF fi # end of overwriting check if test -f 'cgb06.f' then echo shar: will not over-write existing file "'cgb06.f'" else cat << SHAR_EOF > 'cgb06.f' SUBROUTINE CGB06 ( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDB, LDC COMPLEX ALPHA, BETA * .. Array Arguments .. COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * CGB06 (CSYR2K) performs one of the symmetric rank 2k operations * * C := alpha*A*B' + alpha*B*A' + beta*C, * * or * * C := alpha*A'*B + alpha*B'*A + beta*C, * * where alpha and beta are scalars, C is an n by n symmetric matrix * and A and B are n by k matrices in the first case and k by n * matrices 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*B' + alpha*B*A' + * beta*C. * * TRANS = 'T' or 't' C := alpha*A'*B + alpha*B'*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 matrices A and B, and on entry with * TRANS = 'T' or 't', K specifies the number of rows of the * matrices A and B. K 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 * 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. * * B - COMPLEX array of DIMENSION ( LDB, kb ), where kb 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 B must contain the matrix B, otherwise * the leading k 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. When TRANS = 'N' or 'n' * then LDB must be at least max( 1, n ), otherwise LDB must * be at least max( 1, k ). * Unchanged on exit. * * BETA - COMPLEX . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX 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. * * -- Rewritten in May-1994. * GEMM-Based Level 3 BLAS. * Per Ling, Institute of Information Processing, * University of Umea, Sweden. * * * .. Local Scalars .. INTEGER INFO, NROWA INTEGER I, II, IX, ISEC, JJ, JX, JSEC LOGICAL UPPER, NOTR * .. Intrinsic Functions .. INTRINSIC MIN, MAX * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA EXTERNAL CGEMM, CAXPY, CSCAL * .. Parameters .. COMPLEX ONE, ZERO PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), $ ZERO = ( 0.0E+0, 0.0E+0 ) ) * .. User specified parameters for CGB06 .. INTEGER RCB, CB PARAMETER ( RCB = 128, CB = 64 ) * .. Local Arrays .. COMPLEX T1( RCB, RCB ) * .. * .. Executable Statements .. * * Test the input parameters. * UPPER = LSAME( UPLO, 'U' ) NOTR = LSAME( TRANS, 'N' ) IF( NOTR )THEN NROWA = N ELSE NROWA = K END IF INFO = 0 IF( ( .NOT.UPPER ).AND.( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTR ).AND. ( .NOT.LSAME( TRANS, 'T' ) ) )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( LDB.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 12 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'CGB06 ', 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 or k.eq.0. * IF( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) )THEN IF( UPPER )THEN DO 10, I = 1, N CALL CSCAL ( I, BETA, C( 1, I ), 1 ) 10 CONTINUE ELSE DO 20, I = 1, N CALL CSCAL ( N-I+1, BETA, C( I, I ), 1 ) 20 CONTINUE END IF RETURN END IF * * Start the operations. * IF( UPPER )THEN IF( NOTR )THEN * * Form C := alpha*A*B' + alpha*B*A' + beta*C. Upper, Notr. * DO 70, II = 1, N, RCB ISEC = MIN( RCB, N-II+1 ) * * T1 := alpha*A*B', general matrix multiply on rectangular * blocks of A and B. T1 is square. * CALL CGEMM ( 'N', 'T', ISEC, ISEC, K, ALPHA, A( II, 1 ), $ LDA, B( II, 1 ), LDB, ZERO, T1( 1, 1 ), RCB ) * * C := beta*C, a upper triangular diagonal block of C is * updated with beta. * IF( BETA.NE.ONE )THEN DO 30, I = II, II+ISEC-1 CALL CSCAL ( I-II+1, BETA, C( II, I ), 1 ) 30 CONTINUE END IF * * C := T1 + C, the upper triangular part of T1 is added to * the upper triangular diagonal block of C. * DO 40, I = II, II+ISEC-1 CALL CAXPY ( I-II+1, ONE, T1( 1, I-II+1 ), 1, $ C( II, I ), 1 ) 40 CONTINUE * * C := T1' + C, the transpose of the lower triangular part * of T1 is added to the upper triangular diagonal block * of C. Notice that T1 is referenced by row and that the * maximum length of a vector referenced by CAXPY is CB. * DO 60, JJ = II, II+ISEC-1, CB JSEC = MIN( CB, II+ISEC-JJ ) DO 50, I = JJ, II+ISEC-1 CALL CAXPY ( MIN( JSEC, I-JJ+1 ), ONE, $ T1( I-II+1, JJ-II+1 ), RCB, C( JJ, I ), 1 ) 50 CONTINUE 60 CONTINUE * * C := alpha*A*B' + beta*C and C := alpha*B*A' + C, * general matrix multiply on upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'N', 'T', II-1, ISEC, K, ALPHA, $ A( 1, 1 ), LDA, B( II, 1 ), LDB, BETA, $ C( 1, II ), LDC ) CALL CGEMM ( 'N', 'T', II-1, ISEC, K, ALPHA, $ B( 1, 1 ), LDB, A( II, 1 ), LDA, ONE, $ C( 1, II ), LDC ) END IF 70 CONTINUE ELSE * * Form C := alpha*A'*B + alpha*B'*A + beta*C. Upper, Trans. * DO 120, II = 1, N, RCB ISEC = MIN( RCB, N-II+1 ) * * T1 := alpha*A'*B, general matrix multiply on rectangular * blocks of A and B. T1 is square. * CALL CGEMM ( 'T', 'N', ISEC, ISEC, K, ALPHA, A( 1, II ), $ LDA, B( 1, II ), LDB, ZERO, T1( 1, 1 ), RCB ) * * C := beta*C, a upper triangular diagonal block of C is * updated with beta. * IF( BETA.NE.ONE )THEN DO 80, I = II, II+ISEC-1 CALL CSCAL ( I-II+1, BETA, C( II, I ), 1 ) 80 CONTINUE END IF * * C := T1 + C, the upper triangular part of T1 is added to * the upper triangular diagonal block of C. * DO 90, I = II, II+ISEC-1 CALL CAXPY ( I-II+1, ONE, T1( 1, I-II+1 ), 1, $ C( II, I ), 1 ) 90 CONTINUE * * C := T1' + C, the transpose of the lower triangular part * of T1 is added to the upper triangular diagonal block * of C. Notice that T1 is referenced by row and that the * maximum length of a vector referenced by CAXPY is CB. * DO 110, JJ = II, II+ISEC-1, CB JSEC = MIN( CB, II+ISEC-JJ ) DO 100, I = JJ, II+ISEC-1 CALL CAXPY ( MIN( JSEC, I-JJ+1 ), ONE, $ T1( I-II+1, JJ-II+1 ), RCB, C( JJ, I ), 1 ) 100 CONTINUE 110 CONTINUE * * C := alpha*A'*B + beta*C and C := alpha*B'*A + C, * general matrix multiply on upper vertical blocks of C. * IF( II.GT.1 )THEN CALL CGEMM ( 'T', 'N', II-1, ISEC, K, ALPHA, $ A( 1, 1 ), LDA, B( 1, II ), LDB, BETA, $ C( 1, II ), LDC ) CALL CGEMM ( 'T', 'N', II-1, ISEC, K, ALPHA, $ B( 1, 1 ), LDB, A( 1, II ), LDA, ONE, $ C( 1, II ), LDC ) END IF 120 CONTINUE END IF ELSE IF( NOTR )THEN * * Form C := alpha*A*B' + alpha*B*A' + beta*C. Lower, Notr. * DO 170, IX = N, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 * * T1 := alpha*A*B', general matrix multiply on rectangular * blocks of A and B. T1 is square. * CALL CGEMM ( 'N', 'T', ISEC, ISEC, K, ALPHA, A( II, 1 ), $ LDA, B( II, 1 ), LDB, ZERO, T1( 1, 1 ), RCB ) * * C := beta*C, a lower triangular diagonal block of C is * updated with beta. * IF( BETA.NE.ONE )THEN DO 130, I = II, II+ISEC-1 CALL CSCAL ( II+ISEC-I, BETA, C( I, I ), 1 ) 130 CONTINUE END IF * * C := T1 + C, the lower triangular part of T1 is added to * the lower triangular diagonal block of C. * DO 140, I = II, II+ISEC-1 CALL CAXPY ( II+ISEC-I, ONE, T1( I-II+1, I-II+1 ), 1, $ C( I, I ), 1 ) 140 CONTINUE * * C := T1' + C, the transpose of the upper triangular part * of T1 is added to the lower triangular diagonal block * of C. Notice that T1 is referenced by row and that the * maximum length of a vector referenced by CAXPY is CB. * DO 160, JX = II+ISEC-1, II, -CB JJ = MAX( II, JX-CB+1 ) JSEC = JX-JJ+1 DO 150, I = II, JJ+JSEC-1 CALL CAXPY ( MIN( JSEC, JJ+JSEC-I ), ONE, $ T1( I-II+1, MAX( JJ-II+1, I-II+1 ) ), RCB, $ C( MAX( JJ, I ), I ), 1 ) 150 CONTINUE 160 CONTINUE * * C := alpha*A*B' + beta*C and C := alpha*B*A' + C, * general matrix multiply on lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'N', 'T', N-II-ISEC+1, ISEC, K, ALPHA, $ A( II+ISEC, 1 ), LDA, B( II, 1 ), LDB, $ BETA, C( II+ISEC, II ), LDC ) CALL CGEMM ( 'N', 'T', N-II-ISEC+1, ISEC, K, ALPHA, $ B( II+ISEC, 1 ), LDB, A( II, 1 ), LDA, $ ONE, C( II+ISEC, II ), LDC ) END IF 170 CONTINUE ELSE * * Form C := alpha*A'*B + alpha*B'*A + beta*C. Lower, Trans. * DO 220, IX = N, 1, -RCB II = MAX( 1, IX-RCB+1 ) ISEC = IX-II+1 * * T1 := alpha*A*B', general matrix multiply on rectangular * blocks of A and B. T1 is square. * CALL CGEMM ( 'T', 'N', ISEC, ISEC, K, ALPHA, A( 1, II ), $ LDA, B( 1, II ), LDB, ZERO, T1( 1, 1 ), RCB ) * * C := beta*C, a lower triangular diagonal block of C is * updated with beta. * IF( BETA.NE.ONE )THEN DO 180, I = II, II+ISEC-1 CALL CSCAL ( II+ISEC-I, BETA, C( I, I ), 1 ) 180 CONTINUE END IF * * C := T1 + C, the lower triangular part of T1 is added to * the lower triangular diagonal block of C. * DO 190, I = II, II+ISEC-1 CALL CAXPY ( II+ISEC-I, ONE, T1( I-II+1, I-II+1 ), 1, $ C( I, I ), 1 ) 190 CONTINUE * * C := T1' + C, the transpose of the upper triangular part * of T1 is added to the lower triangular diagonal block * of C. Notice that T1 is referenced by row and that the * maximum length of a vector referenced by CAXPY is CB. * DO 210, JX = II+ISEC-1, II, -CB JJ = MAX( II, JX-CB+1 ) JSEC = JX-JJ+1 DO 200, I = II, JJ+JSEC-1 CALL CAXPY ( MIN( JSEC, JJ+JSEC-I ), ONE, $ T1( I-II+1, MAX( JJ-II+1, I-II+1 ) ), RCB, $ C( MAX( JJ, I ), I ), 1 ) 200 CONTINUE 210 CONTINUE * * C := alpha*A'*B + beta*C and C := alpha*B'*A + C, * general matrix multiply on lower vertical blocks of C. * IF( II+ISEC.LE.N )THEN CALL CGEMM ( 'T', 'N', N-II-ISEC+1, ISEC, K, ALPHA, $ A( 1, II+ISEC ), LDA, B( 1, II ), LDB, $ BETA, C( II+ISEC, II ), LDC ) CALL CGEMM ( 'T', 'N', N-II-ISEC+1, ISEC, K, ALPHA, $ B( 1, II+ISEC ), LDB, A( 1, II ), LDA, $ ONE, C( II+ISEC, II ), LDC ) END IF 220 CONTINUE END IF END IF * RETURN * * End of CGB06. * END SHAR_EOF fi # end of overwriting check if test -f 'cgb07.f' then echo shar: will not over-write existing file "'cgb07.f'" else cat << SHAR_EOF > 'cgb07.f' SUBROUTINE CGB07 ( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDB, LDC REAL BETA COMPLEX ALPHA * .. Array Arguments .. COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * CGB07 (CHER2K) performs one of the hermitian rank 2k operations * * C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + beta*C, * * or * * C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + beta*C, * * where alpha and beta are scalars with beta real, C is an n by n * hermitian matrix and A and B are n by k matrices in the first case * and k by n matrices 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( B' ) + * conjg( alpha )*B*conjg( A' ) + * beta*C. * * TRANS = 'C' or 'c' C := alpha*conjg( A' )*B + * conjg( alpha )*conjg( B' )*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 matrices A and B, and on entry with * TRANS = 'C' or 'c', K specifies the number of rows of the * matrices A and B. K 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 * 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. * * B - COMPLEX array of DIMENSION ( LDB, kb ), where kb 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 B must contain the matrix B, otherwise * the leading k 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. When TRANS = 'N' or 'n' * then LDB must be at least max( 1, n ), otherwise LDB must * be at least max( 1, k ). * Unchanged on exit. * * BETA - REAL. * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX 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. * * -- Rewritten in May-1994. * GEMM-Based Level 3 BLAS. * Per Ling, Institute of Information Processing, * University of Umea, Sweden. * * * .. Local Scalars .. INTEGER INFO, NROWA INTEGER I, II, IX, ISEC, J, JJ, JX, JSEC LOGICAL UPPER, NOTR COMPLEX CBETA * .. Intrinsic Functions .. INTRINSIC MIN, MAX, REAL, CMPLX, CONJG * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA EXTERNAL CGEMM, CAXPY, CSCAL * .. Parameters .. REAL ONE, ZERO COMPLEX CONE, CZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0, $ CONE = ( 1.0E+0, 0.0E+0 ), $ CZERO = ( 0.0E+0, 0.0E+0 ) ) * .. User specified parameters for CGB07 .. INTEGER RCB, CB PARAMETER ( RCB = 128, CB = 64 ) * .. Local Arrays .. COMPLEX T1( RCB, RCB ) * .. * .. Executable Statements .. * * Test the input parameters. * UPPER = LSAME( UPLO, 'U' ) NOTR = LSAME( TRANS, 'N' ) IF( NOTR )THEN NROWA = N ELSE NROWA = K END IF INFO = 0 IF( ( .NOT.UPPER ).AND.( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTR ).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( LDB.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 12 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'CGB07 ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.CZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * CBET