double precision function dcabs1(z) double complex z,zz double precision t(2) equivalence (zz,t(1)) zz = z dcabs1 = dabs(t(1)) + dabs(t(2)) return end double precision function dzasum(n,zx,incx) c c takes the sum of the absolute values. c jack dongarra, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double complex zx(*) double precision stemp,dcabs1 integer i,incx,ix,n c dzasum = 0.0d0 stemp = 0.0d0 if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 do 10 i = 1,n stemp = stemp + dcabs1(zx(ix)) ix = ix + incx 10 continue dzasum = stemp return c c code for increment equal to 1 c 20 do 30 i = 1,n stemp = stemp + dcabs1(zx(i)) 30 continue dzasum = stemp return end DOUBLE PRECISION FUNCTION DZNRM2( N, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, N * .. Array Arguments .. COMPLEX*16 X( * ) * .. * * DZNRM2 returns the euclidean norm of a vector via the function * name, so that * * DZNRM2 := sqrt( conjg( x' )*x ) * * * * -- This version written on 25-October-1982. * Modified on 14-October-1993 to inline the call to ZLASSQ. * Sven Hammarling, Nag Ltd. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. INTEGER IX DOUBLE PRECISION NORM, SCALE, SSQ, TEMP * .. Intrinsic Functions .. INTRINSIC ABS, DIMAG, DBLE, SQRT * .. * .. Executable Statements .. IF( N.LT.1 .OR. INCX.LT.1 )THEN NORM = ZERO ELSE SCALE = ZERO SSQ = ONE * The following loop is equivalent to this call to the LAPACK * auxiliary routine: * CALL ZLASSQ( N, X, INCX, SCALE, SSQ ) * DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX IF( DBLE( X( IX ) ).NE.ZERO )THEN TEMP = ABS( DBLE( X( IX ) ) ) IF( SCALE.LT.TEMP )THEN SSQ = ONE + SSQ*( SCALE/TEMP )**2 SCALE = TEMP ELSE SSQ = SSQ + ( TEMP/SCALE )**2 END IF END IF IF( DIMAG( X( IX ) ).NE.ZERO )THEN TEMP = ABS( DIMAG( X( IX ) ) ) IF( SCALE.LT.TEMP )THEN SSQ = ONE + SSQ*( SCALE/TEMP )**2 SCALE = TEMP ELSE SSQ = SSQ + ( TEMP/SCALE )**2 END IF END IF 10 CONTINUE NORM = SCALE * SQRT( SSQ ) END IF * DZNRM2 = NORM RETURN * * End of DZNRM2. * END integer function izamax(n,zx,incx) c c finds the index of element having max. absolute value. c jack dongarra, 1/15/85. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double complex zx(*) double precision smax integer i,incx,ix,n double precision dcabs1 c izamax = 0 if( n.lt.1 .or. incx.le.0 )return izamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 smax = dcabs1(zx(1)) ix = ix + incx do 10 i = 2,n if(dcabs1(zx(ix)).le.smax) go to 5 izamax = i smax = dcabs1(zx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 smax = dcabs1(zx(1)) do 30 i = 2,n if(dcabs1(zx(i)).le.smax) go to 30 izamax = i smax = dcabs1(zx(i)) 30 continue return end subroutine zaxpy(n,za,zx,incx,zy,incy) c c constant times a vector plus a vector. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double complex zx(*),zy(*),za integer i,incx,incy,ix,iy,n double precision dcabs1 if(n.le.0)return if (dcabs1(za) .eq. 0.0d0) return if (incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n zy(iy) = zy(iy) + za*zx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n zy(i) = zy(i) + za*zx(i) 30 continue return end subroutine zcopy(n,zx,incx,zy,incy) c c copies a vector, x, to a vector, y. c jack dongarra, linpack, 4/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double complex zx(*),zy(*) integer i,incx,incy,ix,iy,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n zy(iy) = zx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n zy(i) = zx(i) 30 continue return end double complex function zdotc(n,zx,incx,zy,incy) c c forms the dot product of a vector. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n intrinsic dconjg ztemp = (0.0d0,0.0d0) zdotc = (0.0d0,0.0d0) if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n ztemp = ztemp + dconjg(zx(ix))*zy(iy) ix = ix + incx iy = iy + incy 10 continue zdotc = ztemp return c c code for both increments equal to 1 c 20 do 30 i = 1,n ztemp = ztemp + dconjg(zx(i))*zy(i) 30 continue zdotc = ztemp return end double complex function zdotu(n,zx,incx,zy,incy) c c forms the dot product of two vectors. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n ztemp = (0.0d0,0.0d0) zdotu = (0.0d0,0.0d0) if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n ztemp = ztemp + zx(ix)*zy(iy) ix = ix + incx iy = iy + incy 10 continue zdotu = ztemp return c c code for both increments equal to 1 c 20 do 30 i = 1,n ztemp = ztemp + zx(i)*zy(i) 30 continue zdotu = ztemp return end subroutine zdscal(n,da,zx,incx) c c scales a vector by a constant. c jack dongarra, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c intrinsic dcmplx double complex zx(*) double precision da integer i,incx,ix,n c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 do 10 i = 1,n zx(ix) = dcmplx(da,0.0d0)*zx(ix) ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 do 30 i = 1,n zx(i) = dcmplx(da,0.0d0)*zx(i) 30 continue return end SUBROUTINE ZGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. COMPLEX*16 ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or * * y := alpha*conjg( A' )*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - COMPLEX*16 array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - COMPLEX*16 array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY LOGICAL NOCONJ * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZGEMV ', 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 * NOCONJ = LSAME( TRANS, 'T' ) * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE c END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE c END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = ZERO IF( NOCONJ )THEN DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE ELSE DO 100, I = 1, M TEMP = TEMP + DCONJG( A( I, J ) )*X( I ) 100 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 110 CONTINUE ELSE DO 140, J = 1, N TEMP = ZERO IX = KX IF( NOCONJ )THEN DO 120, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 120 CONTINUE ELSE DO 130, I = 1, M TEMP = TEMP + DCONJG( A( I, J ) )*X( IX ) IX = IX + INCX 130 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 140 CONTINUE END IF END IF * RETURN * * End of ZGEMV . * END SUBROUTINE ZGERC ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. COMPLEX*16 ALPHA INTEGER INCX, INCY, LDA, M, N * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZGERC performs the rank 1 operation * * A := alpha*x*conjg( y' ) + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZGERC ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N c IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( Y( JY ) ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE c END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N c IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( Y( JY ) ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE c END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of ZGERC . * END SUBROUTINE ZHEMV ( UPLO, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. COMPLEX*16 ALPHA, BETA INTEGER INCX, INCY, LDA, N CHARACTER*1 UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZHEMV performs the matrix-vector operation * * y := alpha*A*x + beta*y, * * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n hermitian matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array A is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of A * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of A * is to be referenced. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * 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, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular part of the hermitian matrix and the strictly * lower triangular part of A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular 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 and 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. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. On exit, Y is overwritten by the updated * vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, DBLE * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 5 ELSE IF( INCX.EQ.0 )THEN INFO = 7 ELSE IF( INCY.EQ.0 )THEN INFO = 10 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZHEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set up the start points in X and Y. * IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( N - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( N - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, N Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, N Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, N Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, N Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( UPLO, 'U' ) )THEN * * Form y when A is stored in upper triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO DO 50, I = 1, J - 1 Y( I ) = Y( I ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( I ) 50 CONTINUE Y( J ) = Y( J ) + TEMP1*DBLE( A( J, J ) ) + ALPHA*TEMP2 60 CONTINUE ELSE JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY DO 70, I = 1, J - 1 Y( IY ) = Y( IY ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( IX ) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y( JY ) = Y( JY ) + TEMP1*DBLE( A( J, J ) ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 80 CONTINUE END IF ELSE * * Form y when A is stored in lower triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 100, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO Y( J ) = Y( J ) + TEMP1*DBLE( A( J, J ) ) DO 90, I = J + 1, N Y( I ) = Y( I ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( I ) 90 CONTINUE Y( J ) = Y( J ) + ALPHA*TEMP2 100 CONTINUE ELSE JX = KX JY = KY DO 120, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO Y( JY ) = Y( JY ) + TEMP1*DBLE( A( J, J ) ) IX = JX IY = JY DO 110, I = J + 1, N IX = IX + INCX IY = IY + INCY Y( IY ) = Y( IY ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + DCONJG( A( I, J ) )*X( IX ) 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of ZHEMV . * END SUBROUTINE ZHER2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. COMPLEX*16 ALPHA INTEGER INCX, INCY, LDA, N CHARACTER*1 UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZHER2 performs the hermitian rank 2 operation * * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, * * where alpha is a scalar, x and y are n element vectors and A is an n * by n hermitian matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array A is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of A * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of A * is to be referenced. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular part of the hermitian matrix and the strictly * lower triangular part of A is not referenced. On exit, the * upper triangular part of the array A 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 A must contain the lower * triangular part of the hermitian matrix and the strictly * upper triangular part of A is not referenced. On exit, the * lower triangular part of the array A 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. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, DBLE * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZHER2 ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Set up the start points in X and Y if the increments are not both * unity. * IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( N - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( N - 1 )*INCY END IF JX = KX JY = KY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. * IF( LSAME( UPLO, 'U' ) )THEN * * Form A when A is stored in the upper triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 20, J = 1, N c IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( J ) ) TEMP2 = DCONJG( ALPHA*X( J ) ) DO 10, I = 1, J - 1 A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 10 CONTINUE A( J, J ) = DBLE( A( J, J ) ) + $ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 ) c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF 20 CONTINUE ELSE DO 40, J = 1, N c IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( JY ) ) TEMP2 = DCONJG( ALPHA*X( JX ) ) IX = KX IY = KY DO 30, I = 1, J - 1 A( I, J ) = A( I, J ) + X( IX )*TEMP1 $ + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE A( J, J ) = DBLE( A( J, J ) ) + $ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 ) c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF JX = JX + INCX JY = JY + INCY 40 CONTINUE END IF ELSE * * Form A when A is stored in the lower triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N c IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( J ) ) TEMP2 = DCONJG( ALPHA*X( J ) ) A( J, J ) = DBLE( A( J, J ) ) + $ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 ) DO 50, I = J + 1, N A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 50 CONTINUE c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF 60 CONTINUE ELSE DO 80, J = 1, N c IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( JY ) ) TEMP2 = DCONJG( ALPHA*X( JX ) ) A( J, J ) = DBLE( A( J, J ) ) + $ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 ) IX = JX IY = JY DO 70, I = J + 1, N IX = IX + INCX IY = IY + INCY A( I, J ) = A( I, J ) + X( IX )*TEMP1 $ + Y( IY )*TEMP2 70 CONTINUE c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF JX = JX + INCX JY = JY + INCY 80 CONTINUE END IF END IF * RETURN * * End of ZHER2 . * END SUBROUTINE ZHER2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, BETA, $ C, LDC ) * .. Scalar Arguments .. CHARACTER TRANS, UPLO INTEGER K, LDA, LDB, LDC, N DOUBLE PRECISION BETA COMPLEX*16 ALPHA * .. * .. Array Arguments .. COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * ZHER2K 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*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is * k when TRANS = 'N' or 'n', and is n otherwise. * Before entry with TRANS = 'N' or 'n', the leading n by k * part of the array A must contain the matrix A, otherwise * the leading k by n part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANS = 'N' or 'n' * then LDA must be at least max( 1, n ), otherwise LDA must * be at least max( 1, k ). * Unchanged on exit. * * B - COMPLEX*16 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 - DOUBLE PRECISION . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX*16 array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the hermitian matrix and the strictly * lower triangular part of C is not referenced. On exit, the * upper triangular part of the array C is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array C must contain the lower * triangular part of the hermitian matrix and the strictly * upper triangular part of C is not referenced. On exit, the * lower triangular part of the array C is overwritten by the * lower triangular part of the updated matrix. * Note that the imaginary parts of the diagonal elements need * not be set, they are assumed to be zero, and on exit they * are set to zero. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, n ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. * Ed Anderson, Cray Research Inc. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC DBLE, DCONJG, MAX * .. * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA COMPLEX*16 TEMP1, TEMP2 * .. * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * Test the input parameters. * IF( LSAME( TRANS, 'N' ) ) THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND. $ ( .NOT.LSAME( TRANS, '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( 'ZHER2K', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. $ ( BETA.EQ.ONE ) ) )RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO ) THEN IF( UPPER ) THEN IF( BETA.EQ.DBLE( ZERO ) ) THEN DO 20 J = 1, N DO 10 I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1, N DO 30 I = 1, J - 1 C( I, J ) = BETA*C( I, J ) 30 CONTINUE C( J, J ) = BETA*DBLE( C( J, J ) ) 40 CONTINUE END IF ELSE IF( BETA.EQ.DBLE( ZERO ) ) THEN DO 60 J = 1, N DO 50 I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80 J = 1, N C( J, J ) = BETA*DBLE( C( J, J ) ) DO 70 I = J + 1, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * Start the operations. * IF( LSAME( TRANS, 'N' ) ) THEN * * Form C := alpha*A*conjg( B' ) + conjg( alpha )*B*conjg( A' ) + * C. * IF( UPPER ) THEN DO 130 J = 1, N IF( BETA.EQ.DBLE( ZERO ) ) THEN DO 90 I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE ) THEN DO 100 I = 1, J - 1 C( I, J ) = BETA*C( I, J ) 100 CONTINUE C( J, J ) = BETA*DBLE( C( J, J ) ) ELSE C( J, J ) = DBLE( C( J, J ) ) END IF DO 120 L = 1, K c IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) ) $ THEN TEMP1 = ALPHA*DCONJG( B( J, L ) ) TEMP2 = DCONJG( ALPHA*A( J, L ) ) DO 110 I = 1, J - 1 C( I, J ) = C( I, J ) + A( I, L )*TEMP1 + $ B( I, L )*TEMP2 110 CONTINUE C( J, J ) = DBLE( C( J, J ) ) + $ DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 ) c END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1, N IF( BETA.EQ.DBLE( ZERO ) ) THEN DO 140 I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE ) THEN DO 150 I = J + 1, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE C( J, J ) = BETA*DBLE( C( J, J ) ) ELSE C( J, J ) = DBLE( C( J, J ) ) END IF DO 170 L = 1, K c IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) ) $ THEN TEMP1 = ALPHA*DCONJG( B( J, L ) ) TEMP2 = DCONJG( ALPHA*A( J, L ) ) DO 160 I = J + 1, N C( I, J ) = C( I, J ) + A( I, L )*TEMP1 + $ B( I, L )*TEMP2 160 CONTINUE C( J, J ) = DBLE( C( J, J ) ) + $ DBLE( A( J, L )*TEMP1+B( J, L )*TEMP2 ) c END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*conjg( A' )*B + conjg( alpha )*conjg( B' )*A + * C. * IF( UPPER ) THEN DO 210 J = 1, N DO 200 I = 1, J TEMP1 = ZERO TEMP2 = ZERO DO 190 L = 1, K TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J ) TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J ) 190 CONTINUE IF( I.EQ.J ) THEN IF( BETA.EQ.DBLE( ZERO ) ) THEN C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* $ TEMP2 ) ELSE C( J, J ) = BETA*DBLE( C( J, J ) ) + $ DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* $ TEMP2 ) END IF ELSE IF( BETA.EQ.DBLE( ZERO ) ) THEN C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2 ELSE C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 + $ DCONJG( ALPHA )*TEMP2 END IF END IF 200 CONTINUE 210 CONTINUE ELSE DO 240 J = 1, N DO 230 I = J, N TEMP1 = ZERO TEMP2 = ZERO DO 220 L = 1, K TEMP1 = TEMP1 + DCONJG( A( L, I ) )*B( L, J ) TEMP2 = TEMP2 + DCONJG( B( L, I ) )*A( L, J ) 220 CONTINUE IF( I.EQ.J ) THEN IF( BETA.EQ.DBLE( ZERO ) ) THEN C( J, J ) = DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* $ TEMP2 ) ELSE C( J, J ) = BETA*DBLE( C( J, J ) ) + $ DBLE( ALPHA*TEMP1+DCONJG( ALPHA )* $ TEMP2 ) END IF ELSE IF( BETA.EQ.DBLE( ZERO ) ) THEN C( I, J ) = ALPHA*TEMP1 + DCONJG( ALPHA )*TEMP2 ELSE C( I, J ) = BETA*C( I, J ) + ALPHA*TEMP1 + $ DCONJG( ALPHA )*TEMP2 END IF END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * End of ZHER2K. * END subroutine zscal(n,za,zx,incx) c c scales a vector by a constant. c jack dongarra, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double complex za,zx(*) integer i,incx,ix,n c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 do 10 i = 1,n zx(ix) = za*zx(ix) ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 do 30 i = 1,n zx(i) = za*zx(i) 30 continue return end subroutine zswap (n,zx,incx,zy,incy) c c interchanges two vectors. c jack dongarra, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double complex zx(*),zy(*),ztemp integer i,incx,incy,ix,iy,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n ztemp = zx(ix) zx(ix) = zy(iy) zy(iy) = ztemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 20 do 30 i = 1,n ztemp = zx(i) zx(i) = zy(i) zy(i) = ztemp 30 continue return end SUBROUTINE ZTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB COMPLEX*16 ALPHA * .. Array Arguments .. COMPLEX*16 A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * ZTRMM performs one of the matrix-matrix operations * * B := alpha*op( A )*B, or B := alpha*B*op( A ) * * where alpha is a scalar, B is an m by n matrix, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A' or op( A ) = conjg( A' ). * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) multiplies B from * the left or right as follows: * * SIDE = 'L' or 'l' B := alpha*op( A )*B. * * SIDE = 'R' or 'r' B := alpha*B*op( A ). * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = conjg( A' ). * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - COMPLEX*16 array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the matrix B, and on exit is overwritten by the * transformed matrix. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. Local Scalars .. LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA COMPLEX*16 TEMP * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOCONJ = LSAME( TRANSA, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTRMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*A*B. * IF( UPPER )THEN DO 50, J = 1, N DO 40, K = 1, M c IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) DO 30, I = 1, K - 1 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 30 CONTINUE IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) B( K, J ) = TEMP c END IF 40 CONTINUE 50 CONTINUE ELSE DO 80, J = 1, N DO 70 K = M, 1, -1 c IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) B( K, J ) = TEMP IF( NOUNIT ) $ B( K, J ) = B( K, J )*A( K, K ) DO 60, I = K + 1, M B( I, J ) = B( I, J ) + TEMP*A( I, K ) 60 CONTINUE c END IF 70 CONTINUE 80 CONTINUE END IF ELSE * * Form B := alpha*A'*B or B := alpha*conjg( A' )*B. * IF( UPPER )THEN DO 120, J = 1, N DO 110, I = M, 1, -1 TEMP = B( I, J ) IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 90, K = 1, I - 1 TEMP = TEMP + A( K, I )*B( K, J ) 90 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( I, I ) ) DO 100, K = 1, I - 1 TEMP = TEMP + DCONJG( A( K, I ) )*B( K, J ) 100 CONTINUE END IF B( I, J ) = ALPHA*TEMP 110 CONTINUE 120 CONTINUE ELSE DO 160, J = 1, N DO 150, I = 1, M TEMP = B( I, J ) IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 130, K = I + 1, M TEMP = TEMP + A( K, I )*B( K, J ) 130 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( I, I ) ) DO 140, K = I + 1, M TEMP = TEMP + DCONJG( A( K, I ) )*B( K, J ) 140 CONTINUE END IF B( I, J ) = ALPHA*TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*A. * IF( UPPER )THEN DO 200, J = N, 1, -1 TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 170, I = 1, M B( I, J ) = TEMP*B( I, J ) 170 CONTINUE DO 190, K = 1, J - 1 c IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 180, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 180 CONTINUE c END IF 190 CONTINUE 200 CONTINUE ELSE DO 240, J = 1, N TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 210, I = 1, M B( I, J ) = TEMP*B( I, J ) 210 CONTINUE DO 230, K = J + 1, N c IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 220, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 220 CONTINUE c END IF 230 CONTINUE 240 CONTINUE END IF ELSE * * Form B := alpha*B*A' or B := alpha*B*conjg( A' ). * IF( UPPER )THEN DO 280, K = 1, N DO 260, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN IF( NOCONJ )THEN TEMP = ALPHA*A( J, K ) ELSE TEMP = ALPHA*DCONJG( A( J, K ) ) END IF DO 250, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 250 CONTINUE END IF 260 CONTINUE TEMP = ALPHA IF( NOUNIT )THEN IF( NOCONJ )THEN TEMP = TEMP*A( K, K ) ELSE TEMP = TEMP*DCONJG( A( K, K ) ) END IF END IF IF( TEMP.NE.ONE )THEN DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF 280 CONTINUE ELSE DO 320, K = N, 1, -1 DO 300, J = K + 1, N c IF( A( J, K ).NE.ZERO )THEN IF( NOCONJ )THEN TEMP = ALPHA*A( J, K ) ELSE TEMP = ALPHA*DCONJG( A( J, K ) ) END IF DO 290, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 290 CONTINUE c END IF 300 CONTINUE TEMP = ALPHA IF( NOUNIT )THEN IF( NOCONJ )THEN TEMP = TEMP*A( K, K ) ELSE TEMP = TEMP*DCONJG( A( K, K ) ) END IF END IF IF( TEMP.NE.ONE )THEN DO 310, I = 1, M B( I, K ) = TEMP*B( I, K ) 310 CONTINUE END IF 320 CONTINUE END IF END IF END IF * RETURN * * End of ZTRMM . * END SUBROUTINE ZTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ) * .. * * Purpose * ======= * * ZTRMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, or x := conjg( A' )*x, * * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' x := A*x. * * TRANS = 'T' or 't' x := A'*x. * * TRANS = 'C' or 'c' x := conjg( A' )*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOCONJ, NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTRMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOCONJ = LSAME( TRANS, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := A*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = 1, N c IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*A( I, J ) 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) c END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 30, I = 1, J - 1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) c END IF JX = JX + INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 c IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*A( I, J ) 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) c END IF 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 c IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 70, I = N, J + 1, -1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) c END IF JX = JX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A'*x or x := conjg( A' )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 110, J = N, 1, -1 TEMP = X( J ) IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 90, I = J - 1, 1, -1 TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( J, J ) ) DO 100, I = J - 1, 1, -1 TEMP = TEMP + DCONJG( A( I, J ) )*X( I ) 100 CONTINUE END IF X( J ) = TEMP 110 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 140, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 120, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + A( I, J )*X( IX ) 120 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( J, J ) ) DO 130, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + DCONJG( A( I, J ) )*X( IX ) 130 CONTINUE END IF X( JX ) = TEMP JX = JX - INCX 140 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 170, J = 1, N TEMP = X( J ) IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = J + 1, N TEMP = TEMP + A( I, J )*X( I ) 150 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( J, J ) ) DO 160, I = J + 1, N TEMP = TEMP + DCONJG( A( I, J ) )*X( I ) 160 CONTINUE END IF X( J ) = TEMP 170 CONTINUE ELSE JX = KX DO 200, J = 1, N TEMP = X( JX ) IX = JX IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 180, I = J + 1, N IX = IX + INCX TEMP = TEMP + A( I, J )*X( IX ) 180 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( J, J ) ) DO 190, I = J + 1, N IX = IX + INCX TEMP = TEMP + DCONJG( A( I, J ) )*X( IX ) 190 CONTINUE END IF X( JX ) = TEMP JX = JX + INCX 200 CONTINUE END IF END IF END IF * RETURN * * End of ZTRMV . * END SUBROUTINE ZTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB COMPLEX*16 ALPHA * .. Array Arguments .. COMPLEX*16 A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * ZTRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A' or op( A ) = conjg( A' ). * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = conjg( A' ). * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - COMPLEX*16 array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. Local Scalars .. LOGICAL LSIDE, NOCONJ, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA COMPLEX*16 TEMP * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOCONJ = LSAME( TRANSA, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 c IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE c END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M c IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE c END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B * or B := alpha*inv( conjg( A' ) )*B. * IF( UPPER )THEN DO 140, J = 1, N DO 130, I = 1, M TEMP = ALPHA*B( I, J ) IF( NOCONJ )THEN DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) ELSE DO 120, K = 1, I - 1 TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J ) 120 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( I, I ) ) END IF B( I, J ) = TEMP 130 CONTINUE 140 CONTINUE ELSE DO 180, J = 1, N DO 170, I = M, 1, -1 TEMP = ALPHA*B( I, J ) IF( NOCONJ )THEN DO 150, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) ELSE DO 160, K = I + 1, M TEMP = TEMP - DCONJG( A( K, I ) )*B( K, J ) 160 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( I, I ) ) END IF B( I, J ) = TEMP 170 CONTINUE 180 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 230, J = 1, N IF( ALPHA.NE.ONE )THEN DO 190, I = 1, M B( I, J ) = ALPHA*B( I, J ) 190 CONTINUE END IF DO 210, K = 1, J - 1 c IF( A( K, J ).NE.ZERO )THEN DO 200, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 200 CONTINUE c END IF 210 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 220, I = 1, M B( I, J ) = TEMP*B( I, J ) 220 CONTINUE END IF 230 CONTINUE ELSE DO 280, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 240, I = 1, M B( I, J ) = ALPHA*B( I, J ) 240 CONTINUE END IF DO 260, K = J + 1, N c IF( A( K, J ).NE.ZERO )THEN DO 250, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 250 CONTINUE c END IF 260 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 270, I = 1, M B( I, J ) = TEMP*B( I, J ) 270 CONTINUE END IF 280 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ) * or B := alpha*B*inv( conjg( A' ) ). * IF( UPPER )THEN DO 330, K = N, 1, -1 IF( NOUNIT )THEN IF( NOCONJ )THEN TEMP = ONE/A( K, K ) ELSE TEMP = ONE/DCONJG( A( K, K ) ) END IF DO 290, I = 1, M B( I, K ) = TEMP*B( I, K ) 290 CONTINUE END IF DO 310, J = 1, K - 1 c IF( A( J, K ).NE.ZERO )THEN IF( NOCONJ )THEN TEMP = A( J, K ) ELSE TEMP = DCONJG( A( J, K ) ) c END IF DO 300, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 300 CONTINUE END IF 310 CONTINUE IF( ALPHA.NE.ONE )THEN DO 320, I = 1, M B( I, K ) = ALPHA*B( I, K ) 320 CONTINUE END IF 330 CONTINUE ELSE DO 380, K = 1, N IF( NOUNIT )THEN IF( NOCONJ )THEN TEMP = ONE/A( K, K ) ELSE TEMP = ONE/DCONJG( A( K, K ) ) END IF DO 340, I = 1, M B( I, K ) = TEMP*B( I, K ) 340 CONTINUE END IF DO 360, J = K + 1, N c IF( A( J, K ).NE.ZERO )THEN IF( NOCONJ )THEN TEMP = A( J, K ) ELSE TEMP = DCONJG( A( J, K ) ) END IF DO 350, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 350 CONTINUE c END IF 360 CONTINUE IF( ALPHA.NE.ONE )THEN DO 370, I = 1, M B( I, K ) = ALPHA*B( I, K ) 370 CONTINUE END IF 380 CONTINUE END IF END IF END IF * RETURN * * End of ZTRSM . * END SUBROUTINE ZTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ) * .. * * Purpose * ======= * * ZTRSV solves one of the systems of equations * * A*x = b, or A'*x = b, or conjg( A' )*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' or 'n' A*x = b. * * TRANS = 'T' or 't' A'*x = b. * * TRANS = 'C' or 'c' conjg( A' )*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOCONJ, NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTRSV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOCONJ = LSAME( TRANS, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := inv( A )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 c IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*A( I, J ) 10 CONTINUE c END IF 20 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 c IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 30, I = J - 1, 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 30 CONTINUE c END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N c IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*A( I, J ) 50 CONTINUE c END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N c IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 70, I = J + 1, N IX = IX + INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 70 CONTINUE c END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form x := inv( A' )*x or x := inv( conjg( A' ) )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = X( J ) IF( NOCONJ )THEN DO 90, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) ELSE DO 100, I = 1, J - 1 TEMP = TEMP - DCONJG( A( I, J ) )*X( I ) 100 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( J, J ) ) END IF X( J ) = TEMP 110 CONTINUE ELSE JX = KX DO 140, J = 1, N IX = KX TEMP = X( JX ) IF( NOCONJ )THEN DO 120, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX + INCX 120 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) ELSE DO 130, I = 1, J - 1 TEMP = TEMP - DCONJG( A( I, J ) )*X( IX ) IX = IX + INCX 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( J, J ) ) END IF X( JX ) = TEMP JX = JX + INCX 140 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 170, J = N, 1, -1 TEMP = X( J ) IF( NOCONJ )THEN DO 150, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( I ) 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) ELSE DO 160, I = N, J + 1, -1 TEMP = TEMP - DCONJG( A( I, J ) )*X( I ) 160 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( J, J ) ) END IF X( J ) = TEMP 170 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 200, J = N, 1, -1 IX = KX TEMP = X( JX ) IF( NOCONJ )THEN DO 180, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX - INCX 180 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) ELSE DO 190, I = N, J + 1, -1 TEMP = TEMP - DCONJG( A( I, J ) )*X( IX ) IX = IX - INCX 190 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( J, J ) ) END IF X( JX ) = TEMP JX = JX - INCX 200 CONTINUE END IF END IF END IF * RETURN * * End of ZTRSV . * END subroutine zdrot (n,zx,incx,zy,incy,c,s) c c applies a plane rotation, where the cos and sin (c and s) are c double precision and the vectors zx and zy are double complex. c jack dongarra, linpack, 3/11/78. c double complex zx(1),zy(1),ztemp double precision c,s integer i,incx,incy,ix,iy,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n ztemp = c*zx(ix) + s*zy(iy) zy(iy) = c*zy(iy) - s*zx(ix) zx(ix) = ztemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c 20 do 30 i = 1,n ztemp = c*zx(i) + s*zy(i) zy(i) = c*zy(i) - s*zx(i) zx(i) = ztemp 30 continue return end SUBROUTINE ZGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. COMPLEX*16 ALPHA, BETA INTEGER INCX, INCY, KL, KU, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZGBMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, or * * y := alpha*conjg( A' )*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n band matrix, with kl sub-diagonals and ku super-diagonals. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*conjg( A' )*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * KL - INTEGER. * On entry, KL specifies the number of sub-diagonals of the * matrix A. KL must satisfy 0 .le. KL. * Unchanged on exit. * * KU - INTEGER. * On entry, KU specifies the number of super-diagonals of the * matrix A. KU must satisfy 0 .le. KU. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry, the leading ( kl + ku + 1 ) by n part of the * array A must contain the matrix of coefficients, supplied * column by column, with the leading diagonal of the matrix in * row ( ku + 1 ) of the array, the first super-diagonal * starting at position 2 in row ku, the first sub-diagonal * starting at position 1 in row ( ku + 2 ), and so on. * Elements in the array A that do not correspond to elements * in the band matrix (such as the top left ku by ku triangle) * are not referenced. * The following program segment will transfer a band matrix * from conventional full matrix storage to band storage: * * DO 20, J = 1, N * K = KU + 1 - J * DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL ) * A( K + I, J ) = matrix( I, J ) * 10 CONTINUE * 20 CONTINUE * * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * ( kl + ku + 1 ). * Unchanged on exit. * * X - COMPLEX*16 array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - COMPLEX*16 array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry, the incremented array Y must contain the * vector y. On exit, Y is overwritten by the updated vector y. * * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY, $ LENX, LENY LOGICAL NOCONJ * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( KL.LT.0 )THEN INFO = 4 ELSE IF( KU.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN INFO = 8 ELSE IF( INCX.EQ.0 )THEN INFO = 10 ELSE IF( INCY.EQ.0 )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZGBMV ', 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 * NOCONJ = LSAME( TRANS, 'T' ) * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the band part of A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN KUP1 = KU + 1 IF( LSAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) K = KUP1 - J DO 50, I = MAX( 1, J - KU ), MIN( M, J + KL ) Y( I ) = Y( I ) + TEMP*A( K + I, J ) 50 CONTINUE c END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY K = KUP1 - J DO 70, I = MAX( 1, J - KU ), MIN( M, J + KL ) Y( IY ) = Y( IY ) + TEMP*A( K + I, J ) IY = IY + INCY 70 CONTINUE c END IF JX = JX + INCX IF( J.GT.KU ) $ KY = KY + INCY 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y or y := alpha*conjg( A' )*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = ZERO K = KUP1 - J IF( NOCONJ )THEN DO 90, I = MAX( 1, J - KU ), MIN( M, J + KL ) TEMP = TEMP + A( K + I, J )*X( I ) 90 CONTINUE ELSE DO 100, I = MAX( 1, J - KU ), MIN( M, J + KL ) TEMP = TEMP + DCONJG( A( K + I, J ) )*X( I ) 100 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 110 CONTINUE ELSE DO 140, J = 1, N TEMP = ZERO IX = KX K = KUP1 - J IF( NOCONJ )THEN DO 120, I = MAX( 1, J - KU ), MIN( M, J + KL ) TEMP = TEMP + A( K + I, J )*X( IX ) IX = IX + INCX 120 CONTINUE ELSE DO 130, I = MAX( 1, J - KU ), MIN( M, J + KL ) TEMP = TEMP + DCONJG( A( K + I, J ) )*X( IX ) IX = IX + INCX 130 CONTINUE END IF Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY IF( J.GT.KU ) $ KX = KX + INCX 140 CONTINUE END IF END IF * RETURN * * End of ZGBMV . * END SUBROUTINE ZGERU ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. COMPLEX*16 ALPHA INTEGER INCX, INCY, LDA, M, N * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZGERU performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZGERU ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N c IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE c END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N c IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE c END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of ZGERU . * END SUBROUTINE ZHBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. COMPLEX*16 ALPHA, BETA INTEGER INCX, INCY, K, LDA, N CHARACTER*1 UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZHBMV performs the matrix-vector operation * * y := alpha*A*x + beta*y, * * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n hermitian band matrix, with k super-diagonals. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the band matrix A is being supplied as * follows: * * UPLO = 'U' or 'u' The upper triangular part of A is * being supplied. * * UPLO = 'L' or 'l' The lower triangular part of A is * being supplied. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of super-diagonals of the * matrix A. K must satisfy 0 .le. K. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) * by n part of the array A must contain the upper triangular * band part of the hermitian matrix, supplied column by * column, with the leading diagonal of the matrix in row * ( k + 1 ) of the array, the first super-diagonal starting at * position 2 in row k, and so on. The top left k by k triangle * of the array A is not referenced. * The following program segment will transfer the upper * triangular part of a hermitian band matrix from conventional * full matrix storage to band storage: * * DO 20, J = 1, N * M = K + 1 - J * DO 10, I = MAX( 1, J - K ), J * A( M + I, J ) = matrix( I, J ) * 10 CONTINUE * 20 CONTINUE * * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) * by n part of the array A must contain the lower triangular * band part of the hermitian matrix, supplied column by * column, with the leading diagonal of the matrix in row 1 of * the array, the first sub-diagonal starting at position 1 in * row 2, and so on. The bottom right k by k triangle of the * array A is not referenced. * The following program segment will transfer the lower * triangular part of a hermitian band matrix from conventional * full matrix storage to band storage: * * DO 20, J = 1, N * M = 1 - J * DO 10, I = J, MIN( N, J + K ) * A( M + I, J ) = matrix( I, J ) * 10 CONTINUE * 20 CONTINUE * * Note that the imaginary parts of the diagonal elements need * not be set and 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. LDA must be at least * ( k + 1 ). * Unchanged on exit. * * X - COMPLEX*16 array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * Y - COMPLEX*16 array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the * vector y. On exit, Y is overwritten by the updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN, DBLE * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( K.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZHBMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set up the start points in X and Y. * IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( N - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( N - 1 )*INCY END IF * * Start the operations. In this version the elements of the array A * are accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, N Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, N Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, N Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, N Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LSAME( UPLO, 'U' ) )THEN * * Form y when upper triangle of A is stored. * KPLUS1 = K + 1 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO L = KPLUS1 - J DO 50, I = MAX( 1, J - K ), J - 1 Y( I ) = Y( I ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( I ) 50 CONTINUE Y( J ) = Y( J ) + TEMP1*DBLE( A( KPLUS1, J ) ) $ + ALPHA*TEMP2 60 CONTINUE ELSE JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY L = KPLUS1 - J DO 70, I = MAX( 1, J - K ), J - 1 Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( IX ) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y( JY ) = Y( JY ) + TEMP1*DBLE( A( KPLUS1, J ) ) $ + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY IF( J.GT.K )THEN KX = KX + INCX KY = KY + INCY END IF 80 CONTINUE END IF ELSE * * Form y when lower triangle of A is stored. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 100, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO Y( J ) = Y( J ) + TEMP1*DBLE( A( 1, J ) ) L = 1 - J DO 90, I = J + 1, MIN( N, J + K ) Y( I ) = Y( I ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( I ) 90 CONTINUE Y( J ) = Y( J ) + ALPHA*TEMP2 100 CONTINUE ELSE JX = KX JY = KY DO 120, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO Y( JY ) = Y( JY ) + TEMP1*DBLE( A( 1, J ) ) L = 1 - J IX = JX IY = JY DO 110, I = J + 1, MIN( N, J + K ) IX = IX + INCX IY = IY + INCY Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + DCONJG( A( L + I, J ) )*X( IX ) 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of ZHBMV . * END SUBROUTINE ZHEMM ( 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*16 ALPHA, BETA * .. Array Arguments .. COMPLEX*16 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. * * 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*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 DCONJG, MAX, DBLE * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, K, NROWA COMPLEX*16 TEMP1, TEMP2 * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * 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 SUBROUTINE ZHER ( UPLO, N, ALPHA, X, INCX, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, LDA, N CHARACTER*1 UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ) * .. * * Purpose * ======= * * ZHER performs the hermitian rank 1 operation * * A := alpha*x*conjg( x' ) + A, * * where alpha is a real scalar, x is an n element vector and A is an * n by n hermitian matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array A is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of A * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of A * is to be referenced. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular part of the hermitian matrix and the strictly * lower triangular part of A is not referenced. On exit, the * upper triangular part of the array A 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 A must contain the lower * triangular part of the hermitian matrix and the strictly * upper triangular part of A is not referenced. On exit, the * lower triangular part of the array A 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. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, KX * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, DBLE * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZHER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ALPHA.EQ.DBLE( ZERO ) ) ) $ RETURN * * Set the start point in X if the increment is not unity. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. * IF( LSAME( UPLO, 'U' ) )THEN * * Form A when A is stored in upper triangle. * IF( INCX.EQ.1 )THEN DO 20, J = 1, N c IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( J ) ) DO 10, I = 1, J - 1 A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( J )*TEMP ) c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( JX ) ) IX = KX DO 30, I = 1, J - 1 A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( JX )*TEMP ) c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF JX = JX + INCX 40 CONTINUE END IF ELSE * * Form A when A is stored in lower triangle. * IF( INCX.EQ.1 )THEN DO 60, J = 1, N c IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( J ) ) A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( J ) ) DO 50, I = J + 1, N A( I, J ) = A( I, J ) + X( I )*TEMP 50 CONTINUE c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( JX ) ) A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( JX ) ) IX = JX DO 70, I = J + 1, N IX = IX + INCX A( I, J ) = A( I, J ) + X( IX )*TEMP 70 CONTINUE c ELSE c A( J, J ) = DBLE( A( J, J ) ) c END IF JX = JX + INCX 80 CONTINUE END IF END IF * RETURN * * End of ZHER . * END SUBROUTINE ZHERK( UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER TRANS, UPLO INTEGER K, LDA, LDC, N DOUBLE PRECISION ALPHA, BETA * .. * .. Array Arguments .. COMPLEX*16 A( LDA, * ), C( LDC, * ) * .. * * Purpose * ======= * * ZHERK performs one of the hermitian rank k operations * * C := alpha*A*conjg( A' ) + beta*C, * * or * * C := alpha*conjg( A' )*A + beta*C, * * where alpha and beta are real scalars, C is an n by n hermitian * matrix and A is an n by k matrix in the first case and a k by n * matrix in the second case. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array C is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of C * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of C * is to be referenced. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' C := alpha*A*conjg( A' ) + beta*C. * * TRANS = 'C' or 'c' C := alpha*conjg( A' )*A + beta*C. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with TRANS = 'N' or 'n', K specifies the number * of columns of the matrix A, and on entry with * TRANS = 'C' or 'c', K specifies the number of rows of the * matrix A. K must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is * k when TRANS = 'N' or 'n', and is n otherwise. * Before entry with TRANS = 'N' or 'n', the leading n by k * part of the array A must contain the matrix A, otherwise * the leading k by n part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANS = 'N' or 'n' * then LDA must be at least max( 1, n ), otherwise LDA must * be at least max( 1, k ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX*16 array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the hermitian matrix and the strictly * lower triangular part of C is not referenced. On exit, the * upper triangular part of the array C is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array C must contain the lower * triangular part of the hermitian matrix and the strictly * upper triangular part of C is not referenced. On exit, the * lower triangular part of the array C is overwritten by the * lower triangular part of the updated matrix. * Note that the imaginary parts of the diagonal elements need * not be set, they are assumed to be zero, and on exit they * are set to zero. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, n ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * -- Modified 8-Nov-93 to set C(J,J) to DBLE( C(J,J) ) when BETA = 1. * Ed Anderson, Cray Research Inc. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC DBLE, DCMPLX, DCONJG, MAX * .. * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION RTEMP COMPLEX*16 TEMP * .. * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * IF( LSAME( TRANS, 'N' ) ) THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ) .AND. ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ) .AND. $ ( .NOT.LSAME( TRANS, '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( 'ZHERK ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ) .OR. ( ( ( ALPHA.EQ.ZERO ) .OR. ( K.EQ.0 ) ) .AND. $ ( BETA.EQ.ONE ) ) )RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO ) THEN IF( UPPER ) THEN IF( BETA.EQ.ZERO ) THEN DO 20 J = 1, N DO 10 I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40 J = 1, N DO 30 I = 1, J - 1 C( I, J ) = BETA*C( I, J ) 30 CONTINUE C( J, J ) = BETA*DBLE( C( J, J ) ) 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO ) THEN DO 60 J = 1, N DO 50 I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80 J = 1, N C( J, J ) = BETA*DBLE( C( J, J ) ) DO 70 I = J + 1, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * Start the operations. * IF( LSAME( TRANS, 'N' ) ) THEN * * Form C := alpha*A*conjg( A' ) + beta*C. * IF( UPPER ) THEN DO 130 J = 1, N IF( BETA.EQ.ZERO ) THEN DO 90 I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE ) THEN DO 100 I = 1, J - 1 C( I, J ) = BETA*C( I, J ) 100 CONTINUE C( J, J ) = BETA*DBLE( C( J, J ) ) ELSE C( J, J ) = DBLE( C( J, J ) ) END IF DO 120 L = 1, K IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN TEMP = ALPHA*DCONJG( A( J, L ) ) DO 110 I = 1, J - 1 C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE C( J, J ) = DBLE( C( J, J ) ) + $ DBLE( TEMP*A( I, L ) ) END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1, N IF( BETA.EQ.ZERO ) THEN DO 140 I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE ) THEN C( J, J ) = BETA*DBLE( C( J, J ) ) DO 150 I = J + 1, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE ELSE C( J, J ) = DBLE( C( J, J ) ) END IF DO 170 L = 1, K IF( A( J, L ).NE.DCMPLX( ZERO ) ) THEN TEMP = ALPHA*DCONJG( A( J, L ) ) C( J, J ) = DBLE( C( J, J ) ) + $ DBLE( TEMP*A( J, L ) ) DO 160 I = J + 1, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*conjg( A' )*A + beta*C. * IF( UPPER ) THEN DO 220 J = 1, N DO 200 I = 1, J - 1 TEMP = ZERO DO 190 L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO ) THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE RTEMP = ZERO DO 210 L = 1, K RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J ) 210 CONTINUE IF( BETA.EQ.ZERO ) THEN C( J, J ) = ALPHA*RTEMP ELSE C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) ) END IF 220 CONTINUE ELSE DO 260 J = 1, N RTEMP = ZERO DO 230 L = 1, K RTEMP = RTEMP + DCONJG( A( L, J ) )*A( L, J ) 230 CONTINUE IF( BETA.EQ.ZERO ) THEN C( J, J ) = ALPHA*RTEMP ELSE C( J, J ) = ALPHA*RTEMP + BETA*DBLE( C( J, J ) ) END IF DO 250 I = J + 1, N TEMP = ZERO DO 240 L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*A( L, J ) 240 CONTINUE IF( BETA.EQ.ZERO ) THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 250 CONTINUE 260 CONTINUE END IF END IF * RETURN * * End of ZHERK . * END SUBROUTINE ZHPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) * .. Scalar Arguments .. COMPLEX*16 ALPHA, BETA INTEGER INCX, INCY, N CHARACTER*1 UPLO * .. Array Arguments .. COMPLEX*16 AP( * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZHPMV performs the matrix-vector operation * * y := alpha*A*x + beta*y, * * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n hermitian matrix, supplied in packed form. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' or 'u' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' or 'l' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * AP - COMPLEX*16 array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular part of the hermitian matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular part of the hermitian matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. * Note that the imaginary parts of the diagonal elements need * not be set and are assumed to be zero. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. On exit, Y is overwritten by the updated * vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, DBLE * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 6 ELSE IF( INCY.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZHPMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set up the start points in X and Y. * IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( N - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( N - 1 )*INCY END IF * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, N Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, N Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, N Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, N Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN KK = 1 IF( LSAME( UPLO, 'U' ) )THEN * * Form y when AP contains the upper triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO K = KK DO 50, I = 1, J - 1 Y( I ) = Y( I ) + TEMP1*AP( K ) TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( I ) K = K + 1 50 CONTINUE Y( J ) = Y( J ) + TEMP1*DBLE( AP( KK + J - 1 ) ) $ + ALPHA*TEMP2 KK = KK + J 60 CONTINUE ELSE JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY DO 70, K = KK, KK + J - 2 Y( IY ) = Y( IY ) + TEMP1*AP( K ) TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( IX ) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y( JY ) = Y( JY ) + TEMP1*DBLE( AP( KK + J - 1 ) ) $ + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY KK = KK + J 80 CONTINUE END IF ELSE * * Form y when AP contains the lower triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 100, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO Y( J ) = Y( J ) + TEMP1*DBLE( AP( KK ) ) K = KK + 1 DO 90, I = J + 1, N Y( I ) = Y( I ) + TEMP1*AP( K ) TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( I ) K = K + 1 90 CONTINUE Y( J ) = Y( J ) + ALPHA*TEMP2 KK = KK + ( N - J + 1 ) 100 CONTINUE ELSE JX = KX JY = KY DO 120, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO Y( JY ) = Y( JY ) + TEMP1*DBLE( AP( KK ) ) IX = JX IY = JY DO 110, K = KK + 1, KK + N - J IX = IX + INCX IY = IY + INCY Y( IY ) = Y( IY ) + TEMP1*AP( K ) TEMP2 = TEMP2 + DCONJG( AP( K ) )*X( IX ) 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY KK = KK + ( N - J + 1 ) 120 CONTINUE END IF END IF * RETURN * * End of ZHPMV . * END SUBROUTINE ZHPR ( UPLO, N, ALPHA, X, INCX, AP ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, N CHARACTER*1 UPLO * .. Array Arguments .. COMPLEX*16 AP( * ), X( * ) * .. * * Purpose * ======= * * ZHPR performs the hermitian rank 1 operation * * A := alpha*x*conjg( x' ) + A, * * where alpha is a real scalar, x is an n element vector and A is an * n by n hermitian matrix, supplied in packed form. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' or 'u' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' or 'l' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * AP - COMPLEX*16 array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular part of the hermitian matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. On exit, the array * AP is overwritten by the upper triangular part of the * updated matrix. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular part of the hermitian matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. On exit, the array * AP 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. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, DBLE * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZHPR ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ALPHA.EQ.DBLE( ZERO ) ) ) $ RETURN * * Set the start point in X if the increment is not unity. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * KK = 1 IF( LSAME( UPLO, 'U' ) )THEN * * Form A when upper triangle is stored in AP. * IF( INCX.EQ.1 )THEN DO 20, J = 1, N c IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( J ) ) K = KK DO 10, I = 1, J - 1 AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 10 CONTINUE AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) $ + DBLE( X( J )*TEMP ) c ELSE c AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) c END IF KK = KK + J 20 CONTINUE ELSE JX = KX DO 40, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( JX ) ) IX = KX DO 30, K = KK, KK + J - 2 AP( K ) = AP( K ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) $ + DBLE( X( JX )*TEMP ) c ELSE c AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) c END IF JX = JX + INCX KK = KK + J 40 CONTINUE END IF ELSE * * Form A when lower triangle is stored in AP. * IF( INCX.EQ.1 )THEN DO 60, J = 1, N c IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( J ) ) AP( KK ) = DBLE( AP( KK ) ) + DBLE( TEMP*X( J ) ) K = KK + 1 DO 50, I = J + 1, N AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 50 CONTINUE c ELSE c AP( KK ) = DBLE( AP( KK ) ) c END IF KK = KK + N - J + 1 60 CONTINUE ELSE JX = KX DO 80, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( X( JX ) ) AP( KK ) = DBLE( AP( KK ) ) + DBLE( TEMP*X( JX ) ) IX = JX DO 70, K = KK + 1, KK + N - J IX = IX + INCX AP( K ) = AP( K ) + X( IX )*TEMP 70 CONTINUE c ELSE c AP( KK ) = DBLE( AP( KK ) ) c END IF JX = JX + INCX KK = KK + N - J + 1 80 CONTINUE END IF END IF * RETURN * * End of ZHPR . * END SUBROUTINE ZHPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) * .. Scalar Arguments .. COMPLEX*16 ALPHA INTEGER INCX, INCY, N CHARACTER*1 UPLO * .. Array Arguments .. COMPLEX*16 AP( * ), X( * ), Y( * ) * .. * * Purpose * ======= * * ZHPR2 performs the hermitian rank 2 operation * * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, * * where alpha is a scalar, x and y are n element vectors and A is an * n by n hermitian matrix, supplied in packed form. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' or 'u' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' or 'l' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * AP - COMPLEX*16 array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular part of the hermitian matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. On exit, the array * AP is overwritten by the upper triangular part of the * updated matrix. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular part of the hermitian matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. On exit, the array * AP 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. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, DBLE * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZHPR2 ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Set up the start points in X and Y if the increments are not both * unity. * IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( N - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( N - 1 )*INCY END IF JX = KX JY = KY END IF * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * KK = 1 IF( LSAME( UPLO, 'U' ) )THEN * * Form A when upper triangle is stored in AP. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 20, J = 1, N c IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( J ) ) TEMP2 = DCONJG( ALPHA*X( J ) ) K = KK DO 10, I = 1, J - 1 AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 10 CONTINUE AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) + $ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 ) c ELSE c AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) c END IF KK = KK + J 20 CONTINUE ELSE DO 40, J = 1, N c IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( JY ) ) TEMP2 = DCONJG( ALPHA*X( JX ) ) IX = KX IY = KY DO 30, K = KK, KK + J - 2 AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) + $ DBLE( X( JX )*TEMP1 + $ Y( JY )*TEMP2 ) c ELSE c AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) ) c END IF JX = JX + INCX JY = JY + INCY KK = KK + J 40 CONTINUE END IF ELSE * * Form A when lower triangle is stored in AP. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N c IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( J ) ) TEMP2 = DCONJG( ALPHA*X( J ) ) AP( KK ) = DBLE( AP( KK ) ) + $ DBLE( X( J )*TEMP1 + Y( J )*TEMP2 ) K = KK + 1 DO 50, I = J + 1, N AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 50 CONTINUE c ELSE c AP( KK ) = DBLE( AP( KK ) ) c END IF KK = KK + N - J + 1 60 CONTINUE ELSE DO 80, J = 1, N c IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN TEMP1 = ALPHA*DCONJG( Y( JY ) ) TEMP2 = DCONJG( ALPHA*X( JX ) ) AP( KK ) = DBLE( AP( KK ) ) + $ DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 ) IX = JX IY = JY DO 70, K = KK + 1, KK + N - J IX = IX + INCX IY = IY + INCY AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 70 CONTINUE c ELSE c AP( KK ) = DBLE( AP( KK ) ) c END IF JX = JX + INCX JY = JY + INCY KK = KK + N - J + 1 80 CONTINUE END IF END IF * RETURN * * End of ZHPR2 . * END subroutine zrotg(ca,cb,c,s) double complex ca,cb,s double precision c double precision norm,scale double complex alpha intrinsic dconjg, dcmplx if (cdabs(ca) .ne. 0.0d0) go to 10 c = 0.0d0 s = (1.0d0,0.0d0) ca = cb go to 20 10 continue scale = cdabs(ca) + cdabs(cb) norm = scale*dsqrt((cdabs(ca/dcmplx(scale,0.0d0)))**2 + * (cdabs(cb/dcmplx(scale,0.0d0)))**2) alpha = ca /cdabs(ca) c = cdabs(ca) / norm s = alpha * dconjg(cb) / norm ca = alpha * norm 20 continue return end SUBROUTINE ZSYMM ( 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*16 ALPHA, BETA * .. Array Arguments .. COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * ZSYMM 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*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 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*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 MAX * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, K, NROWA COMPLEX*16 TEMP1, TEMP2 * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * 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( 'ZSYMM ', 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 ZSYMM . * END SUBROUTINE ZSYR2K( 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*16 ALPHA, BETA * .. Array Arguments .. COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * ZSYR2K 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*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is * k when TRANS = 'N' or 'n', and is n otherwise. * Before entry with TRANS = 'N' or 'n', the leading n by k * part of the array A must contain the matrix A, otherwise * the leading k by n part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANS = 'N' or 'n' * then LDA must be at least max( 1, n ), otherwise LDA must * be at least max( 1, k ). * Unchanged on exit. * * B - COMPLEX*16 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*16 . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX*16 array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of C is not referenced. On exit, the * upper triangular part of the array C is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array C must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of C is not referenced. On exit, the * lower triangular part of the array C is overwritten by the * lower triangular part of the updated matrix. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, n ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA COMPLEX*16 TEMP1, TEMP2 * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * Test the input parameters. * IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. $ ( .NOT.LSAME( TRANS, 'T' ) ) )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( 'ZSYR2K', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * Start the operations. * IF( LSAME( TRANS, 'N' ) )THEN * * Form C := alpha*A*B' + alpha*B*A' + C. * IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K c IF( ( A( J, L ).NE.ZERO ).OR. c $ ( B( J, L ).NE.ZERO ) )THEN TEMP1 = ALPHA*B( J, L ) TEMP2 = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + A( I, L )*TEMP1 + $ B( I, L )*TEMP2 110 CONTINUE c END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K c IF( ( A( J, L ).NE.ZERO ).OR. c $ ( B( J, L ).NE.ZERO ) )THEN TEMP1 = ALPHA*B( J, L ) TEMP2 = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + A( I, L )*TEMP1 + $ B( I, L )*TEMP2 160 CONTINUE c END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*A'*B + alpha*B'*A + C. * IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP1 = ZERO TEMP2 = ZERO DO 190, L = 1, K TEMP1 = TEMP1 + A( L, I )*B( L, J ) TEMP2 = TEMP2 + B( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ ALPHA*TEMP1 + ALPHA*TEMP2 END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP1 = ZERO TEMP2 = ZERO DO 220, L = 1, K TEMP1 = TEMP1 + A( L, I )*B( L, J ) TEMP2 = TEMP2 + B( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ ALPHA*TEMP1 + ALPHA*TEMP2 END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * End of ZSYR2K. * END SUBROUTINE ZSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC COMPLEX*16 ALPHA, BETA * .. Array Arguments .. COMPLEX*16 A( LDA, * ), C( LDC, * ) * .. * * Purpose * ======= * * ZSYRK 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*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is * k when TRANS = 'N' or 'n', and is n otherwise. * Before entry with TRANS = 'N' or 'n', the leading n by k * part of the array A must contain the matrix A, otherwise * the leading k by n part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANS = 'N' or 'n' * then LDA must be at least max( 1, n ), otherwise LDA must * be at least max( 1, k ). * Unchanged on exit. * * BETA - COMPLEX*16 . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - COMPLEX*16 array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of C is not referenced. On exit, the * upper triangular part of the array C is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array C must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of C is not referenced. On exit, the * lower triangular part of the array C is overwritten by the * lower triangular part of the updated matrix. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, n ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA COMPLEX*16 TEMP * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * Test the input parameters. * IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. $ ( .NOT.LSAME( TRANS, 'T' ) ) )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( 'ZSYRK ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * Start the operations. * IF( LSAME( TRANS, 'N' ) )THEN * * Form C := alpha*A*A' + beta*C. * IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K c IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE c END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K c IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE c END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*A'*A + beta*C. * IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP = ZERO DO 190, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP = ZERO DO 220, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * End of ZSYRK . * END SUBROUTINE ZTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, K, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ) * .. * * Purpose * ======= * * ZTBMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, or x := conjg( A' )*x, * * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular band matrix, with ( k + 1 ) diagonals. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' x := A*x. * * TRANS = 'T' or 't' x := A'*x. * * TRANS = 'C' or 'c' x := conjg( A' )*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with UPLO = 'U' or 'u', K specifies the number of * super-diagonals of the matrix A. * On entry with UPLO = 'L' or 'l', K specifies the number of * sub-diagonals of the matrix A. * K must satisfy 0 .le. K. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) * by n part of the array A must contain the upper triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row * ( k + 1 ) of the array, the first super-diagonal starting at * position 2 in row k, and so on. The top left k by k triangle * of the array A is not referenced. * The following program segment will transfer an upper * triangular band matrix from conventional full matrix storage * to band storage: * * DO 20, J = 1, N * M = K + 1 - J * DO 10, I = MAX( 1, J - K ), J * A( M + I, J ) = matrix( I, J ) * 10 CONTINUE * 20 CONTINUE * * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) * by n part of the array A must contain the lower triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row 1 of * the array, the first sub-diagonal starting at position 1 in * row 2, and so on. The bottom right k by k triangle of the * array A is not referenced. * The following program segment will transfer a lower * triangular band matrix from conventional full matrix storage * to band storage: * * DO 20, J = 1, N * M = 1 - J * DO 10, I = J, MIN( N, J + K ) * A( M + I, J ) = matrix( I, J ) * 10 CONTINUE * 20 CONTINUE * * Note that when DIAG = 'U' or 'u' the elements of the array A * corresponding to the diagonal elements of the matrix are not * referenced, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * ( k + 1 ). * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L LOGICAL NOCONJ, NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( K.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 7 ELSE IF( INCX.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTBMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOCONJ = LSAME( TRANS, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := A*x. * IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 20, J = 1, N c IF( X( J ).NE.ZERO )THEN TEMP = X( J ) L = KPLUS1 - J DO 10, I = MAX( 1, J - K ), J - 1 X( I ) = X( I ) + TEMP*A( L + I, J ) 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( KPLUS1, J ) c END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX L = KPLUS1 - J DO 30, I = MAX( 1, J - K ), J - 1 X( IX ) = X( IX ) + TEMP*A( L + I, J ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( KPLUS1, J ) c END IF JX = JX + INCX IF( J.GT.K ) $ KX = KX + INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 c IF( X( J ).NE.ZERO )THEN TEMP = X( J ) L = 1 - J DO 50, I = MIN( N, J + K ), J + 1, -1 X( I ) = X( I ) + TEMP*A( L + I, J ) 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( 1, J ) c END IF 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 c IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX L = 1 - J DO 70, I = MIN( N, J + K ), J + 1, -1 X( IX ) = X( IX ) + TEMP*A( L + I, J ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( 1, J ) c END IF JX = JX - INCX IF( ( N - J ).GE.K ) $ KX = KX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A'*x or x := conjg( A' )*x. * IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 110, J = N, 1, -1 TEMP = X( J ) L = KPLUS1 - J IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( KPLUS1, J ) DO 90, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + A( L + I, J )*X( I ) 90 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( KPLUS1, J ) ) DO 100, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + DCONJG( A( L + I, J ) )*X( I ) 100 CONTINUE END IF X( J ) = TEMP 110 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 140, J = N, 1, -1 TEMP = X( JX ) KX = KX - INCX IX = KX L = KPLUS1 - J IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( KPLUS1, J ) DO 120, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + A( L + I, J )*X( IX ) IX = IX - INCX 120 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( KPLUS1, J ) ) DO 130, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + DCONJG( A( L + I, J ) )*X( IX ) IX = IX - INCX 130 CONTINUE END IF X( JX ) = TEMP JX = JX - INCX 140 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 170, J = 1, N TEMP = X( J ) L = 1 - J IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( 1, J ) DO 150, I = J + 1, MIN( N, J + K ) TEMP = TEMP + A( L + I, J )*X( I ) 150 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( 1, J ) ) DO 160, I = J + 1, MIN( N, J + K ) TEMP = TEMP + DCONJG( A( L + I, J ) )*X( I ) 160 CONTINUE END IF X( J ) = TEMP 170 CONTINUE ELSE JX = KX DO 200, J = 1, N TEMP = X( JX ) KX = KX + INCX IX = KX L = 1 - J IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*A( 1, J ) DO 180, I = J + 1, MIN( N, J + K ) TEMP = TEMP + A( L + I, J )*X( IX ) IX = IX + INCX 180 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( A( 1, J ) ) DO 190, I = J + 1, MIN( N, J + K ) TEMP = TEMP + DCONJG( A( L + I, J ) )*X( IX ) IX = IX + INCX 190 CONTINUE END IF X( JX ) = TEMP JX = JX + INCX 200 CONTINUE END IF END IF END IF * RETURN * * End of ZTBMV . * END SUBROUTINE ZTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, K, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. COMPLEX*16 A( LDA, * ), X( * ) * .. * * Purpose * ======= * * ZTBSV solves one of the systems of equations * * A*x = b, or A'*x = b, or conjg( A' )*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular band matrix, with ( k + 1 ) * diagonals. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' or 'n' A*x = b. * * TRANS = 'T' or 't' A'*x = b. * * TRANS = 'C' or 'c' conjg( A' )*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with UPLO = 'U' or 'u', K specifies the number of * super-diagonals of the matrix A. * On entry with UPLO = 'L' or 'l', K specifies the number of * sub-diagonals of the matrix A. * K must satisfy 0 .le. K. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) * by n part of the array A must contain the upper triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row * ( k + 1 ) of the array, the first super-diagonal starting at * position 2 in row k, and so on. The top left k by k triangle * of the array A is not referenced. * The following program segment will transfer an upper * triangular band matrix from conventional full matrix storage * to band storage: * * DO 20, J = 1, N * M = K + 1 - J * DO 10, I = MAX( 1, J - K ), J * A( M + I, J ) = matrix( I, J ) * 10 CONTINUE * 20 CONTINUE * * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) * by n part of the array A must contain the lower triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row 1 of * the array, the first sub-diagonal starting at position 1 in * row 2, and so on. The bottom right k by k triangle of the * array A is not referenced. * The following program segment will transfer a lower * triangular band matrix from conventional full matrix storage * to band storage: * * DO 20, J = 1, N * M = 1 - J * DO 10, I = J, MIN( N, J + K ) * A( M + I, J ) = matrix( I, J ) * 10 CONTINUE * 20 CONTINUE * * Note that when DIAG = 'U' or 'u' the elements of the array A * corresponding to the diagonal elements of the matrix are not * referenced, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * ( k + 1 ). * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L LOGICAL NOCONJ, NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( K.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 7 ELSE IF( INCX.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTBSV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOCONJ = LSAME( TRANS, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed by sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := inv( A )*x. * IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 c IF( X( J ).NE.ZERO )THEN L = KPLUS1 - J IF( NOUNIT ) $ X( J ) = X( J )/A( KPLUS1, J ) TEMP = X( J ) DO 10, I = J - 1, MAX( 1, J - K ), -1 X( I ) = X( I ) - TEMP*A( L + I, J ) 10 CONTINUE c END IF 20 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 40, J = N, 1, -1 KX = KX - INCX c IF( X( JX ).NE.ZERO )THEN IX = KX L = KPLUS1 - J IF( NOUNIT ) $ X( JX ) = X( JX )/A( KPLUS1, J ) TEMP = X( JX ) DO 30, I = J - 1, MAX( 1, J - K ), -1 X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX - INCX 30 CONTINUE c END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N c IF( X( J ).NE.ZERO )THEN L = 1 - J IF( NOUNIT ) $ X( J ) = X( J )/A( 1, J ) TEMP = X( J ) DO 50, I = J + 1, MIN( N, J + K ) X( I ) = X( I ) - TEMP*A( L + I, J ) 50 CONTINUE c END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N KX = KX + INCX c IF( X( JX ).NE.ZERO )THEN IX = KX L = 1 - J IF( NOUNIT ) $ X( JX ) = X( JX )/A( 1, J ) TEMP = X( JX ) DO 70, I = J + 1, MIN( N, J + K ) X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX + INCX 70 CONTINUE c END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form x := inv( A' )*x or x := inv( conjg( A') )*x. * IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = X( J ) L = KPLUS1 - J IF( NOCONJ )THEN DO 90, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( KPLUS1, J ) ELSE DO 100, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I ) 100 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( KPLUS1, J ) ) END IF X( J ) = TEMP 110 CONTINUE ELSE JX = KX DO 140, J = 1, N TEMP = X( JX ) IX = KX L = KPLUS1 - J IF( NOCONJ )THEN DO 120, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX + INCX 120 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( KPLUS1, J ) ELSE DO 130, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX ) IX = IX + INCX 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( KPLUS1, J ) ) END IF X( JX ) = TEMP JX = JX + INCX IF( J.GT.K ) $ KX = KX + INCX 140 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 170, J = N, 1, -1 TEMP = X( J ) L = 1 - J IF( NOCONJ )THEN DO 150, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( I ) 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( 1, J ) ELSE DO 160, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( I ) 160 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( 1, J ) ) END IF X( J ) = TEMP 170 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 200, J = N, 1, -1 TEMP = X( JX ) IX = KX L = 1 - J IF( NOCONJ )THEN DO 180, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX - INCX 180 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( 1, J ) ELSE DO 190, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - DCONJG( A( L + I, J ) )*X( IX ) IX = IX - INCX 190 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( A( 1, J ) ) END IF X( JX ) = TEMP JX = JX - INCX IF( ( N - J ).GE.K ) $ KX = KX - INCX 200 CONTINUE END IF END IF END IF * RETURN * * End of ZTBSV . * END SUBROUTINE ZTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. COMPLEX*16 AP( * ), X( * ) * .. * * Purpose * ======= * * ZTPMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, or x := conjg( A' )*x, * * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix, supplied in packed form. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' x := A*x. * * TRANS = 'T' or 't' x := A'*x. * * TRANS = 'C' or 'c' x := conjg( A' )*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * AP - COMPLEX*16 array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) * respectively, and so on. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) * respectively, and so on. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced, but are assumed to be unity. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX LOGICAL NOCONJ, NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( INCX.EQ.0 )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTPMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOCONJ = LSAME( TRANS, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of AP are * accessed sequentially with one pass through AP. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x:= A*x. * IF( LSAME( UPLO, 'U' ) )THEN KK = 1 IF( INCX.EQ.1 )THEN DO 20, J = 1, N c IF( X( J ).NE.ZERO )THEN TEMP = X( J ) K = KK DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*AP( K ) K = K + 1 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*AP( KK + J - 1 ) c END IF KK = KK + J 20 CONTINUE ELSE JX = KX DO 40, J = 1, N c IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 30, K = KK, KK + J - 2 X( IX ) = X( IX ) + TEMP*AP( K ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*AP( KK + J - 1 ) c END IF JX = JX + INCX KK = KK + J 40 CONTINUE END IF ELSE KK = ( N*( N + 1 ) )/2 IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 c IF( X( J ).NE.ZERO )THEN TEMP = X( J ) K = KK DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*AP( K ) K = K - 1 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*AP( KK - N + J ) c END IF KK = KK - ( N - J + 1 ) 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 c IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 70, K = KK, KK - ( N - ( J + 1 ) ), -1 X( IX ) = X( IX ) + TEMP*AP( K ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*AP( KK - N + J ) c END IF JX = JX - INCX KK = KK - ( N - J + 1 ) 80 CONTINUE END IF END IF ELSE * * Form x := A'*x or x := conjg( A' )*x. * IF( LSAME( UPLO, 'U' ) )THEN KK = ( N*( N + 1 ) )/2 IF( INCX.EQ.1 )THEN DO 110, J = N, 1, -1 TEMP = X( J ) K = KK - 1 IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) DO 90, I = J - 1, 1, -1 TEMP = TEMP + AP( K )*X( I ) K = K - 1 90 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( AP( KK ) ) DO 100, I = J - 1, 1, -1 TEMP = TEMP + DCONJG( AP( K ) )*X( I ) K = K - 1 100 CONTINUE END IF X( J ) = TEMP KK = KK - J 110 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 140, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) DO 120, K = KK - 1, KK - J + 1, -1 IX = IX - INCX TEMP = TEMP + AP( K )*X( IX ) 120 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( AP( KK ) ) DO 130, K = KK - 1, KK - J + 1, -1 IX = IX - INCX TEMP = TEMP + DCONJG( AP( K ) )*X( IX ) 130 CONTINUE END IF X( JX ) = TEMP JX = JX - INCX KK = KK - J 140 CONTINUE END IF ELSE KK = 1 IF( INCX.EQ.1 )THEN DO 170, J = 1, N TEMP = X( J ) K = KK + 1 IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) DO 150, I = J + 1, N TEMP = TEMP + AP( K )*X( I ) K = K + 1 150 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( AP( KK ) ) DO 160, I = J + 1, N TEMP = TEMP + DCONJG( AP( K ) )*X( I ) K = K + 1 160 CONTINUE END IF X( J ) = TEMP KK = KK + ( N - J + 1 ) 170 CONTINUE ELSE JX = KX DO 200, J = 1, N TEMP = X( JX ) IX = JX IF( NOCONJ )THEN IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) DO 180, K = KK + 1, KK + N - J IX = IX + INCX TEMP = TEMP + AP( K )*X( IX ) 180 CONTINUE ELSE IF( NOUNIT ) $ TEMP = TEMP*DCONJG( AP( KK ) ) DO 190, K = KK + 1, KK + N - J IX = IX + INCX TEMP = TEMP + DCONJG( AP( K ) )*X( IX ) 190 CONTINUE END IF X( JX ) = TEMP JX = JX + INCX KK = KK + ( N - J + 1 ) 200 CONTINUE END IF END IF END IF * RETURN * * End of ZTPMV . * END SUBROUTINE ZTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. COMPLEX*16 AP( * ), X( * ) * .. * * Purpose * ======= * * ZTPSV solves one of the systems of equations * * A*x = b, or A'*x = b, or conjg( A' )*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix, supplied in packed form. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' or 'n' A*x = b. * * TRANS = 'T' or 't' A'*x = b. * * TRANS = 'C' or 'c' conjg( A' )*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * AP - COMPLEX*16 array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) * respectively, and so on. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) * respectively, and so on. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced, but are assumed to be unity. * Unchanged on exit. * * X - COMPLEX*16 array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. COMPLEX*16 TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX LOGICAL NOCONJ, NOUNIT * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( INCX.EQ.0 )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZTPSV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOCONJ = LSAME( TRANS, 'T' ) NOUNIT = LSAME( DIAG , 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of AP are * accessed sequentially with one pass through AP. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := inv( A )*x. * IF( LSAME( UPLO, 'U' ) )THEN KK = ( N*( N + 1 ) )/2 IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 c IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/AP( KK ) TEMP = X( J ) K = KK - 1 DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*AP( K ) K = K - 1 10 CONTINUE c END IF KK = KK - J 20 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 c IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/AP( KK ) TEMP = X( JX ) IX = JX DO 30, K = KK - 1, KK - J + 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*AP( K ) 30 CONTINUE c END IF JX = JX - INCX KK = KK - J 40 CONTINUE END IF ELSE KK = 1 IF( INCX.EQ.1 )THEN DO 60, J = 1, N c IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/AP( KK ) TEMP = X( J ) K = KK + 1 DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*AP( K ) K = K + 1 50 CONTINUE c END IF KK = KK + ( N - J + 1 ) 60 CONTINUE ELSE JX = KX DO 80, J = 1, N c IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/AP( KK ) TEMP = X( JX ) IX = JX DO 70, K = KK + 1, KK + N - J IX = IX + INCX X( IX ) = X( IX ) - TEMP*AP( K ) 70 CONTINUE c END IF JX = JX + INCX KK = KK + ( N - J + 1 ) 80 CONTINUE END IF END IF ELSE * * Form x := inv( A' )*x or x := inv( conjg( A' ) )*x. * IF( LSAME( UPLO, 'U' ) )THEN KK = 1 IF( INCX.EQ.1 )THEN DO 110, J = 1, N TEMP = X( J ) K = KK IF( NOCONJ )THEN DO 90, I = 1, J - 1 TEMP = TEMP - AP( K )*X( I ) K = K + 1 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK + J - 1 ) ELSE DO 100, I = 1, J - 1 TEMP = TEMP - DCONJG( AP( K ) )*X( I ) K = K + 1 100 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( AP( KK + J - 1 ) ) END IF X( J ) = TEMP KK = KK + J 110 CONTINUE ELSE JX = KX DO 140, J = 1, N TEMP = X( JX ) IX = KX IF( NOCONJ )THEN DO 120, K = KK, KK + J - 2 TEMP = TEMP - AP( K )*X( IX ) IX = IX + INCX 120 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK + J - 1 ) ELSE DO 130, K = KK, KK + J - 2 TEMP = TEMP - DCONJG( AP( K ) )*X( IX ) IX = IX + INCX 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( AP( KK + J - 1 ) ) END IF X( JX ) = TEMP JX = JX + INCX KK = KK + J 140 CONTINUE END IF ELSE KK = ( N*( N + 1 ) )/2 IF( INCX.EQ.1 )THEN DO 170, J = N, 1, -1 TEMP = X( J ) K = KK IF( NOCONJ )THEN DO 150, I = N, J + 1, -1 TEMP = TEMP - AP( K )*X( I ) K = K - 1 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK - N + J ) ELSE DO 160, I = N, J + 1, -1 TEMP = TEMP - DCONJG( AP( K ) )*X( I ) K = K - 1 160 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( AP( KK - N + J ) ) END IF X( J ) = TEMP KK = KK - ( N - J + 1 ) 170 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 200, J = N, 1, -1 TEMP = X( JX ) IX = KX IF( NOCONJ )THEN DO 180, K = KK, KK - ( N - ( J + 1 ) ), -1 TEMP = TEMP - AP( K )*X( IX ) IX = IX - INCX 180 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK - N + J ) ELSE DO 190, K = KK, KK - ( N - ( J + 1 ) ), -1 TEMP = TEMP - DCONJG( AP( K ) )*X( IX ) IX = IX - INCX 190 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/DCONJG( AP( KK - N + J ) ) END IF X( JX ) = TEMP JX = JX - INCX KK = KK - ( N - J + 1 ) 200 CONTINUE END IF END IF END IF * RETURN * * End of ZTPSV . * END SUBROUTINE ZGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC COMPLEX*16 ALPHA, BETA * .. Array Arguments .. COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * ZGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X' or op( X ) = conjg( X' ), * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = conjg( A' ). * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = conjg( B' ). * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - COMPLEX*16 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - COMPLEX*16 array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - COMPLEX*16 array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - 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 matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX * .. Local Scalars .. LOGICAL CONJA, CONJB, NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB COMPLEX*16 TEMP * .. Parameters .. COMPLEX*16 ONE PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * conjugated or transposed, set CONJA and CONJB as true if A and * B respectively are to be transposed but not conjugated and set * NROWA, NCOLA and NROWB as the number of rows and columns of A * and the number of rows of B respectively. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) CONJA = LSAME( TRANSA, 'C' ) CONJB = LSAME( TRANSB, 'C' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.CONJA ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.CONJB ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'ZGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And 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( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K c IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE c END IF 80 CONTINUE 90 CONTINUE ELSE IF( CONJA )THEN * * Form C := alpha*conjg( A' )*B + beta*C. * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 150, J = 1, N DO 140, I = 1, M TEMP = ZERO DO 130, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 130 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 140 CONTINUE 150 CONTINUE END IF ELSE IF( NOTA )THEN IF( CONJB )THEN * * Form C := alpha*A*conjg( B' ) + beta*C. * DO 200, J = 1, N IF( BETA.EQ.ZERO )THEN DO 160, I = 1, M C( I, J ) = ZERO 160 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 170, I = 1, M C( I, J ) = BETA*C( I, J ) 170 CONTINUE END IF DO 190, L = 1, K c IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*DCONJG( B( J, L ) ) DO 180, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 180 CONTINUE c END IF 190 CONTINUE 200 CONTINUE ELSE * * Form C := alpha*A*B' + beta*C * DO 250, J = 1, N IF( BETA.EQ.ZERO )THEN DO 210, I = 1, M C( I, J ) = ZERO 210 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 220, I = 1, M C( I, J ) = BETA*C( I, J ) 220 CONTINUE END IF DO 240, L = 1, K c IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 230, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 230 CONTINUE c END IF 240 CONTINUE 250 CONTINUE END IF ELSE IF( CONJA )THEN IF( CONJB )THEN * * Form C := alpha*conjg( A' )*conjg( B' ) + beta*C. * DO 280, J = 1, N DO 270, I = 1, M TEMP = ZERO DO 260, L = 1, K TEMP = TEMP + $ DCONJG( A( L, I ) )*DCONJG( B( J, L ) ) 260 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 270 CONTINUE 280 CONTINUE ELSE * * Form C := alpha*conjg( A' )*B' + beta*C * DO 310, J = 1, N DO 300, I = 1, M TEMP = ZERO DO 290, L = 1, K TEMP = TEMP + DCONJG( A( L, I ) )*B( J, L ) 290 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 300 CONTINUE 310 CONTINUE END IF ELSE IF( CONJB )THEN * * Form C := alpha*A'*conjg( B' ) + beta*C * DO 340, J = 1, N DO 330, I = 1, M TEMP = ZERO DO 320, L = 1, K TEMP = TEMP + A( L, I )*DCONJG( B( J, L ) ) 320 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 330 CONTINUE 340 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 370, J = 1, N DO 360, I = 1, M TEMP = ZERO DO 350, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 350 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 360 CONTINUE 370 CONTINUE END IF END IF * RETURN * * End of ZGEMM . * END