subroutine dgemv ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) ! !******************************************************************************* ! !! DGEMV ! ! .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N character TRANS ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) ! .. ! ! Purpose ! ======= ! ! DGEMV performs one of the matrix-vector operations ! ! y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, ! ! where alpha and beta are scalars, x and y are vectors and A is an ! m by n matrix. ! ! Parameters ! ========== ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' y := alpha*A*x + beta*y. ! ! TRANS = 'T' or 't' y := alpha*A'*x + beta*y. ! ! TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix A. ! M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry, the leading m by n part of the array A must ! contain the matrix of coefficients. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, m ). ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of DIMENSION at least ! ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. ! Before entry, the incremented array X must contain the ! vector x. ! Unchanged on exit. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! BETA - 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 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. ! Before entry with BETA non-zero, the incremented array Y ! must contain the vector y. On exit, Y is overwritten by the ! updated vector y. ! ! INCY - INTEGER. ! On entry, INCY specifies the increment for the elements of ! Y. INCY must not be zero. ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 1 else if ( M < 0 ) then INFO = 2 else if ( N < 0 ) then INFO = 3 else if ( LDA < MAX( 1, M ) ) then INFO = 6 else if ( INCX == 0 ) then INFO = 8 else if ( INCY == 0 ) then INFO = 11 end if if ( INFO /= 0 ) then CALL XERBLA( 'DGEMV ', INFO ) return end if ! ! Quick return if possible. ! if ( ( M == 0 ) .or. ( N == 0 ) .or. & ( ( ALPHA == ZERO ).AND.( BETA == ONE ) ) ) then return end if ! ! Set LENX and LENY, the lengths of the vectors x and y, and set ! up the start points in X and Y. ! if ( LSAME( TRANS, 'N' ) ) then LENX = N LENY = M else LENX = M LENY = N end if if ( INCX > 0 ) then KX = 1 else KX = 1 - ( LENX - 1 )*INCX end if if ( INCY > 0 ) then KY = 1 else KY = 1 - ( LENY - 1 )*INCY end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through A. ! ! First form y := beta*y. ! if ( BETA /= ONE ) then if ( INCY == 1 ) then if ( BETA == ZERO ) then DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE else DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE end if else IY = KY if ( BETA == ZERO ) then DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE else DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE end if end if end if if ( ALPHA == ZERO ) then return end if if ( LSAME( TRANS, 'N' ) ) then ! ! Form y := alpha*A*x + y. ! JX = KX if ( INCY == 1 ) then DO 60, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE end if JX = JX + INCX 60 CONTINUE else DO 80, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE end if JX = JX + INCX 80 CONTINUE end if else ! ! Form y := alpha*A'*x + y. ! JY = KY if ( INCX == 1 ) then DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE else DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE end if end if return end subroutine dgbmv ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) ! !******************************************************************************* ! !! DGBMV ! DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, KL, KU, LDA, M, N character TRANS ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) ! .. ! ! Purpose ! ======= ! ! DGBMV performs one of the matrix-vector operations ! ! y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, ! ! where alpha and beta are scalars, x and y are vectors and A is an ! m by n band matrix, with kl sub-diagonals and ku super-diagonals. ! ! Parameters ! ========== ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' y := alpha*A*x + beta*y. ! ! TRANS = 'T' or 't' y := alpha*A'*x + beta*y. ! ! TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix A. ! M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! KL - INTEGER. ! On entry, KL specifies the number of sub-diagonals of the ! matrix A. KL must satisfy 0 .le. KL. ! Unchanged on exit. ! ! KU - INTEGER. ! On entry, KU specifies the number of super-diagonals of the ! matrix A. KU must satisfy 0 .le. KU. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry, the leading ( kl + ku + 1 ) by n part of the ! array A must contain the matrix of coefficients, supplied ! column by column, with the leading diagonal of the matrix in ! row ( ku + 1 ) of the array, the first super-diagonal ! starting at position 2 in row ku, the first sub-diagonal ! starting at position 1 in row ( ku + 2 ), and so on. ! Elements in the array A that do not correspond to elements ! in the band matrix (such as the top left ku by ku triangle) ! are not referenced. ! The following program segment will transfer a band matrix ! from conventional full matrix storage to band storage: ! ! DO 20, J = 1, N ! K = KU + 1 - J ! DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL ) ! A( K + I, J ) = matrix( I, J ) ! 10 CONTINUE ! 20 CONTINUE ! ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! ( kl + ku + 1 ). ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of DIMENSION at least ! ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. ! Before entry, the incremented array X must contain the ! vector x. ! Unchanged on exit. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! BETA - 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 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. ! Before entry, the incremented array Y must contain the ! vector y. On exit, Y is overwritten by the updated vector y. ! ! INCY - INTEGER. ! On entry, INCY specifies the increment for the elements of ! Y. INCY must not be zero. ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY integer LENX, LENY ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 1 else if ( M < 0 ) then INFO = 2 else if ( N < 0 ) then INFO = 3 else if ( KL < 0 ) then INFO = 4 else if ( KU < 0 ) then INFO = 5 else if ( LDA < ( KL + KU + 1 ) ) then INFO = 8 else if ( INCX == 0 ) then INFO = 10 else if ( INCY == 0 ) then INFO = 13 end if if ( INFO /= 0 ) then CALL XERBLA( 'DGBMV ', INFO ) return end if ! ! Quick return if possible. ! if ( ( M == 0 ) .or. ( N == 0 ) .or. & ( ( ALPHA == ZERO ).AND.( BETA == ONE ) ) ) then return end if ! ! Set LENX and LENY, the lengths of the vectors x and y, and set ! up the start points in X and Y. ! if ( LSAME( TRANS, 'N' ) ) then LENX = N LENY = M else LENX = M LENY = N end if if ( INCX > 0 ) then KX = 1 else KX = 1 - ( LENX - 1 )*INCX end if if ( INCY > 0 ) then KY = 1 else KY = 1 - ( LENY - 1 )*INCY end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through the band part of A. ! ! First form y := beta*y. ! if ( BETA /= ONE ) then if ( INCY == 1 ) then if ( BETA == ZERO ) then DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE else DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE end if else IY = KY if ( BETA == ZERO ) then DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE else DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE end if end if end if if ( ALPHA == ZERO ) then return end if KUP1 = KU + 1 if ( LSAME( TRANS, 'N' ) ) then ! ! Form y := alpha*A*x + y. ! JX = KX if ( INCY == 1 ) then DO 60, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) K = KUP1 - J DO 50, I = MAX( 1, J - KU ), MIN( M, J + KL ) Y( I ) = Y( I ) + TEMP*A( K + I, J ) 50 CONTINUE end if JX = JX + INCX 60 CONTINUE else DO 80, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) IY = KY K = KUP1 - J DO 70, I = MAX( 1, J - KU ), MIN( M, J + KL ) Y( IY ) = Y( IY ) + TEMP*A( K + I, J ) IY = IY + INCY 70 CONTINUE end if JX = JX + INCX if ( J > KU ) then KY = KY + INCY end if 80 CONTINUE end if else ! ! Form y := alpha*A'*x + y. ! JY = KY if ( INCX == 1 ) then DO 100, J = 1, N TEMP = ZERO K = KUP1 - J DO 90, I = MAX( 1, J - KU ), MIN( M, J + KL ) TEMP = TEMP + A( K + I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE else DO 120, J = 1, N TEMP = ZERO IX = KX K = KUP1 - J DO 110, I = MAX( 1, J - KU ), MIN( M, J + KL ) TEMP = TEMP + A( K + I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY if ( J > KU ) then KX = KX + INCX end if 120 CONTINUE end if end if return end subroutine dsymv ( UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) ! !******************************************************************************* ! !! DSYMV ! DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, N character UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) ! .. ! ! Purpose ! ======= ! ! DSYMV 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. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array A is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of A ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of A ! is to be referenced. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array A must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of A is not referenced. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array A must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of A is not referenced. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, n ). ! 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. ! ! ! 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, KX, KY ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. & .NOT.LSAME( UPLO, 'L' ) ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( LDA < MAX( 1, N ) ) then INFO = 5 else if ( INCX == 0 ) then INFO = 7 else if ( INCY == 0 ) then INFO = 10 end if if ( INFO /= 0 ) then CALL XERBLA( 'DSYMV ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. ( ( ALPHA == ZERO ).AND.( BETA == ONE ) ) ) then return end if ! ! Set up the start points in X and Y. ! if ( INCX > 0 ) then KX = 1 else KX = 1 - ( N - 1 )*INCX end if if ( INCY > 0 ) then KY = 1 else KY = 1 - ( N - 1 )*INCY end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through the triangular part ! of A. ! ! First form y := beta*y. ! if ( BETA /= ONE ) then if ( INCY == 1 ) then if ( BETA == 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 == 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 == ZERO ) then return end if if ( LSAME( UPLO, 'U' ) ) then ! ! Form y when A is stored in upper triangle. ! if ( ( INCX == 1 ).AND.( INCY == 1 ) ) then DO 60, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO DO 50, I = 1, J - 1 Y( I ) = Y( I ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + A( I, J )*X( I ) 50 CONTINUE Y( J ) = Y( J ) + TEMP1*A( J, J ) + ALPHA*TEMP2 60 CONTINUE else JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY DO 70, I = 1, J - 1 Y( IY ) = Y( IY ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + A( I, J )*X( IX ) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y( JY ) = Y( JY ) + TEMP1*A( J, J ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 80 CONTINUE end if else ! ! Form y when A is stored in lower triangle. ! if ( ( INCX == 1 ).AND.( INCY == 1 ) ) then DO 100, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO Y( J ) = Y( J ) + TEMP1*A( J, J ) DO 90, I = J + 1, N Y( I ) = Y( I ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + A( I, J )*X( I ) 90 CONTINUE Y( J ) = Y( J ) + ALPHA*TEMP2 100 CONTINUE else JX = KX JY = KY DO 120, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO Y( JY ) = Y( JY ) + TEMP1*A( J, J ) IX = JX IY = JY DO 110, I = J + 1, N IX = IX + INCX IY = IY + INCY Y( IY ) = Y( IY ) + TEMP1*A( I, J ) TEMP2 = TEMP2 + A( I, J )*X( IX ) 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 120 CONTINUE end if end if return end subroutine dsbmv ( UPLO, N, K, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) ! !******************************************************************************* ! !! DSBMV ! 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. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the band matrix A is being supplied as ! follows: ! ! UPLO = 'U' or 'u' The upper triangular part of A is ! being supplied. ! ! UPLO = 'L' or 'l' The lower triangular part of A is ! being supplied. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry, K specifies the number of super-diagonals of the ! matrix A. K must satisfy 0 .le. K. ! Unchanged on exit. ! ! ALPHA - 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 ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. & .NOT.LSAME( UPLO, 'L' ) ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( K < 0 ) then INFO = 3 else if ( LDA < ( K + 1 ) ) then INFO = 6 else if ( INCX == 0 ) then INFO = 8 else if ( INCY == 0 ) then INFO = 11 end if if ( INFO /= 0 ) then CALL XERBLA( 'DSBMV ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. ( ( ALPHA == ZERO ).AND.( BETA == ONE ) ) ) then return end if ! ! Set up the start points in X and Y. ! if ( INCX > 0 ) then KX = 1 else KX = 1 - ( N - 1 )*INCX end if if ( INCY > 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 /= ONE ) then if ( INCY == 1 ) then if ( BETA == 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 == 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 == ZERO ) then return end if if ( LSAME( UPLO, 'U' ) ) then ! ! Form y when upper triangle of A is stored. ! KPLUS1 = K + 1 if ( ( INCX == 1 ).AND.( INCY == 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 > 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 == 1 ).AND.( INCY == 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 subroutine dspmv ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) ! !******************************************************************************* ! !! DSPMV ! 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. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the matrix A is supplied in the packed ! array AP as follows: ! ! UPLO = 'U' or 'u' The upper triangular part of A is ! supplied in AP. ! ! UPLO = 'L' or 'l' The lower triangular part of A is ! supplied in AP. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! 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. ! ! ! 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 ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. & .NOT.LSAME( UPLO, 'L' ) ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( INCX == 0 ) then INFO = 6 else if ( INCY == 0 ) then INFO = 9 end if if ( INFO /= 0 ) then CALL XERBLA( 'DSPMV ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. ( ( ALPHA == ZERO ).AND.( BETA.EQ.ONE ) ) ) then return end if ! ! Set up the start points in X and Y. ! if ( INCX > 0 ) then KX = 1 else KX = 1 - ( N - 1 )*INCX end if if ( INCY > 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 /= ONE ) then if ( INCY == 1 ) then if ( BETA == 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 == 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 == ZERO ) then return end if KK = 1 if ( LSAME( UPLO, 'U' ) ) then ! ! Form y when AP contains the upper triangle. ! if ( ( INCX == 1 ).AND.( INCY == 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 == 1 ).AND.( INCY == 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 subroutine dtrmv ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) ! !******************************************************************************* ! !! DTRMV ! INTEGER INCX, LDA, N character DIAG, TRANS, UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DTRMV 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 matrix. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' x := A*x. ! ! TRANS = 'T' or 't' x := A'*x. ! ! TRANS = 'C' or 'c' x := A'*x. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit ! triangular as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array A must contain the upper ! triangular matrix and the strictly lower triangular part of ! A is not referenced. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array A must contain the lower ! triangular matrix and the strictly upper triangular part of ! A is not referenced. ! Note that when DIAG = 'U' or 'u', the diagonal elements of ! A are not referenced either, but are assumed to be unity. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, n ). ! Unchanged on exit. ! ! X - 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. ! ! ! 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, KX LOGICAL NOUNIT ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. & .NOT.LSAME( UPLO , 'L' ) ) then INFO = 1 else if ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 2 else if ( .NOT.LSAME( DIAG , 'U' ).AND. & .NOT.LSAME( DIAG , 'N' ) ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( LDA < MAX( 1, N ) ) then INFO = 6 else if ( INCX == 0 ) then INFO = 8 end if if ( INFO /= 0 ) then CALL XERBLA( 'DTRMV ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if 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 /= 1 ) then KX = 1 end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through A. ! if ( LSAME( TRANS, 'N' ) ) then ! ! Form x := A*x. ! if ( LSAME( UPLO, 'U' ) ) then if ( INCX == 1 ) then DO 20, J = 1, N if ( X( J ) /= ZERO ) then TEMP = X( J ) DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*A( I, J ) 10 CONTINUE if ( NOUNIT ) then X( J ) = X( J )*A( J, J ) end if end if 20 CONTINUE else JX = KX DO 40, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = X( JX ) IX = KX DO 30, I = 1, J - 1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX + INCX 30 CONTINUE if ( NOUNIT ) then X( JX ) = X( JX )*A( J, J ) end if end if JX = JX + INCX 40 CONTINUE end if else if ( INCX == 1 ) then DO 60, J = N, 1, -1 if ( X( J ) /= ZERO ) then TEMP = X( J ) DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*A( I, J ) 50 CONTINUE if ( NOUNIT ) then X( J ) = X( J )*A( J, J ) end if end if 60 CONTINUE else KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 if ( X( JX ) /= ZERO ) then TEMP = X( JX ) IX = KX DO 70, I = N, J + 1, -1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX - INCX 70 CONTINUE if ( NOUNIT ) then X( JX ) = X( JX )*A( J, J ) end if end if JX = JX - INCX 80 CONTINUE end if end if else ! ! Form x := A'*x. ! if ( LSAME( UPLO, 'U' ) ) then if ( INCX == 1 ) then DO 100, J = N, 1, -1 TEMP = X( J ) if ( NOUNIT ) then TEMP = TEMP*A( J, J ) end if DO 90, I = J - 1, 1, -1 TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE else JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX if ( NOUNIT ) then TEMP = TEMP*A( J, J ) end if DO 110, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + A( I, J )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE end if else if ( INCX == 1 ) then DO 140, J = 1, N TEMP = X( J ) if ( NOUNIT ) then TEMP = TEMP*A( J, J ) end if DO 130, I = J + 1, N TEMP = TEMP + A( I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE else JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX if ( NOUNIT ) then TEMP = TEMP*A( J, J ) end if DO 150, I = J + 1, N IX = IX + INCX TEMP = TEMP + A( I, J )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE end if end if end if return end subroutine dtbmv ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) ! !******************************************************************************* ! !! DTBMV ! 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. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' x := A*x. ! ! TRANS = 'T' or 't' x := A'*x. ! ! TRANS = 'C' or 'c' x := 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. ! ! ! 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 ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. & .NOT.LSAME( UPLO , 'L' ) ) then INFO = 1 else if ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 2 else if ( .NOT.LSAME( DIAG , 'U' ).AND. & .NOT.LSAME( DIAG , 'N' ) ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( K < 0 ) then INFO = 5 else if ( LDA < ( K + 1 ) ) then INFO = 7 else if ( INCX == 0 ) then INFO = 9 end if if ( INFO /= 0 ) then CALL XERBLA( 'DTBMV ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if ! 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 /= 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 == 1 ) then DO 20, J = 1, N if ( X( J ) /= 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 ) then X( J ) = X( J )*A( KPLUS1, J ) end if end if 20 CONTINUE else JX = KX DO 40, J = 1, N if ( X( JX ) /= 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 ) then X( JX ) = X( JX )*A( KPLUS1, J ) end if end if JX = JX + INCX if ( J > K ) then KX = KX + INCX end if 40 CONTINUE end if else if ( INCX == 1 ) then DO 60, J = N, 1, -1 if ( X( J ) /= 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 ) then X( J ) = X( J )*A( 1, J ) end if end if 60 CONTINUE else KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 if ( X( JX ) /= 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 ) then X( JX ) = X( JX )*A( 1, J ) end if end if JX = JX - INCX if ( ( N - J ).GE.K ) then KX = KX - INCX end if 80 CONTINUE end if end if else ! ! Form x := A'*x. ! if ( LSAME( UPLO, 'U' ) ) then KPLUS1 = K + 1 if ( INCX == 1 ) then DO 100, J = N, 1, -1 TEMP = X( J ) L = KPLUS1 - J if ( NOUNIT ) then TEMP = TEMP*A( KPLUS1, J ) end if 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 ) then TEMP = TEMP*A( KPLUS1, J ) end if 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 == 1 ) then DO 140, J = 1, N TEMP = X( J ) L = 1 - J if ( NOUNIT ) then TEMP = TEMP*A( 1, J ) end if 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 ) then TEMP = TEMP*A( 1, J ) end if 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 subroutine dtpmv ( UPLO, TRANS, DIAG, N, AP, X, INCX ) ! !******************************************************************************* ! !! DTPMV ! INTEGER INCX, N character DIAG, TRANS, UPLO ! .. Array Arguments .. DOUBLE PRECISION AP( * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DTPMV 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 matrix, supplied in packed form. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' x := A*x. ! ! TRANS = 'T' or 't' x := A'*x. ! ! TRANS = 'C' or 'c' x := A'*x. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit ! triangular as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! AP - 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 matrix packed sequentially, ! column by column, so that AP( 1 ) contains a( 1, 1 ), ! AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) ! respectively, and so on. ! Before entry with UPLO = 'L' or 'l', the array AP must ! contain the lower triangular matrix packed sequentially, ! column by column, so that AP( 1 ) contains a( 1, 1 ), ! AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) ! respectively, and so on. ! Note that when DIAG = 'U' or 'u', the diagonal elements of ! A are not referenced, but are assumed to be unity. ! Unchanged on exit. ! ! X - 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. ! ! ! 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, K, KK, KX LOGICAL NOUNIT ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. & .NOT.LSAME( UPLO , 'L' ) ) then INFO = 1 else if ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 2 else if ( .NOT.LSAME( DIAG , 'U' ).AND. & .NOT.LSAME( DIAG , 'N' ) ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( INCX == 0 ) then INFO = 7 end if if ( INFO /= 0 ) then CALL XERBLA( 'DTPMV ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if ! 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 /= 1 ) then KX = 1 end if ! ! Start the operations. In this version the elements of AP are ! accessed sequentially with one pass through AP. ! if ( LSAME( TRANS, 'N' ) ) then ! ! Form x:= A*x. ! if ( LSAME( UPLO, 'U' ) ) then KK =1 if ( INCX == 1 ) then DO 20, J = 1, N if ( X( J ) /= ZERO ) then TEMP = X( J ) K = KK DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*AP( K ) K = K + 1 10 CONTINUE if ( NOUNIT ) then X( J ) = X( J )*AP( KK + J - 1 ) end if end if KK = KK + J 20 CONTINUE else JX = KX DO 40, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = X( JX ) IX = KX DO 30, K = KK, KK + J - 2 X( IX ) = X( IX ) + TEMP*AP( K ) IX = IX + INCX 30 CONTINUE if ( NOUNIT ) then X( JX ) = X( JX )*AP( KK + J - 1 ) end if end if JX = JX + INCX KK = KK + J 40 CONTINUE end if else KK = ( N*( N + 1 ) )/2 if ( INCX == 1 ) then DO 60, J = N, 1, -1 if ( X( J ) /= ZERO ) then TEMP = X( J ) K = KK DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*AP( K ) K = K - 1 50 CONTINUE if ( NOUNIT ) then X( J ) = X( J )*AP( KK - N + J ) end if end if KK = KK - ( N - J + 1 ) 60 CONTINUE else KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 if ( X( JX ) /= ZERO ) then TEMP = X( JX ) IX = KX DO 70, K = KK, KK - ( N - ( J + 1 ) ), -1 X( IX ) = X( IX ) + TEMP*AP( K ) IX = IX - INCX 70 CONTINUE if ( NOUNIT ) then X( JX ) = X( JX )*AP( KK - N + J ) end if end if JX = JX - INCX KK = KK - ( N - J + 1 ) 80 CONTINUE end if end if else ! ! Form x := A'*x. ! if ( LSAME( UPLO, 'U' ) ) then KK = ( N*( N + 1 ) )/2 if ( INCX == 1 ) then DO 100, J = N, 1, -1 TEMP = X( J ) if ( NOUNIT ) then TEMP = TEMP*AP( KK ) end if K = KK - 1 DO 90, I = J - 1, 1, -1 TEMP = TEMP + AP( K )*X( I ) K = K - 1 90 CONTINUE X( J ) = TEMP KK = KK - J 100 CONTINUE else JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX if ( NOUNIT ) then TEMP = TEMP*AP( KK ) end if DO 110, K = KK - 1, KK - J + 1, -1 IX = IX - INCX TEMP = TEMP + AP( K )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX KK = KK - J 120 CONTINUE end if else KK = 1 if ( INCX == 1 ) then DO 140, J = 1, N TEMP = X( J ) if ( NOUNIT ) then TEMP = TEMP*AP( KK ) end if K = KK + 1 DO 130, I = J + 1, N TEMP = TEMP + AP( K )*X( I ) K = K + 1 130 CONTINUE X( J ) = TEMP KK = KK + ( N - J + 1 ) 140 CONTINUE else JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX if ( NOUNIT ) then TEMP = TEMP*AP( KK ) end if DO 150, K = KK + 1, KK + N - J IX = IX + INCX TEMP = TEMP + AP( K )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX KK = KK + ( N - J + 1 ) 160 CONTINUE end if end if end if return end subroutine dtrsv ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) ! !******************************************************************************* ! !! DTRSV ! INTEGER INCX, LDA, N character DIAG, TRANS, UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DTRSV solves one of the systems of equations ! ! A*x = b, or A'*x = b, ! ! where b and x are n element vectors and A is an n by n unit, or ! non-unit, upper or lower triangular matrix. ! ! No test for singularity or near-singularity is included in this ! routine. Such tests must be performed before calling this routine. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the equations to be solved as ! follows: ! ! TRANS = 'N' or 'n' A*x = b. ! ! TRANS = 'T' or 't' A'*x = b. ! ! TRANS = 'C' or 'c' A'*x = b. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit ! triangular as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array A must contain the upper ! triangular matrix and the strictly lower triangular part of ! A is not referenced. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array A must contain the lower ! triangular matrix and the strictly upper triangular part of ! A is not referenced. ! Note that when DIAG = 'U' or 'u', the diagonal elements of ! A are not referenced either, but are assumed to be unity. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, n ). ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of dimension at least ! ( 1 + ( n - 1 )*abs( INCX ) ). ! Before entry, the incremented array X must contain the n ! element right-hand side vector b. On exit, X is overwritten ! with the solution vector x. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. & .NOT.LSAME( UPLO , 'L' ) ) then INFO = 1 else if ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 2 else if ( .NOT.LSAME( DIAG , 'U' ).AND. & .NOT.LSAME( DIAG , 'N' ) ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( LDA < MAX( 1, N ) ) then INFO = 6 else if ( INCX == 0 ) then INFO = 8 end if if ( INFO /= 0 ) then CALL XERBLA( 'DTRSV ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if ! 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 /= 1 ) then KX = 1 end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through A. ! if ( LSAME( TRANS, 'N' ) ) then ! ! Form x := inv( A )*x. ! if ( LSAME( UPLO, 'U' ) ) then if ( INCX == 1 ) then DO 20, J = N, 1, -1 if ( X( J ) /= ZERO ) then if ( NOUNIT ) then X( J ) = X( J )/A( J, J ) end if TEMP = X( J ) DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*A( I, J ) 10 CONTINUE end if 20 CONTINUE else JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 if ( X( JX ) /= ZERO ) then if ( NOUNIT ) then X( JX ) = X( JX )/A( J, J ) end if TEMP = X( JX ) IX = JX DO 30, I = J - 1, 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 30 CONTINUE end if JX = JX - INCX 40 CONTINUE end if else if ( INCX == 1 ) then DO 60, J = 1, N if ( X( J ) /= ZERO ) then if ( NOUNIT ) then X( J ) = X( J )/A( J, J ) end if TEMP = X( J ) DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*A( I, J ) 50 CONTINUE end if 60 CONTINUE else JX = KX DO 80, J = 1, N if ( X( JX ) /= ZERO ) then if ( NOUNIT ) then X( JX ) = X( JX )/A( J, J ) end if TEMP = X( JX ) IX = JX DO 70, I = J + 1, N IX = IX + INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 70 CONTINUE end if JX = JX + INCX 80 CONTINUE end if end if else ! ! Form x := inv( A' )*x. ! if ( LSAME( UPLO, 'U' ) ) then if ( INCX == 1 ) then DO 100, J = 1, N TEMP = X( J ) DO 90, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( I ) 90 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( J, J ) end if X( J ) = TEMP 100 CONTINUE else JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX DO 110, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( J, J ) end if X( JX ) = TEMP JX = JX + INCX 120 CONTINUE end if else if ( INCX == 1 ) then DO 140, J = N, 1, -1 TEMP = X( J ) DO 130, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( I ) 130 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( J, J ) end if X( J ) = TEMP 140 CONTINUE else KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX DO 150, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX - INCX 150 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( J, J ) end if X( JX ) = TEMP JX = JX - INCX 160 CONTINUE end if end if end if return end subroutine dtbsv ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) ! !******************************************************************************* ! !! DTBSV ! INTEGER INCX, K, LDA, N character DIAG, TRANS, UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DTBSV solves one of the systems of equations ! ! A*x = b, or A'*x = b, ! ! where b and x are n element vectors and A is an n by n unit, or ! non-unit, upper or lower triangular band matrix, with ( k + 1 ) ! diagonals. ! ! No test for singularity or near-singularity is included in this ! routine. Such tests must be performed before calling this routine. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the equations to be solved as ! follows: ! ! TRANS = 'N' or 'n' A*x = b. ! ! TRANS = 'T' or 't' A'*x = b. ! ! TRANS = 'C' or 'c' A'*x = b. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit ! triangular as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with UPLO = 'U' or 'u', K specifies the number of ! super-diagonals of the matrix A. ! On entry with UPLO = 'L' or 'l', K specifies the number of ! sub-diagonals of the matrix A. ! K must satisfy 0 .le. K. ! Unchanged on exit. ! ! A - 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 right-hand side vector b. On exit, X is overwritten ! with the solution vector x. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. 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 ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. & .NOT.LSAME( UPLO , 'L' ) ) then INFO = 1 else if ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 2 else if ( .NOT.LSAME( DIAG , 'U' ).AND. & .NOT.LSAME( DIAG , 'N' ) ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( K < 0 ) then INFO = 5 else if ( LDA < ( K + 1 ) ) then INFO = 7 else if ( INCX == 0 ) then INFO = 9 end if if ( INFO /= 0 ) then CALL XERBLA( 'DTBSV ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if ! 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 /= 1 ) then KX = 1 end if ! ! Start the operations. In this version the elements of A are ! accessed by sequentially with one pass through A. ! if ( LSAME( TRANS, 'N' ) ) then ! ! Form x := inv( A )*x. ! if ( LSAME( UPLO, 'U' ) ) then KPLUS1 = K + 1 if ( INCX == 1 ) then DO 20, J = N, 1, -1 if ( X( J ) /= ZERO ) then L = KPLUS1 - J if ( NOUNIT ) then X( J ) = X( J )/A( KPLUS1, J ) end if TEMP = X( J ) DO 10, I = J - 1, MAX( 1, J - K ), -1 X( I ) = X( I ) - TEMP*A( L + I, J ) 10 CONTINUE end if 20 CONTINUE else KX = KX + ( N - 1 )*INCX JX = KX DO 40, J = N, 1, -1 KX = KX - INCX if ( X( JX ) /= ZERO ) then IX = KX L = KPLUS1 - J if ( NOUNIT ) then X( JX ) = X( JX )/A( KPLUS1, J ) end if TEMP = X( JX ) DO 30, I = J - 1, MAX( 1, J - K ), -1 X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX - INCX 30 CONTINUE end if JX = JX - INCX 40 CONTINUE end if else if ( INCX == 1 ) then DO 60, J = 1, N if ( X( J ) /= ZERO ) then L = 1 - J if ( NOUNIT ) then X( J ) = X( J )/A( 1, J ) end if TEMP = X( J ) DO 50, I = J + 1, MIN( N, J + K ) X( I ) = X( I ) - TEMP*A( L + I, J ) 50 CONTINUE end if 60 CONTINUE else JX = KX DO 80, J = 1, N KX = KX + INCX if ( X( JX ) /= ZERO ) then IX = KX L = 1 - J if ( NOUNIT ) then X( JX ) = X( JX )/A( 1, J ) end if TEMP = X( JX ) DO 70, I = J + 1, MIN( N, J + K ) X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX + INCX 70 CONTINUE end if JX = JX + INCX 80 CONTINUE end if end if else ! ! Form x := inv( A')*x. ! if ( LSAME( UPLO, 'U' ) ) then KPLUS1 = K + 1 if ( INCX == 1 ) then DO 100, J = 1, N TEMP = X( J ) L = KPLUS1 - J DO 90, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( I ) 90 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( KPLUS1, J ) end if X( J ) = TEMP 100 CONTINUE else JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX L = KPLUS1 - J DO 110, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX + INCX 110 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( KPLUS1, J ) end if X( JX ) = TEMP JX = JX + INCX if ( J > K ) then KX = KX + INCX end if 120 CONTINUE end if else if ( INCX == 1 ) then DO 140, J = N, 1, -1 TEMP = X( J ) L = 1 - J DO 130, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( I ) 130 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( 1, J ) end if X( J ) = TEMP 140 CONTINUE else KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX L = 1 - J DO 150, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX - INCX 150 CONTINUE if ( NOUNIT ) then TEMP = TEMP/A( 1, J ) end if X( JX ) = TEMP JX = JX - INCX if ( ( N - J ).GE.K ) then KX = KX - INCX end if 160 CONTINUE end if end if end if return end subroutine dtpsv ( UPLO, TRANS, DIAG, N, AP, X, INCX ) ! !******************************************************************************* ! !! DTPSV ! INTEGER INCX, N character DIAG, TRANS, UPLO ! .. Array Arguments .. DOUBLE PRECISION AP( * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DTPSV solves one of the systems of equations ! ! A*x = b, or A'*x = b, ! ! where b and x are n element vectors and A is an n by n unit, or ! non-unit, upper or lower triangular matrix, supplied in packed form. ! ! No test for singularity or near-singularity is included in this ! routine. Such tests must be performed before calling this routine. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the equations to be solved as ! follows: ! ! TRANS = 'N' or 'n' A*x = b. ! ! TRANS = 'T' or 't' A'*x = b. ! ! TRANS = 'C' or 'c' A'*x = b. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit ! triangular as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! AP - 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 matrix packed sequentially, ! column by column, so that AP( 1 ) contains a( 1, 1 ), ! AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) ! respectively, and so on. ! Before entry with UPLO = 'L' or 'l', the array AP must ! contain the lower triangular matrix packed sequentially, ! column by column, so that AP( 1 ) contains a( 1, 1 ), ! AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) ! respectively, and so on. ! Note that when DIAG = 'U' or 'u', the diagonal elements of ! A are not referenced, but are assumed to be unity. ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of dimension at least ! ( 1 + ( n - 1 )*abs( INCX ) ). ! Before entry, the incremented array X must contain the n ! element right-hand side vector b. On exit, X is overwritten ! with the solution vector x. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX LOGICAL NOUNIT ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. & .NOT.LSAME( UPLO , 'L' ) ) then INFO = 1 else if ( .NOT.LSAME( TRANS, 'N' ).AND. & .NOT.LSAME( TRANS, 'T' ).AND. & .NOT.LSAME( TRANS, 'C' ) ) then INFO = 2 else if ( .NOT.LSAME( DIAG , 'U' ).AND. & .NOT.LSAME( DIAG , 'N' ) ) then INFO = 3 else if ( N < 0 ) then INFO = 4 else if ( INCX == 0 ) then INFO = 7 end if if ( INFO /= 0 ) then CALL XERBLA( 'DTPSV ', INFO ) return end if ! ! Quick return if possible. ! if ( N == 0 ) then return end if 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 /= 1 ) then KX = 1 end if ! ! Start the operations. In this version the elements of AP are ! accessed sequentially with one pass through AP. ! if ( LSAME( TRANS, 'N' ) ) then ! ! Form x := inv( A )*x. ! if ( LSAME( UPLO, 'U' ) ) then KK = ( N*( N + 1 ) )/2 if ( INCX == 1 ) then DO 20, J = N, 1, -1 if ( X( J ) /= ZERO ) then if ( NOUNIT ) then X( J ) = X( J )/AP( KK ) end if TEMP = X( J ) K = KK - 1 DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*AP( K ) K = K - 1 10 CONTINUE end if KK = KK - J 20 CONTINUE else JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 if ( X( JX ) /= ZERO ) then if ( NOUNIT ) then X( JX ) = X( JX )/AP( KK ) end if TEMP = X( JX ) IX = JX DO 30, K = KK - 1, KK - J + 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*AP( K ) 30 CONTINUE end if JX = JX - INCX KK = KK - J 40 CONTINUE end if else KK = 1 if ( INCX == 1 ) then DO 60, J = 1, N if ( X( J ) /= ZERO ) then if ( NOUNIT ) then X( J ) = X( J )/AP( KK ) end if TEMP = X( J ) K = KK + 1 DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*AP( K ) K = K + 1 50 CONTINUE end if KK = KK + ( N - J + 1 ) 60 CONTINUE else JX = KX DO 80, J = 1, N if ( X( JX ) /= ZERO ) then if ( NOUNIT ) then X( JX ) = X( JX )/AP( KK ) end if TEMP = X( JX ) IX = JX DO 70, K = KK + 1, KK + N - J IX = IX + INCX X( IX ) = X( IX ) - TEMP*AP( K ) 70 CONTINUE end if JX = JX + INCX KK = KK + ( N - J + 1 ) 80 CONTINUE end if end if else ! ! Form x := inv( A' )*x. ! if ( LSAME( UPLO, 'U' ) ) then KK = 1 if ( INCX == 1 ) then DO 100, J = 1, N TEMP = X( J ) K = KK DO 90, I = 1, J - 1 TEMP = TEMP - AP( K )*X( I ) K = K + 1 90 CONTINUE if ( NOUNIT ) then TEMP = TEMP/AP( KK + J - 1 ) end if X( J ) = TEMP KK = KK + J 100 CONTINUE else JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX DO 110, K = KK, KK + J - 2 TEMP = TEMP - AP( K )*X( IX ) IX = IX + INCX 110 CONTINUE if ( NOUNIT ) then TEMP = TEMP/AP( KK + J - 1 ) end if X( JX ) = TEMP JX = JX + INCX KK = KK + J 120 CONTINUE end if else KK = ( N*( N + 1 ) )/2 if ( INCX == 1 ) then DO 140, J = N, 1, -1 TEMP = X( J ) K = KK DO 130, I = N, J + 1, -1 TEMP = TEMP - AP( K )*X( I ) K = K - 1 130 CONTINUE if ( NOUNIT ) then TEMP = TEMP/AP( KK - N + J ) end if X( J ) = TEMP KK = KK - ( N - J + 1 ) 140 CONTINUE else KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX DO 150, K = KK, KK - ( N - ( J + 1 ) ), -1 TEMP = TEMP - AP( K )*X( IX ) IX = IX - INCX 150 CONTINUE if ( NOUNIT ) then TEMP = TEMP/AP( KK - N + J ) end if X( JX ) = TEMP JX = JX - INCX KK = KK - (N - J + 1 ) 160 CONTINUE end if end if end if return end subroutine dger ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) ! !******************************************************************************* ! !! DGER ! DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, M, N ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) ! .. ! ! Purpose ! ======= ! ! DGER performs the rank 1 operation ! ! A := alpha*x*y' + A, ! ! where alpha is a scalar, x is an m element vector, y is an n element ! vector and A is an m by n matrix. ! ! Parameters ! ========== ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix A. ! M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of dimension at least ! ( 1 + ( m - 1 )*abs( INCX ) ). ! Before entry, the incremented array X must contain the m ! element vector x. ! Unchanged on exit. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! Y - 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. ! Unchanged on exit. ! ! INCY - INTEGER. ! On entry, INCY specifies the increment for the elements of ! Y. INCY must not be zero. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry, the leading m by n part of the array A must ! contain the matrix of coefficients. On exit, A is ! overwritten by the updated matrix. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, m ). ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JY, KX ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( M < 0 ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( INCX == 0 ) then INFO = 5 else if ( INCY == 0 ) then INFO = 7 else if ( LDA < MAX( 1, M ) ) then INFO = 9 end if if ( INFO /= 0 ) then CALL XERBLA( 'DGER ', INFO ) return end if ! ! Quick return if possible. ! if ( ( M == 0 ) .or. ( N == 0 ) .or. ( ALPHA.EQ.ZERO ) ) then return end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through A. ! if ( INCY > 0 ) then JY = 1 else JY = 1 - ( N - 1 )*INCY end if if ( INCX == 1 ) then DO 20, J = 1, N if ( Y( JY ) /= ZERO ) then TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE end if JY = JY + INCY 20 CONTINUE else if ( INCX > 0 ) then KX = 1 else KX = 1 - ( M - 1 )*INCX end if DO 40, J = 1, N if ( Y( JY ) /= ZERO ) then TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE end if JY = JY + INCY 40 CONTINUE end if return end subroutine dsyr ( UPLO, N, ALPHA, X, INCX, A, LDA ) ! !******************************************************************************* ! !! DSYR ! DOUBLE PRECISION ALPHA INTEGER INCX, LDA, N character UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DSYR performs the symmetric rank 1 operation ! ! A := alpha*x*x' + A, ! ! where alpha is a real scalar, x is an n element vector and A is an ! n by n symmetric matrix. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array A is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of A ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of A ! is to be referenced. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! X - 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. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array A must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of A is not referenced. On exit, the ! upper triangular part of the array A is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array A must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of A is not referenced. On exit, the ! lower triangular part of the array A is overwritten by the ! lower triangular part of the updated matrix. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, n ). ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. & .NOT.LSAME( UPLO, 'L' ) ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( INCX == 0 ) then INFO = 5 else if ( LDA < MAX( 1, N ) ) then INFO = 7 end if if ( INFO /= 0 ) then CALL XERBLA( 'DSYR ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. ( ALPHA == ZERO ) ) then return end if ! ! Set the start point in X if the increment is not unity. ! if ( INCX.LE.0 ) then KX = 1 - ( N - 1 )*INCX else if ( INCX /= 1 ) then KX = 1 end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through the triangular part ! of A. ! if ( LSAME( UPLO, 'U' ) ) then ! ! Form A when A is stored in upper triangle. ! if ( INCX == 1 ) then DO 20, J = 1, N if ( X( J ) /= ZERO ) then TEMP = ALPHA*X( J ) DO 10, I = 1, J A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE end if 20 CONTINUE else JX = KX DO 40, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) IX = KX DO 30, I = 1, J A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE end if JX = JX + INCX 40 CONTINUE end if else ! ! Form A when A is stored in lower triangle. ! if ( INCX == 1 ) then DO 60, J = 1, N if ( X( J ) /= ZERO ) then TEMP = ALPHA*X( J ) DO 50, I = J, N A( I, J ) = A( I, J ) + X( I )*TEMP 50 CONTINUE end if 60 CONTINUE else JX = KX DO 80, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) IX = JX DO 70, I = J, N A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 70 CONTINUE end if JX = JX + INCX 80 CONTINUE end if end if return end subroutine dspr ( UPLO, N, ALPHA, X, INCX, AP ) ! !******************************************************************************* ! !! DSPR ! DOUBLE PRECISION ALPHA INTEGER INCX, N character UPLO ! .. Array Arguments .. DOUBLE PRECISION AP( * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DSPR performs the symmetric rank 1 operation ! ! A := alpha*x*x' + A, ! ! where alpha is a real scalar, x is an n element vector and A is an ! n by n symmetric matrix, supplied in packed form. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the matrix A is supplied in the packed ! array AP as follows: ! ! UPLO = 'U' or 'u' The upper triangular part of A is ! supplied in AP. ! ! UPLO = 'L' or 'l' The lower triangular part of A is ! supplied in AP. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! X - 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. ! ! 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. On exit, the array ! AP is overwritten by the upper triangular part of the ! updated matrix. ! Before entry with UPLO = 'L' or 'l', the array AP must ! contain the lower triangular part of the 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. On exit, the array ! AP is overwritten by the lower triangular part of the ! updated matrix. ! ! ! 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, K, KK, KX ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. & .NOT.LSAME( UPLO, 'L' ) ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( INCX == 0 ) then INFO = 5 end if if ( INFO /= 0 ) then CALL XERBLA( 'DSPR ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. ( ALPHA == ZERO ) ) then return end if ! ! Set the start point in X if the increment is not unity. ! if ( INCX.LE.0 ) then KX = 1 - ( N - 1 )*INCX else if ( INCX /= 1 ) then KX = 1 end if ! ! Start the operations. In this version the elements of the array AP ! are accessed sequentially with one pass through AP. ! KK = 1 if ( LSAME( UPLO, 'U' ) ) then ! ! Form A when upper triangle is stored in AP. ! if ( INCX == 1 ) then DO 20, J = 1, N if ( X( J ) /= ZERO ) then TEMP = ALPHA*X( J ) K = KK DO 10, I = 1, J AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 10 CONTINUE end if KK = KK + J 20 CONTINUE else JX = KX DO 40, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) IX = KX DO 30, K = KK, KK + J - 1 AP( K ) = AP( K ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE end if JX = JX + INCX KK = KK + J 40 CONTINUE end if else ! ! Form A when lower triangle is stored in AP. ! if ( INCX == 1 ) then DO 60, J = 1, N if ( X( J ) /= ZERO ) then TEMP = ALPHA*X( J ) K = KK DO 50, I = J, N AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 50 CONTINUE end if KK = KK + N - J + 1 60 CONTINUE else JX = KX DO 80, J = 1, N if ( X( JX ) /= ZERO ) then TEMP = ALPHA*X( JX ) IX = JX DO 70, K = KK, KK + N - J AP( K ) = AP( K ) + X( IX )*TEMP IX = IX + INCX 70 CONTINUE end if JX = JX + INCX KK = KK + N - J + 1 80 CONTINUE end if end if return end subroutine dsyr2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) ! !******************************************************************************* ! !! DSYR2 ! DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, N character UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) ! .. ! ! Purpose ! ======= ! ! DSYR2 performs the symmetric rank 2 operation ! ! A := alpha*x*y' + alpha*y*x' + A, ! ! where alpha is a scalar, x and y are n element vectors and A is an n ! by n symmetric matrix. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array A is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of A ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of A ! is to be referenced. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! X - 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. ! ! 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. ! Unchanged on exit. ! ! INCY - INTEGER. ! On entry, INCY specifies the increment for the elements of ! Y. INCY must not be zero. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array A must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of A is not referenced. On exit, the ! upper triangular part of the array A is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array A must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of A is not referenced. On exit, the ! lower triangular part of the array A is overwritten by the ! lower triangular part of the updated matrix. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, n ). ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. & .NOT.LSAME( UPLO, 'L' ) ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( INCX == 0 ) then INFO = 5 else if ( INCY == 0 ) then INFO = 7 else if ( LDA < MAX( 1, N ) ) then INFO = 9 end if if ( INFO /= 0 ) then CALL XERBLA( 'DSYR2 ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. ( ALPHA == ZERO ) ) then return end if ! ! Set up the start points in X and Y if the increments are not both ! unity. ! if ( ( INCX /= 1 ) .or. ( INCY /= 1 ) ) then if ( INCX > 0 ) then KX = 1 else KX = 1 - ( N - 1 )*INCX end if if ( INCY > 0 ) then KY = 1 else KY = 1 - ( N - 1 )*INCY end if JX = KX JY = KY end if ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through the triangular part ! of A. ! if ( LSAME( UPLO, 'U' ) ) then ! ! Form A when A is stored in the upper triangle. ! if ( ( INCX == 1 ).AND.( INCY == 1 ) ) then DO 20, J = 1, N if ( ( X( J ) /= ZERO ) .or. ( Y( J ) /= ZERO ) ) then TEMP1 = ALPHA*Y( J ) TEMP2 = ALPHA*X( J ) DO 10, I = 1, J A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 10 CONTINUE end if 20 CONTINUE else DO 40, J = 1, N if ( ( X( JX ) /= ZERO ) .or. ( Y( JY ) /= ZERO ) ) then TEMP1 = ALPHA*Y( JY ) TEMP2 = ALPHA*X( JX ) IX = KX IY = KY DO 30, I = 1, J A( I, J ) = A( I, J ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE end if JX = JX + INCX JY = JY + INCY 40 CONTINUE end if else ! ! Form A when A is stored in the lower triangle. ! if ( ( INCX == 1 ).AND.( INCY == 1 ) ) then DO 60, J = 1, N if ( ( X( J ) /= ZERO ) .or. ( Y( J ) /= ZERO ) ) then TEMP1 = ALPHA*Y( J ) TEMP2 = ALPHA*X( J ) DO 50, I = J, N A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2 50 CONTINUE end if 60 CONTINUE else DO 80, J = 1, N if ( ( X( JX ) /= ZERO ) .or. ( Y( JY ) /= ZERO ) ) then TEMP1 = ALPHA*Y( JY ) TEMP2 = ALPHA*X( JX ) IX = JX IY = JY DO 70, I = J, N A( I, J ) = A( I, J ) + X( IX )*TEMP1 & + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 70 CONTINUE end if JX = JX + INCX JY = JY + INCY 80 CONTINUE end if end if return end subroutine dspr2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) ! !******************************************************************************* ! !! DSPR2 ! DOUBLE PRECISION ALPHA INTEGER INCX, INCY, N character UPLO ! .. Array Arguments .. DOUBLE PRECISION AP( * ), X( * ), Y( * ) ! .. ! ! Purpose ! ======= ! ! DSPR2 performs the symmetric rank 2 operation ! ! A := alpha*x*y' + alpha*y*x' + A, ! ! where alpha is a scalar, x and y are n element vectors and A is an ! n by n symmetric matrix, supplied in packed form. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the matrix A is supplied in the packed ! array AP as follows: ! ! UPLO = 'U' or 'u' The upper triangular part of A is ! supplied in AP. ! ! UPLO = 'L' or 'l' The lower triangular part of A is ! supplied in AP. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! X - 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. ! ! 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. ! Unchanged on exit. ! ! INCY - INTEGER. ! On entry, INCY specifies the increment for the elements of ! Y. INCY must not be zero. ! Unchanged on exit. ! ! AP - 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. On exit, the array ! AP is overwritten by the upper triangular part of the ! updated matrix. ! Before entry with UPLO = 'L' or 'l', the array AP must ! contain the lower triangular part of the 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. On exit, the array ! AP is overwritten by the lower triangular part of the ! updated matrix. ! ! ! 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 TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. & .NOT.LSAME( UPLO, 'L' ) ) then INFO = 1 else if ( N < 0 ) then INFO = 2 else if ( INCX == 0 ) then INFO = 5 else if ( INCY == 0 ) then INFO = 7 end if if ( INFO /= 0 ) then CALL XERBLA( 'DSPR2 ', INFO ) return end if ! ! Quick return if possible. ! if ( ( N == 0 ) .or. ( ALPHA == ZERO ) ) then return end if ! ! Set up the start points in X and Y if the increments are not both ! unity. ! if ( ( INCX /= 1 ) .or. ( INCY /= 1 ) ) then if ( INCX > 0 ) then KX = 1 else KX = 1 - ( N - 1 )*INCX end if if ( INCY > 0 ) then KY = 1 else KY = 1 - ( N - 1 )*INCY end if JX = KX JY = KY end if ! ! Start the operations. In this version the elements of the array AP ! are accessed sequentially with one pass through AP. ! KK = 1 if ( LSAME( UPLO, 'U' ) ) then ! ! Form A when upper triangle is stored in AP. ! if ( ( INCX == 1 ).AND.( INCY == 1 ) ) then DO 20, J = 1, N if ( ( X( J ) /= ZERO ) .or. ( Y( J ) /= ZERO ) ) then TEMP1 = ALPHA*Y( J ) TEMP2 = ALPHA*X( J ) K = KK DO 10, I = 1, J AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 10 CONTINUE end if KK = KK + J 20 CONTINUE else DO 40, J = 1, N if ( ( X( JX ) /= ZERO ) .or. ( Y( JY ) /= ZERO ) ) then TEMP1 = ALPHA*Y( JY ) TEMP2 = ALPHA*X( JX ) IX = KX IY = KY DO 30, K = KK, KK + J - 1 AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE end if JX = JX + INCX JY = JY + INCY KK = KK + J 40 CONTINUE end if else ! ! Form A when lower triangle is stored in AP. ! if ( ( INCX == 1 ).AND.( INCY == 1 ) ) then DO 60, J = 1, N if ( ( X( J ) /= ZERO ) .or. ( Y( J ) /= ZERO ) ) then TEMP1 = ALPHA*Y( J ) TEMP2 = ALPHA*X( J ) K = KK DO 50, I = J, N AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 50 CONTINUE end if KK = KK + N - J + 1 60 CONTINUE else DO 80, J = 1, N if ( ( X( JX ) /= ZERO ) .or. ( Y( JY ) /= ZERO ) ) then TEMP1 = ALPHA*Y( JY ) TEMP2 = ALPHA*X( JX ) IX = JX IY = JY DO 70, K = KK, KK + N - J AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 70 CONTINUE end if JX = JX + INCX JY = JY + INCY KK = KK + N - J + 1 80 CONTINUE end if end if return end