mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-03-25 18:50:40 +08:00
Remove unused fortran files
This commit is contained in:
Gael Guennebaud
parent
56ca44ad1a
commit
57ec399ec9
@ -1,310 +0,0 @@
|
||||
SUBROUTINE CHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX ALPHA,BETA
|
||||
INTEGER INCX,INCY,K,LDA,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX A(LDA,*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* CHBMV 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.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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 .
|
||||
* On entry, ALPHA specifies the scalar alpha.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* A - COMPLEX 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 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 .
|
||||
* On entry, BETA specifies the scalar beta.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* Y - COMPLEX 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* 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 ONE
|
||||
PARAMETER (ONE= (1.0E+0,0.0E+0))
|
||||
COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0E+0,0.0E+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
COMPLEX 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 CONJG,MAX,MIN,REAL
|
||||
* ..
|
||||
*
|
||||
* 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('CHBMV ',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 + CONJG(A(L+I,J))*X(I)
|
||||
50 CONTINUE
|
||||
Y(J) = Y(J) + TEMP1*REAL(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 + CONJG(A(L+I,J))*X(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
Y(JY) = Y(JY) + TEMP1*REAL(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*REAL(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 + CONJG(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*REAL(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 + CONJG(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 CHBMV .
|
||||
*
|
||||
END
|
||||
@ -1,272 +0,0 @@
|
||||
SUBROUTINE CHPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
COMPLEX ALPHA,BETA
|
||||
INTEGER INCX,INCY,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX AP(*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* CHPMV 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.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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 .
|
||||
* On entry, ALPHA specifies the scalar alpha.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* AP - COMPLEX 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 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 .
|
||||
* 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 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* 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 ONE
|
||||
PARAMETER (ONE= (1.0E+0,0.0E+0))
|
||||
COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0E+0,0.0E+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
COMPLEX 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 CONJG,REAL
|
||||
* ..
|
||||
*
|
||||
* 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('CHPMV ',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 + CONJG(AP(K))*X(I)
|
||||
K = K + 1
|
||||
50 CONTINUE
|
||||
Y(J) = Y(J) + TEMP1*REAL(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 + CONJG(AP(K))*X(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
Y(JY) = Y(JY) + TEMP1*REAL(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*REAL(AP(KK))
|
||||
K = KK + 1
|
||||
DO 90 I = J + 1,N
|
||||
Y(I) = Y(I) + TEMP1*AP(K)
|
||||
TEMP2 = TEMP2 + CONJG(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*REAL(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 + CONJG(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 CHPMV .
|
||||
*
|
||||
END
|
||||
@ -1,366 +0,0 @@
|
||||
SUBROUTINE CTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,K,LDA,N
|
||||
CHARACTER DIAG,TRANS,UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
COMPLEX A(LDA,*),X(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* CTBMV 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.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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 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 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* 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 ZERO
|
||||
PARAMETER (ZERO= (0.0E+0,0.0E+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
COMPLEX 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 CONJG,MAX,MIN
|
||||
* ..
|
||||
*
|
||||
* 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('CTBMV ',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
|
||||
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)
|
||||
END IF
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 40 J = 1,N
|
||||
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)
|
||||
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
|
||||
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)
|
||||
END IF
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 80 J = N,1,-1
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
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)
|
||||
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*CONJG(A(KPLUS1,J))
|
||||
DO 100 I = J - 1,MAX(1,J-K),-1
|
||||
TEMP = TEMP + CONJG(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*CONJG(A(KPLUS1,J))
|
||||
DO 130 I = J - 1,MAX(1,J-K),-1
|
||||
TEMP = TEMP + CONJG(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*CONJG(A(1,J))
|
||||
DO 160 I = J + 1,MIN(N,J+K)
|
||||
TEMP = TEMP + CONJG(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*CONJG(A(1,J))
|
||||
DO 190 I = J + 1,MIN(N,J+K)
|
||||
TEMP = TEMP + CONJG(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 CTBMV .
|
||||
*
|
||||
END
|
||||
@ -1,147 +0,0 @@
|
||||
SUBROUTINE DROTM(N,DX,INCX,DY,INCY,DPARAM)
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION DPARAM(5),DX(*),DY(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
|
||||
*
|
||||
* (DX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF DX ARE IN
|
||||
* (DY**T)
|
||||
*
|
||||
* DX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE
|
||||
* LX = (-INCX)*N, AND SIMILARLY FOR SY USING LY AND INCY.
|
||||
* WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
|
||||
*
|
||||
* DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0
|
||||
*
|
||||
* (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0)
|
||||
* H=( ) ( ) ( ) ( )
|
||||
* (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0).
|
||||
* SEE DROTMG FOR A DESCRIPTION OF DATA STORAGE IN DPARAM.
|
||||
*
|
||||
* Arguments
|
||||
* =========
|
||||
*
|
||||
* N (input) INTEGER
|
||||
* number of elements in input vector(s)
|
||||
*
|
||||
* DX (input/output) DOUBLE PRECISION array, dimension N
|
||||
* double precision vector with N elements
|
||||
*
|
||||
* INCX (input) INTEGER
|
||||
* storage spacing between elements of DX
|
||||
*
|
||||
* DY (input/output) DOUBLE PRECISION array, dimension N
|
||||
* double precision vector with N elements
|
||||
*
|
||||
* INCY (input) INTEGER
|
||||
* storage spacing between elements of DY
|
||||
*
|
||||
* DPARAM (input/output) DOUBLE PRECISION array, dimension 5
|
||||
* DPARAM(1)=DFLAG
|
||||
* DPARAM(2)=DH11
|
||||
* DPARAM(3)=DH21
|
||||
* DPARAM(4)=DH12
|
||||
* DPARAM(5)=DH22
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,TWO,W,Z,ZERO
|
||||
INTEGER I,KX,KY,NSTEPS
|
||||
* ..
|
||||
* .. Data statements ..
|
||||
DATA ZERO,TWO/0.D0,2.D0/
|
||||
* ..
|
||||
*
|
||||
DFLAG = DPARAM(1)
|
||||
IF (N.LE.0 .OR. (DFLAG+TWO.EQ.ZERO)) GO TO 140
|
||||
IF (.NOT. (INCX.EQ.INCY.AND.INCX.GT.0)) GO TO 70
|
||||
*
|
||||
NSTEPS = N*INCX
|
||||
IF (DFLAG) 50,10,30
|
||||
10 CONTINUE
|
||||
DH12 = DPARAM(4)
|
||||
DH21 = DPARAM(3)
|
||||
DO 20 I = 1,NSTEPS,INCX
|
||||
W = DX(I)
|
||||
Z = DY(I)
|
||||
DX(I) = W + Z*DH12
|
||||
DY(I) = W*DH21 + Z
|
||||
20 CONTINUE
|
||||
GO TO 140
|
||||
30 CONTINUE
|
||||
DH11 = DPARAM(2)
|
||||
DH22 = DPARAM(5)
|
||||
DO 40 I = 1,NSTEPS,INCX
|
||||
W = DX(I)
|
||||
Z = DY(I)
|
||||
DX(I) = W*DH11 + Z
|
||||
DY(I) = -W + DH22*Z
|
||||
40 CONTINUE
|
||||
GO TO 140
|
||||
50 CONTINUE
|
||||
DH11 = DPARAM(2)
|
||||
DH12 = DPARAM(4)
|
||||
DH21 = DPARAM(3)
|
||||
DH22 = DPARAM(5)
|
||||
DO 60 I = 1,NSTEPS,INCX
|
||||
W = DX(I)
|
||||
Z = DY(I)
|
||||
DX(I) = W*DH11 + Z*DH12
|
||||
DY(I) = W*DH21 + Z*DH22
|
||||
60 CONTINUE
|
||||
GO TO 140
|
||||
70 CONTINUE
|
||||
KX = 1
|
||||
KY = 1
|
||||
IF (INCX.LT.0) KX = 1 + (1-N)*INCX
|
||||
IF (INCY.LT.0) KY = 1 + (1-N)*INCY
|
||||
*
|
||||
IF (DFLAG) 120,80,100
|
||||
80 CONTINUE
|
||||
DH12 = DPARAM(4)
|
||||
DH21 = DPARAM(3)
|
||||
DO 90 I = 1,N
|
||||
W = DX(KX)
|
||||
Z = DY(KY)
|
||||
DX(KX) = W + Z*DH12
|
||||
DY(KY) = W*DH21 + Z
|
||||
KX = KX + INCX
|
||||
KY = KY + INCY
|
||||
90 CONTINUE
|
||||
GO TO 140
|
||||
100 CONTINUE
|
||||
DH11 = DPARAM(2)
|
||||
DH22 = DPARAM(5)
|
||||
DO 110 I = 1,N
|
||||
W = DX(KX)
|
||||
Z = DY(KY)
|
||||
DX(KX) = W*DH11 + Z
|
||||
DY(KY) = -W + DH22*Z
|
||||
KX = KX + INCX
|
||||
KY = KY + INCY
|
||||
110 CONTINUE
|
||||
GO TO 140
|
||||
120 CONTINUE
|
||||
DH11 = DPARAM(2)
|
||||
DH12 = DPARAM(4)
|
||||
DH21 = DPARAM(3)
|
||||
DH22 = DPARAM(5)
|
||||
DO 130 I = 1,N
|
||||
W = DX(KX)
|
||||
Z = DY(KY)
|
||||
DX(KX) = W*DH11 + Z*DH12
|
||||
DY(KY) = W*DH21 + Z*DH22
|
||||
KX = KX + INCX
|
||||
KY = KY + INCY
|
||||
130 CONTINUE
|
||||
140 CONTINUE
|
||||
RETURN
|
||||
END
|
||||
@ -1,206 +0,0 @@
|
||||
SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
|
||||
* .. Scalar Arguments ..
|
||||
DOUBLE PRECISION DD1,DD2,DX1,DY1
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION DPARAM(5)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
|
||||
* THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)*
|
||||
* DY2)**T.
|
||||
* WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
|
||||
*
|
||||
* DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0
|
||||
*
|
||||
* (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0)
|
||||
* H=( ) ( ) ( ) ( )
|
||||
* (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0).
|
||||
* LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22
|
||||
* RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE
|
||||
* VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.)
|
||||
*
|
||||
* THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
|
||||
* INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
|
||||
* OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
|
||||
*
|
||||
*
|
||||
* Arguments
|
||||
* =========
|
||||
*
|
||||
* DD1 (input/output) DOUBLE PRECISION
|
||||
*
|
||||
* DD2 (input/output) DOUBLE PRECISION
|
||||
*
|
||||
* DX1 (input/output) DOUBLE PRECISION
|
||||
*
|
||||
* DY1 (input) DOUBLE PRECISION
|
||||
*
|
||||
* DPARAM (input/output) DOUBLE PRECISION array, dimension 5
|
||||
* DPARAM(1)=DFLAG
|
||||
* DPARAM(2)=DH11
|
||||
* DPARAM(3)=DH21
|
||||
* DPARAM(4)=DH12
|
||||
* DPARAM(5)=DH22
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP,
|
||||
+ DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO
|
||||
INTEGER IGO
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC DABS
|
||||
* ..
|
||||
* .. Data statements ..
|
||||
*
|
||||
DATA ZERO,ONE,TWO/0.D0,1.D0,2.D0/
|
||||
DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/
|
||||
* ..
|
||||
|
||||
IF (.NOT.DD1.LT.ZERO) GO TO 10
|
||||
* GO ZERO-H-D-AND-DX1..
|
||||
GO TO 60
|
||||
10 CONTINUE
|
||||
* CASE-DD1-NONNEGATIVE
|
||||
DP2 = DD2*DY1
|
||||
IF (.NOT.DP2.EQ.ZERO) GO TO 20
|
||||
DFLAG = -TWO
|
||||
GO TO 260
|
||||
* REGULAR-CASE..
|
||||
20 CONTINUE
|
||||
DP1 = DD1*DX1
|
||||
DQ2 = DP2*DY1
|
||||
DQ1 = DP1*DX1
|
||||
*
|
||||
IF (.NOT.DABS(DQ1).GT.DABS(DQ2)) GO TO 40
|
||||
DH21 = -DY1/DX1
|
||||
DH12 = DP2/DP1
|
||||
*
|
||||
DU = ONE - DH12*DH21
|
||||
*
|
||||
IF (.NOT.DU.LE.ZERO) GO TO 30
|
||||
* GO ZERO-H-D-AND-DX1..
|
||||
GO TO 60
|
||||
30 CONTINUE
|
||||
DFLAG = ZERO
|
||||
DD1 = DD1/DU
|
||||
DD2 = DD2/DU
|
||||
DX1 = DX1*DU
|
||||
* GO SCALE-CHECK..
|
||||
GO TO 100
|
||||
40 CONTINUE
|
||||
IF (.NOT.DQ2.LT.ZERO) GO TO 50
|
||||
* GO ZERO-H-D-AND-DX1..
|
||||
GO TO 60
|
||||
50 CONTINUE
|
||||
DFLAG = ONE
|
||||
DH11 = DP1/DP2
|
||||
DH22 = DX1/DY1
|
||||
DU = ONE + DH11*DH22
|
||||
DTEMP = DD2/DU
|
||||
DD2 = DD1/DU
|
||||
DD1 = DTEMP
|
||||
DX1 = DY1*DU
|
||||
* GO SCALE-CHECK
|
||||
GO TO 100
|
||||
* PROCEDURE..ZERO-H-D-AND-DX1..
|
||||
60 CONTINUE
|
||||
DFLAG = -ONE
|
||||
DH11 = ZERO
|
||||
DH12 = ZERO
|
||||
DH21 = ZERO
|
||||
DH22 = ZERO
|
||||
*
|
||||
DD1 = ZERO
|
||||
DD2 = ZERO
|
||||
DX1 = ZERO
|
||||
* RETURN..
|
||||
GO TO 220
|
||||
* PROCEDURE..FIX-H..
|
||||
70 CONTINUE
|
||||
IF (.NOT.DFLAG.GE.ZERO) GO TO 90
|
||||
*
|
||||
IF (.NOT.DFLAG.EQ.ZERO) GO TO 80
|
||||
DH11 = ONE
|
||||
DH22 = ONE
|
||||
DFLAG = -ONE
|
||||
GO TO 90
|
||||
80 CONTINUE
|
||||
DH21 = -ONE
|
||||
DH12 = ONE
|
||||
DFLAG = -ONE
|
||||
90 CONTINUE
|
||||
GO TO IGO(120,150,180,210)
|
||||
* PROCEDURE..SCALE-CHECK
|
||||
100 CONTINUE
|
||||
110 CONTINUE
|
||||
IF (.NOT.DD1.LE.RGAMSQ) GO TO 130
|
||||
IF (DD1.EQ.ZERO) GO TO 160
|
||||
ASSIGN 120 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
120 CONTINUE
|
||||
DD1 = DD1*GAM**2
|
||||
DX1 = DX1/GAM
|
||||
DH11 = DH11/GAM
|
||||
DH12 = DH12/GAM
|
||||
GO TO 110
|
||||
130 CONTINUE
|
||||
140 CONTINUE
|
||||
IF (.NOT.DD1.GE.GAMSQ) GO TO 160
|
||||
ASSIGN 150 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
150 CONTINUE
|
||||
DD1 = DD1/GAM**2
|
||||
DX1 = DX1*GAM
|
||||
DH11 = DH11*GAM
|
||||
DH12 = DH12*GAM
|
||||
GO TO 140
|
||||
160 CONTINUE
|
||||
170 CONTINUE
|
||||
IF (.NOT.DABS(DD2).LE.RGAMSQ) GO TO 190
|
||||
IF (DD2.EQ.ZERO) GO TO 220
|
||||
ASSIGN 180 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
180 CONTINUE
|
||||
DD2 = DD2*GAM**2
|
||||
DH21 = DH21/GAM
|
||||
DH22 = DH22/GAM
|
||||
GO TO 170
|
||||
190 CONTINUE
|
||||
200 CONTINUE
|
||||
IF (.NOT.DABS(DD2).GE.GAMSQ) GO TO 220
|
||||
ASSIGN 210 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
210 CONTINUE
|
||||
DD2 = DD2/GAM**2
|
||||
DH21 = DH21*GAM
|
||||
DH22 = DH22*GAM
|
||||
GO TO 200
|
||||
220 CONTINUE
|
||||
IF (DFLAG) 250,230,240
|
||||
230 CONTINUE
|
||||
DPARAM(3) = DH21
|
||||
DPARAM(4) = DH12
|
||||
GO TO 260
|
||||
240 CONTINUE
|
||||
DPARAM(2) = DH11
|
||||
DPARAM(5) = DH22
|
||||
GO TO 260
|
||||
250 CONTINUE
|
||||
DPARAM(2) = DH11
|
||||
DPARAM(3) = DH21
|
||||
DPARAM(4) = DH12
|
||||
DPARAM(5) = DH22
|
||||
260 CONTINUE
|
||||
DPARAM(1) = DFLAG
|
||||
RETURN
|
||||
END
|
||||
@ -1,304 +0,0 @@
|
||||
SUBROUTINE DSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
DOUBLE PRECISION ALPHA,BETA
|
||||
INTEGER INCX,INCY,K,LDA,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION A(LDA,*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* DSBMV 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 symmetric band matrix, with k super-diagonals.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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 - DOUBLE PRECISION.
|
||||
* On entry, ALPHA specifies the scalar alpha.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* A - DOUBLE PRECISION 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 symmetric 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 symmetric 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 symmetric 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 symmetric 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
|
||||
*
|
||||
* 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION.
|
||||
* On entry, BETA specifies the scalar beta.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* Y - DOUBLE PRECISION 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 ..
|
||||
DOUBLE PRECISION ONE,ZERO
|
||||
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
DOUBLE PRECISION 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 MAX,MIN
|
||||
* ..
|
||||
*
|
||||
* 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('DSBMV ',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 + A(L+I,J)*X(I)
|
||||
50 CONTINUE
|
||||
Y(J) = Y(J) + TEMP1*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 + A(L+I,J)*X(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
Y(JY) = Y(JY) + TEMP1*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*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 + 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*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 + 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 DSBMV .
|
||||
*
|
||||
END
|
||||
@ -1,265 +0,0 @@
|
||||
SUBROUTINE DSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
DOUBLE PRECISION ALPHA,BETA
|
||||
INTEGER INCX,INCY,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION AP(*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* DSPMV 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 symmetric matrix, supplied in packed form.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* AP - DOUBLE PRECISION 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 symmetric 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 symmetric 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.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* X - DOUBLE PRECISION array of dimension at least
|
||||
* ( 1 + ( n - 1 )*abs( INCX ) ).
|
||||
* Before entry, the incremented array X must contain the n
|
||||
* element vector x.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* INCX - INTEGER.
|
||||
* On entry, INCX specifies the increment for the elements of
|
||||
* X. INCX must not be zero.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* BETA - DOUBLE PRECISION.
|
||||
* On entry, BETA specifies the scalar beta. When BETA is
|
||||
* supplied as zero then Y need not be set on input.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* Y - DOUBLE PRECISION array of dimension at least
|
||||
* ( 1 + ( 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* Level 2 Blas routine.
|
||||
*
|
||||
* -- Written on 22-October-1986.
|
||||
* Jack Dongarra, Argonne National Lab.
|
||||
* Jeremy Du Croz, Nag Central Office.
|
||||
* Sven Hammarling, Nag Central Office.
|
||||
* Richard Hanson, Sandia National Labs.
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
DOUBLE PRECISION ONE,ZERO
|
||||
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
DOUBLE PRECISION TEMP1,TEMP2
|
||||
INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
*
|
||||
* 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('DSPMV ',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 + AP(K)*X(I)
|
||||
K = K + 1
|
||||
50 CONTINUE
|
||||
Y(J) = Y(J) + TEMP1*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 + AP(K)*X(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
Y(JY) = Y(JY) + TEMP1*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*AP(KK)
|
||||
K = KK + 1
|
||||
DO 90 I = J + 1,N
|
||||
Y(I) = Y(I) + TEMP1*AP(K)
|
||||
TEMP2 = TEMP2 + 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*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 + 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 DSPMV .
|
||||
*
|
||||
END
|
||||
@ -1,335 +0,0 @@
|
||||
SUBROUTINE DTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,K,LDA,N
|
||||
CHARACTER DIAG,TRANS,UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION A(LDA,*),X(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* DTBMV performs one of the matrix-vector operations
|
||||
*
|
||||
* x := A*x, or x := A'*x,
|
||||
*
|
||||
* where x is an n element vector and A is an n by n unit, or non-unit,
|
||||
* upper or lower triangular band matrix, with ( k + 1 ) diagonals.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* UPLO - CHARACTER*1.
|
||||
* On entry, UPLO specifies whether the matrix is an upper or
|
||||
* lower triangular matrix as follows:
|
||||
*
|
||||
* UPLO = 'U' or 'u' A is an upper triangular matrix.
|
||||
*
|
||||
* UPLO = 'L' or 'l' A is a lower triangular matrix.
|
||||
*
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* TRANS - CHARACTER*1.
|
||||
* On entry, TRANS specifies the operation to be performed as
|
||||
* follows:
|
||||
*
|
||||
* TRANS = 'N' or 'n' x := A*x.
|
||||
*
|
||||
* TRANS = 'T' or 't' x := A'*x.
|
||||
*
|
||||
* TRANS = 'C' or 'c' x := A'*x.
|
||||
*
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* DIAG - CHARACTER*1.
|
||||
* On entry, DIAG specifies whether or not A is unit
|
||||
* triangular as follows:
|
||||
*
|
||||
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
|
||||
*
|
||||
* DIAG = 'N' or 'n' A is not assumed to be unit
|
||||
* triangular.
|
||||
*
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* N - INTEGER.
|
||||
* On entry, N specifies the order of the matrix A.
|
||||
* N must be at least zero.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* 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 - DOUBLE PRECISION 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 - DOUBLE PRECISION array of dimension at least
|
||||
* ( 1 + ( n - 1 )*abs( INCX ) ).
|
||||
* Before entry, the incremented array X must contain the n
|
||||
* element vector x. On exit, X is overwritten with the
|
||||
* tranformed vector x.
|
||||
*
|
||||
* INCX - INTEGER.
|
||||
* On entry, INCX specifies the increment for the elements of
|
||||
* X. INCX must not be zero.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* Level 2 Blas routine.
|
||||
*
|
||||
* -- Written on 22-October-1986.
|
||||
* Jack Dongarra, Argonne National Lab.
|
||||
* Jeremy Du Croz, Nag Central Office.
|
||||
* Sven Hammarling, Nag Central Office.
|
||||
* Richard Hanson, Sandia National Labs.
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
DOUBLE PRECISION ZERO
|
||||
PARAMETER (ZERO=0.0D+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
DOUBLE PRECISION TEMP
|
||||
INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L
|
||||
LOGICAL NOUNIT
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX,MIN
|
||||
* ..
|
||||
*
|
||||
* 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('DTBMV ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF (N.EQ.0) RETURN
|
||||
*
|
||||
NOUNIT = LSAME(DIAG,'N')
|
||||
*
|
||||
* Set up the start point in X if the increment is not unity. This
|
||||
* will be ( N - 1 )*INCX too small for descending loops.
|
||||
*
|
||||
IF (INCX.LE.0) THEN
|
||||
KX = 1 - (N-1)*INCX
|
||||
ELSE IF (INCX.NE.1) THEN
|
||||
KX = 1
|
||||
END IF
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
*
|
||||
* Form x := A*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
KPLUS1 = K + 1
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 20 J = 1,N
|
||||
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)
|
||||
END IF
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 40 J = 1,N
|
||||
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)
|
||||
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
|
||||
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)
|
||||
END IF
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 80 J = N,1,-1
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
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)
|
||||
END IF
|
||||
JX = JX - INCX
|
||||
IF ((N-J).GE.K) KX = KX - INCX
|
||||
80 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form x := A'*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
KPLUS1 = K + 1
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 100 J = N,1,-1
|
||||
TEMP = X(J)
|
||||
L = KPLUS1 - J
|
||||
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
|
||||
X(J) = TEMP
|
||||
100 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 120 J = N,1,-1
|
||||
TEMP = X(JX)
|
||||
KX = KX - INCX
|
||||
IX = KX
|
||||
L = KPLUS1 - J
|
||||
IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
|
||||
DO 110 I = J - 1,MAX(1,J-K),-1
|
||||
TEMP = TEMP + A(L+I,J)*X(IX)
|
||||
IX = IX - INCX
|
||||
110 CONTINUE
|
||||
X(JX) = TEMP
|
||||
JX = JX - INCX
|
||||
120 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 140 J = 1,N
|
||||
TEMP = X(J)
|
||||
L = 1 - J
|
||||
IF (NOUNIT) TEMP = TEMP*A(1,J)
|
||||
DO 130 I = J + 1,MIN(N,J+K)
|
||||
TEMP = TEMP + A(L+I,J)*X(I)
|
||||
130 CONTINUE
|
||||
X(J) = TEMP
|
||||
140 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 160 J = 1,N
|
||||
TEMP = X(JX)
|
||||
KX = KX + INCX
|
||||
IX = KX
|
||||
L = 1 - J
|
||||
IF (NOUNIT) TEMP = TEMP*A(1,J)
|
||||
DO 150 I = J + 1,MIN(N,J+K)
|
||||
TEMP = TEMP + A(L+I,J)*X(IX)
|
||||
IX = IX + INCX
|
||||
150 CONTINUE
|
||||
X(JX) = TEMP
|
||||
JX = JX + INCX
|
||||
160 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of DTBMV .
|
||||
*
|
||||
END
|
||||
@ -1,85 +0,0 @@
|
||||
LOGICAL FUNCTION LSAME(CA,CB)
|
||||
*
|
||||
* -- LAPACK auxiliary routine (version 3.1) --
|
||||
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
|
||||
* November 2006
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
CHARACTER CA,CB
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* LSAME returns .TRUE. if CA is the same letter as CB regardless of
|
||||
* case.
|
||||
*
|
||||
* Arguments
|
||||
* =========
|
||||
*
|
||||
* CA (input) CHARACTER*1
|
||||
*
|
||||
* CB (input) CHARACTER*1
|
||||
* CA and CB specify the single characters to be compared.
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ICHAR
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
INTEGER INTA,INTB,ZCODE
|
||||
* ..
|
||||
*
|
||||
* Test if the characters are equal
|
||||
*
|
||||
LSAME = CA .EQ. CB
|
||||
IF (LSAME) RETURN
|
||||
*
|
||||
* Now test for equivalence if both characters are alphabetic.
|
||||
*
|
||||
ZCODE = ICHAR('Z')
|
||||
*
|
||||
* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
|
||||
* machines, on which ICHAR returns a value with bit 8 set.
|
||||
* ICHAR('A') on Prime machines returns 193 which is the same as
|
||||
* ICHAR('A') on an EBCDIC machine.
|
||||
*
|
||||
INTA = ICHAR(CA)
|
||||
INTB = ICHAR(CB)
|
||||
*
|
||||
IF (ZCODE.EQ.90 .OR. ZCODE.EQ.122) THEN
|
||||
*
|
||||
* ASCII is assumed - ZCODE is the ASCII code of either lower or
|
||||
* upper case 'Z'.
|
||||
*
|
||||
IF (INTA.GE.97 .AND. INTA.LE.122) INTA = INTA - 32
|
||||
IF (INTB.GE.97 .AND. INTB.LE.122) INTB = INTB - 32
|
||||
*
|
||||
ELSE IF (ZCODE.EQ.233 .OR. ZCODE.EQ.169) THEN
|
||||
*
|
||||
* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
|
||||
* upper case 'Z'.
|
||||
*
|
||||
IF (INTA.GE.129 .AND. INTA.LE.137 .OR.
|
||||
+ INTA.GE.145 .AND. INTA.LE.153 .OR.
|
||||
+ INTA.GE.162 .AND. INTA.LE.169) INTA = INTA + 64
|
||||
IF (INTB.GE.129 .AND. INTB.LE.137 .OR.
|
||||
+ INTB.GE.145 .AND. INTB.LE.153 .OR.
|
||||
+ INTB.GE.162 .AND. INTB.LE.169) INTB = INTB + 64
|
||||
*
|
||||
ELSE IF (ZCODE.EQ.218 .OR. ZCODE.EQ.250) THEN
|
||||
*
|
||||
* ASCII is assumed, on Prime machines - ZCODE is the ASCII code
|
||||
* plus 128 of either lower or upper case 'Z'.
|
||||
*
|
||||
IF (INTA.GE.225 .AND. INTA.LE.250) INTA = INTA - 32
|
||||
IF (INTB.GE.225 .AND. INTB.LE.250) INTB = INTB - 32
|
||||
END IF
|
||||
LSAME = INTA .EQ. INTB
|
||||
*
|
||||
* RETURN
|
||||
*
|
||||
* End of LSAME
|
||||
*
|
||||
END
|
||||
@ -1,148 +0,0 @@
|
||||
SUBROUTINE SROTM(N,SX,INCX,SY,INCY,SPARAM)
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,INCY,N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL SPARAM(5),SX(*),SY(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
|
||||
*
|
||||
* (SX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF SX ARE IN
|
||||
* (DX**T)
|
||||
*
|
||||
* SX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE
|
||||
* LX = (-INCX)*N, AND SIMILARLY FOR SY USING USING LY AND INCY.
|
||||
* WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
|
||||
*
|
||||
* SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0
|
||||
*
|
||||
* (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0)
|
||||
* H=( ) ( ) ( ) ( )
|
||||
* (SH21 SH22), (SH21 1.E0), (-1.E0 SH22), (0.E0 1.E0).
|
||||
* SEE SROTMG FOR A DESCRIPTION OF DATA STORAGE IN SPARAM.
|
||||
*
|
||||
*
|
||||
* Arguments
|
||||
* =========
|
||||
*
|
||||
* N (input) INTEGER
|
||||
* number of elements in input vector(s)
|
||||
*
|
||||
* SX (input/output) REAL array, dimension N
|
||||
* double precision vector with N elements
|
||||
*
|
||||
* INCX (input) INTEGER
|
||||
* storage spacing between elements of SX
|
||||
*
|
||||
* SY (input/output) REAL array, dimension N
|
||||
* double precision vector with N elements
|
||||
*
|
||||
* INCY (input) INTEGER
|
||||
* storage spacing between elements of SY
|
||||
*
|
||||
* SPARAM (input/output) REAL array, dimension 5
|
||||
* SPARAM(1)=SFLAG
|
||||
* SPARAM(2)=SH11
|
||||
* SPARAM(3)=SH21
|
||||
* SPARAM(4)=SH12
|
||||
* SPARAM(5)=SH22
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
REAL SFLAG,SH11,SH12,SH21,SH22,TWO,W,Z,ZERO
|
||||
INTEGER I,KX,KY,NSTEPS
|
||||
* ..
|
||||
* .. Data statements ..
|
||||
DATA ZERO,TWO/0.E0,2.E0/
|
||||
* ..
|
||||
*
|
||||
SFLAG = SPARAM(1)
|
||||
IF (N.LE.0 .OR. (SFLAG+TWO.EQ.ZERO)) GO TO 140
|
||||
IF (.NOT. (INCX.EQ.INCY.AND.INCX.GT.0)) GO TO 70
|
||||
*
|
||||
NSTEPS = N*INCX
|
||||
IF (SFLAG) 50,10,30
|
||||
10 CONTINUE
|
||||
SH12 = SPARAM(4)
|
||||
SH21 = SPARAM(3)
|
||||
DO 20 I = 1,NSTEPS,INCX
|
||||
W = SX(I)
|
||||
Z = SY(I)
|
||||
SX(I) = W + Z*SH12
|
||||
SY(I) = W*SH21 + Z
|
||||
20 CONTINUE
|
||||
GO TO 140
|
||||
30 CONTINUE
|
||||
SH11 = SPARAM(2)
|
||||
SH22 = SPARAM(5)
|
||||
DO 40 I = 1,NSTEPS,INCX
|
||||
W = SX(I)
|
||||
Z = SY(I)
|
||||
SX(I) = W*SH11 + Z
|
||||
SY(I) = -W + SH22*Z
|
||||
40 CONTINUE
|
||||
GO TO 140
|
||||
50 CONTINUE
|
||||
SH11 = SPARAM(2)
|
||||
SH12 = SPARAM(4)
|
||||
SH21 = SPARAM(3)
|
||||
SH22 = SPARAM(5)
|
||||
DO 60 I = 1,NSTEPS,INCX
|
||||
W = SX(I)
|
||||
Z = SY(I)
|
||||
SX(I) = W*SH11 + Z*SH12
|
||||
SY(I) = W*SH21 + Z*SH22
|
||||
60 CONTINUE
|
||||
GO TO 140
|
||||
70 CONTINUE
|
||||
KX = 1
|
||||
KY = 1
|
||||
IF (INCX.LT.0) KX = 1 + (1-N)*INCX
|
||||
IF (INCY.LT.0) KY = 1 + (1-N)*INCY
|
||||
*
|
||||
IF (SFLAG) 120,80,100
|
||||
80 CONTINUE
|
||||
SH12 = SPARAM(4)
|
||||
SH21 = SPARAM(3)
|
||||
DO 90 I = 1,N
|
||||
W = SX(KX)
|
||||
Z = SY(KY)
|
||||
SX(KX) = W + Z*SH12
|
||||
SY(KY) = W*SH21 + Z
|
||||
KX = KX + INCX
|
||||
KY = KY + INCY
|
||||
90 CONTINUE
|
||||
GO TO 140
|
||||
100 CONTINUE
|
||||
SH11 = SPARAM(2)
|
||||
SH22 = SPARAM(5)
|
||||
DO 110 I = 1,N
|
||||
W = SX(KX)
|
||||
Z = SY(KY)
|
||||
SX(KX) = W*SH11 + Z
|
||||
SY(KY) = -W + SH22*Z
|
||||
KX = KX + INCX
|
||||
KY = KY + INCY
|
||||
110 CONTINUE
|
||||
GO TO 140
|
||||
120 CONTINUE
|
||||
SH11 = SPARAM(2)
|
||||
SH12 = SPARAM(4)
|
||||
SH21 = SPARAM(3)
|
||||
SH22 = SPARAM(5)
|
||||
DO 130 I = 1,N
|
||||
W = SX(KX)
|
||||
Z = SY(KY)
|
||||
SX(KX) = W*SH11 + Z*SH12
|
||||
SY(KY) = W*SH21 + Z*SH22
|
||||
KX = KX + INCX
|
||||
KY = KY + INCY
|
||||
130 CONTINUE
|
||||
140 CONTINUE
|
||||
RETURN
|
||||
END
|
||||
@ -1,208 +0,0 @@
|
||||
SUBROUTINE SROTMG(SD1,SD2,SX1,SY1,SPARAM)
|
||||
* .. Scalar Arguments ..
|
||||
REAL SD1,SD2,SX1,SY1
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL SPARAM(5)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
|
||||
* THE SECOND COMPONENT OF THE 2-VECTOR (SQRT(SD1)*SX1,SQRT(SD2)*
|
||||
* SY2)**T.
|
||||
* WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
|
||||
*
|
||||
* SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0
|
||||
*
|
||||
* (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0)
|
||||
* H=( ) ( ) ( ) ( )
|
||||
* (SH21 SH22), (SH21 1.E0), (-1.E0 SH22), (0.E0 1.E0).
|
||||
* LOCATIONS 2-4 OF SPARAM CONTAIN SH11,SH21,SH12, AND SH22
|
||||
* RESPECTIVELY. (VALUES OF 1.E0, -1.E0, OR 0.E0 IMPLIED BY THE
|
||||
* VALUE OF SPARAM(1) ARE NOT STORED IN SPARAM.)
|
||||
*
|
||||
* THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
|
||||
* INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
|
||||
* OF SD1 AND SD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
|
||||
*
|
||||
*
|
||||
* Arguments
|
||||
* =========
|
||||
*
|
||||
*
|
||||
* SD1 (input/output) REAL
|
||||
*
|
||||
* SD2 (input/output) REAL
|
||||
*
|
||||
* SX1 (input/output) REAL
|
||||
*
|
||||
* SY1 (input) REAL
|
||||
*
|
||||
*
|
||||
* SPARAM (input/output) REAL array, dimension 5
|
||||
* SPARAM(1)=SFLAG
|
||||
* SPARAM(2)=SH11
|
||||
* SPARAM(3)=SH21
|
||||
* SPARAM(4)=SH12
|
||||
* SPARAM(5)=SH22
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Local Scalars ..
|
||||
REAL GAM,GAMSQ,ONE,RGAMSQ,SFLAG,SH11,SH12,SH21,SH22,SP1,SP2,SQ1,
|
||||
+ SQ2,STEMP,SU,TWO,ZERO
|
||||
INTEGER IGO
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC ABS
|
||||
* ..
|
||||
* .. Data statements ..
|
||||
*
|
||||
DATA ZERO,ONE,TWO/0.E0,1.E0,2.E0/
|
||||
DATA GAM,GAMSQ,RGAMSQ/4096.E0,1.67772E7,5.96046E-8/
|
||||
* ..
|
||||
|
||||
IF (.NOT.SD1.LT.ZERO) GO TO 10
|
||||
* GO ZERO-H-D-AND-SX1..
|
||||
GO TO 60
|
||||
10 CONTINUE
|
||||
* CASE-SD1-NONNEGATIVE
|
||||
SP2 = SD2*SY1
|
||||
IF (.NOT.SP2.EQ.ZERO) GO TO 20
|
||||
SFLAG = -TWO
|
||||
GO TO 260
|
||||
* REGULAR-CASE..
|
||||
20 CONTINUE
|
||||
SP1 = SD1*SX1
|
||||
SQ2 = SP2*SY1
|
||||
SQ1 = SP1*SX1
|
||||
*
|
||||
IF (.NOT.ABS(SQ1).GT.ABS(SQ2)) GO TO 40
|
||||
SH21 = -SY1/SX1
|
||||
SH12 = SP2/SP1
|
||||
*
|
||||
SU = ONE - SH12*SH21
|
||||
*
|
||||
IF (.NOT.SU.LE.ZERO) GO TO 30
|
||||
* GO ZERO-H-D-AND-SX1..
|
||||
GO TO 60
|
||||
30 CONTINUE
|
||||
SFLAG = ZERO
|
||||
SD1 = SD1/SU
|
||||
SD2 = SD2/SU
|
||||
SX1 = SX1*SU
|
||||
* GO SCALE-CHECK..
|
||||
GO TO 100
|
||||
40 CONTINUE
|
||||
IF (.NOT.SQ2.LT.ZERO) GO TO 50
|
||||
* GO ZERO-H-D-AND-SX1..
|
||||
GO TO 60
|
||||
50 CONTINUE
|
||||
SFLAG = ONE
|
||||
SH11 = SP1/SP2
|
||||
SH22 = SX1/SY1
|
||||
SU = ONE + SH11*SH22
|
||||
STEMP = SD2/SU
|
||||
SD2 = SD1/SU
|
||||
SD1 = STEMP
|
||||
SX1 = SY1*SU
|
||||
* GO SCALE-CHECK
|
||||
GO TO 100
|
||||
* PROCEDURE..ZERO-H-D-AND-SX1..
|
||||
60 CONTINUE
|
||||
SFLAG = -ONE
|
||||
SH11 = ZERO
|
||||
SH12 = ZERO
|
||||
SH21 = ZERO
|
||||
SH22 = ZERO
|
||||
*
|
||||
SD1 = ZERO
|
||||
SD2 = ZERO
|
||||
SX1 = ZERO
|
||||
* RETURN..
|
||||
GO TO 220
|
||||
* PROCEDURE..FIX-H..
|
||||
70 CONTINUE
|
||||
IF (.NOT.SFLAG.GE.ZERO) GO TO 90
|
||||
*
|
||||
IF (.NOT.SFLAG.EQ.ZERO) GO TO 80
|
||||
SH11 = ONE
|
||||
SH22 = ONE
|
||||
SFLAG = -ONE
|
||||
GO TO 90
|
||||
80 CONTINUE
|
||||
SH21 = -ONE
|
||||
SH12 = ONE
|
||||
SFLAG = -ONE
|
||||
90 CONTINUE
|
||||
GO TO IGO(120,150,180,210)
|
||||
* PROCEDURE..SCALE-CHECK
|
||||
100 CONTINUE
|
||||
110 CONTINUE
|
||||
IF (.NOT.SD1.LE.RGAMSQ) GO TO 130
|
||||
IF (SD1.EQ.ZERO) GO TO 160
|
||||
ASSIGN 120 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
120 CONTINUE
|
||||
SD1 = SD1*GAM**2
|
||||
SX1 = SX1/GAM
|
||||
SH11 = SH11/GAM
|
||||
SH12 = SH12/GAM
|
||||
GO TO 110
|
||||
130 CONTINUE
|
||||
140 CONTINUE
|
||||
IF (.NOT.SD1.GE.GAMSQ) GO TO 160
|
||||
ASSIGN 150 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
150 CONTINUE
|
||||
SD1 = SD1/GAM**2
|
||||
SX1 = SX1*GAM
|
||||
SH11 = SH11*GAM
|
||||
SH12 = SH12*GAM
|
||||
GO TO 140
|
||||
160 CONTINUE
|
||||
170 CONTINUE
|
||||
IF (.NOT.ABS(SD2).LE.RGAMSQ) GO TO 190
|
||||
IF (SD2.EQ.ZERO) GO TO 220
|
||||
ASSIGN 180 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
180 CONTINUE
|
||||
SD2 = SD2*GAM**2
|
||||
SH21 = SH21/GAM
|
||||
SH22 = SH22/GAM
|
||||
GO TO 170
|
||||
190 CONTINUE
|
||||
200 CONTINUE
|
||||
IF (.NOT.ABS(SD2).GE.GAMSQ) GO TO 220
|
||||
ASSIGN 210 TO IGO
|
||||
* FIX-H..
|
||||
GO TO 70
|
||||
210 CONTINUE
|
||||
SD2 = SD2/GAM**2
|
||||
SH21 = SH21*GAM
|
||||
SH22 = SH22*GAM
|
||||
GO TO 200
|
||||
220 CONTINUE
|
||||
IF (SFLAG) 250,230,240
|
||||
230 CONTINUE
|
||||
SPARAM(3) = SH21
|
||||
SPARAM(4) = SH12
|
||||
GO TO 260
|
||||
240 CONTINUE
|
||||
SPARAM(2) = SH11
|
||||
SPARAM(5) = SH22
|
||||
GO TO 260
|
||||
250 CONTINUE
|
||||
SPARAM(2) = SH11
|
||||
SPARAM(3) = SH21
|
||||
SPARAM(4) = SH12
|
||||
SPARAM(5) = SH22
|
||||
260 CONTINUE
|
||||
SPARAM(1) = SFLAG
|
||||
RETURN
|
||||
END
|
||||
@ -1,306 +0,0 @@
|
||||
SUBROUTINE SSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
REAL ALPHA,BETA
|
||||
INTEGER INCX,INCY,K,LDA,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL A(LDA,*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* SSBMV 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 symmetric band matrix, with k super-diagonals.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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 - REAL .
|
||||
* On entry, ALPHA specifies the scalar alpha.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* A - REAL 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 symmetric 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 symmetric 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 symmetric 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 symmetric 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
|
||||
*
|
||||
* 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 - REAL 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 - REAL .
|
||||
* On entry, BETA specifies the scalar beta.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* Y - REAL 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* 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 ..
|
||||
REAL ONE,ZERO
|
||||
PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
REAL 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 MAX,MIN
|
||||
* ..
|
||||
*
|
||||
* 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('SSBMV ',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 + A(L+I,J)*X(I)
|
||||
50 CONTINUE
|
||||
Y(J) = Y(J) + TEMP1*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 + A(L+I,J)*X(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
Y(JY) = Y(JY) + TEMP1*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*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 + 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*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 + 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 SSBMV .
|
||||
*
|
||||
END
|
||||
@ -1,265 +0,0 @@
|
||||
SUBROUTINE SSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
REAL ALPHA,BETA
|
||||
INTEGER INCX,INCY,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL AP(*),X(*),Y(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* SSPMV 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 symmetric matrix, supplied in packed form.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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 - REAL .
|
||||
* On entry, ALPHA specifies the scalar alpha.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* AP - REAL 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 symmetric 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 symmetric 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.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* X - REAL 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 - REAL .
|
||||
* 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 - REAL 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* 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 ..
|
||||
REAL ONE,ZERO
|
||||
PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
REAL TEMP1,TEMP2
|
||||
INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
*
|
||||
* 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('SSPMV ',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 + AP(K)*X(I)
|
||||
K = K + 1
|
||||
50 CONTINUE
|
||||
Y(J) = Y(J) + TEMP1*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 + AP(K)*X(IX)
|
||||
IX = IX + INCX
|
||||
IY = IY + INCY
|
||||
70 CONTINUE
|
||||
Y(JY) = Y(JY) + TEMP1*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*AP(KK)
|
||||
K = KK + 1
|
||||
DO 90 I = J + 1,N
|
||||
Y(I) = Y(I) + TEMP1*AP(K)
|
||||
TEMP2 = TEMP2 + 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*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 + 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 SSPMV .
|
||||
*
|
||||
END
|
||||
@ -1,335 +0,0 @@
|
||||
SUBROUTINE STBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,K,LDA,N
|
||||
CHARACTER DIAG,TRANS,UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
REAL A(LDA,*),X(*)
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* STBMV performs one of the matrix-vector operations
|
||||
*
|
||||
* x := A*x, or x := A'*x,
|
||||
*
|
||||
* where x is an n element vector and A is an n by n unit, or non-unit,
|
||||
* upper or lower triangular band matrix, with ( k + 1 ) diagonals.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* UPLO - CHARACTER*1.
|
||||
* On entry, UPLO specifies whether the matrix is an upper or
|
||||
* lower triangular matrix as follows:
|
||||
*
|
||||
* UPLO = 'U' or 'u' A is an upper triangular matrix.
|
||||
*
|
||||
* UPLO = 'L' or 'l' A is a lower triangular matrix.
|
||||
*
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* TRANS - CHARACTER*1.
|
||||
* On entry, TRANS specifies the operation to be performed as
|
||||
* follows:
|
||||
*
|
||||
* TRANS = 'N' or 'n' x := A*x.
|
||||
*
|
||||
* TRANS = 'T' or 't' x := A'*x.
|
||||
*
|
||||
* TRANS = 'C' or 'c' x := A'*x.
|
||||
*
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* DIAG - CHARACTER*1.
|
||||
* On entry, DIAG specifies whether or not A is unit
|
||||
* triangular as follows:
|
||||
*
|
||||
* DIAG = 'U' or 'u' A is assumed to be unit triangular.
|
||||
*
|
||||
* DIAG = 'N' or 'n' A is not assumed to be unit
|
||||
* triangular.
|
||||
*
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* N - INTEGER.
|
||||
* On entry, N specifies the order of the matrix A.
|
||||
* N must be at least zero.
|
||||
* Unchanged on exit.
|
||||
*
|
||||
* 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 - REAL 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 - REAL 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* 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 ..
|
||||
REAL ZERO
|
||||
PARAMETER (ZERO=0.0E+0)
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
REAL TEMP
|
||||
INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L
|
||||
LOGICAL NOUNIT
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX,MIN
|
||||
* ..
|
||||
*
|
||||
* 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('STBMV ',INFO)
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible.
|
||||
*
|
||||
IF (N.EQ.0) RETURN
|
||||
*
|
||||
NOUNIT = LSAME(DIAG,'N')
|
||||
*
|
||||
* Set up the start point in X if the increment is not unity. This
|
||||
* will be ( N - 1 )*INCX too small for descending loops.
|
||||
*
|
||||
IF (INCX.LE.0) THEN
|
||||
KX = 1 - (N-1)*INCX
|
||||
ELSE IF (INCX.NE.1) THEN
|
||||
KX = 1
|
||||
END IF
|
||||
*
|
||||
* Start the operations. In this version the elements of A are
|
||||
* accessed sequentially with one pass through A.
|
||||
*
|
||||
IF (LSAME(TRANS,'N')) THEN
|
||||
*
|
||||
* Form x := A*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
KPLUS1 = K + 1
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 20 J = 1,N
|
||||
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)
|
||||
END IF
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 40 J = 1,N
|
||||
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)
|
||||
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
|
||||
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)
|
||||
END IF
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 80 J = N,1,-1
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
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)
|
||||
END IF
|
||||
JX = JX - INCX
|
||||
IF ((N-J).GE.K) KX = KX - INCX
|
||||
80 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
ELSE
|
||||
*
|
||||
* Form x := A'*x.
|
||||
*
|
||||
IF (LSAME(UPLO,'U')) THEN
|
||||
KPLUS1 = K + 1
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 100 J = N,1,-1
|
||||
TEMP = X(J)
|
||||
L = KPLUS1 - J
|
||||
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
|
||||
X(J) = TEMP
|
||||
100 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 120 J = N,1,-1
|
||||
TEMP = X(JX)
|
||||
KX = KX - INCX
|
||||
IX = KX
|
||||
L = KPLUS1 - J
|
||||
IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
|
||||
DO 110 I = J - 1,MAX(1,J-K),-1
|
||||
TEMP = TEMP + A(L+I,J)*X(IX)
|
||||
IX = IX - INCX
|
||||
110 CONTINUE
|
||||
X(JX) = TEMP
|
||||
JX = JX - INCX
|
||||
120 CONTINUE
|
||||
END IF
|
||||
ELSE
|
||||
IF (INCX.EQ.1) THEN
|
||||
DO 140 J = 1,N
|
||||
TEMP = X(J)
|
||||
L = 1 - J
|
||||
IF (NOUNIT) TEMP = TEMP*A(1,J)
|
||||
DO 130 I = J + 1,MIN(N,J+K)
|
||||
TEMP = TEMP + A(L+I,J)*X(I)
|
||||
130 CONTINUE
|
||||
X(J) = TEMP
|
||||
140 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 160 J = 1,N
|
||||
TEMP = X(JX)
|
||||
KX = KX + INCX
|
||||
IX = KX
|
||||
L = 1 - J
|
||||
IF (NOUNIT) TEMP = TEMP*A(1,J)
|
||||
DO 150 I = J + 1,MIN(N,J+K)
|
||||
TEMP = TEMP + A(L+I,J)*X(IX)
|
||||
IX = IX + INCX
|
||||
150 CONTINUE
|
||||
X(JX) = TEMP
|
||||
JX = JX + INCX
|
||||
160 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of STBMV .
|
||||
*
|
||||
END
|
||||
@ -1,310 +0,0 @@
|
||||
SUBROUTINE ZHBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
DOUBLE COMPLEX ALPHA,BETA
|
||||
INTEGER INCX,INCY,K,LDA,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE COMPLEX 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.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* Level 2 Blas routine.
|
||||
*
|
||||
* -- Written on 22-October-1986.
|
||||
* Jack Dongarra, Argonne National Lab.
|
||||
* Jeremy Du Croz, Nag Central Office.
|
||||
* Sven Hammarling, Nag Central Office.
|
||||
* Richard Hanson, Sandia National Labs.
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
DOUBLE COMPLEX ONE
|
||||
PARAMETER (ONE= (1.0D+0,0.0D+0))
|
||||
DOUBLE COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0D+0,0.0D+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
DOUBLE COMPLEX 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 DBLE,DCONJG,MAX,MIN
|
||||
* ..
|
||||
*
|
||||
* 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
|
||||
@ -1,272 +0,0 @@
|
||||
SUBROUTINE ZHPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
|
||||
* .. Scalar Arguments ..
|
||||
DOUBLE COMPLEX ALPHA,BETA
|
||||
INTEGER INCX,INCY,N
|
||||
CHARACTER UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE COMPLEX 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.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* Level 2 Blas routine.
|
||||
*
|
||||
* -- Written on 22-October-1986.
|
||||
* Jack Dongarra, Argonne National Lab.
|
||||
* Jeremy Du Croz, Nag Central Office.
|
||||
* Sven Hammarling, Nag Central Office.
|
||||
* Richard Hanson, Sandia National Labs.
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
DOUBLE COMPLEX ONE
|
||||
PARAMETER (ONE= (1.0D+0,0.0D+0))
|
||||
DOUBLE COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0D+0,0.0D+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
DOUBLE COMPLEX 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 DBLE,DCONJG
|
||||
* ..
|
||||
*
|
||||
* 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
|
||||
@ -1,366 +0,0 @@
|
||||
SUBROUTINE ZTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
|
||||
* .. Scalar Arguments ..
|
||||
INTEGER INCX,K,LDA,N
|
||||
CHARACTER DIAG,TRANS,UPLO
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE COMPLEX 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.
|
||||
*
|
||||
* Arguments
|
||||
* ==========
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* Further Details
|
||||
* ===============
|
||||
*
|
||||
* Level 2 Blas routine.
|
||||
*
|
||||
* -- Written on 22-October-1986.
|
||||
* Jack Dongarra, Argonne National Lab.
|
||||
* Jeremy Du Croz, Nag Central Office.
|
||||
* Sven Hammarling, Nag Central Office.
|
||||
* Richard Hanson, Sandia National Labs.
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
DOUBLE COMPLEX ZERO
|
||||
PARAMETER (ZERO= (0.0D+0,0.0D+0))
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
DOUBLE COMPLEX 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
|
||||
* ..
|
||||
*
|
||||
* 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
|
||||
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)
|
||||
END IF
|
||||
20 CONTINUE
|
||||
ELSE
|
||||
JX = KX
|
||||
DO 40 J = 1,N
|
||||
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)
|
||||
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
|
||||
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)
|
||||
END IF
|
||||
60 CONTINUE
|
||||
ELSE
|
||||
KX = KX + (N-1)*INCX
|
||||
JX = KX
|
||||
DO 80 J = N,1,-1
|
||||
IF (X(JX).NE.ZERO) THEN
|
||||
TEMP = X(JX)
|
||||
IX = KX
|
||||
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)
|
||||
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
|
||||
Loading…
x
Reference in New Issue
Block a user