*> \brief \b DASUM
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DASUM takes the sum of the absolute values.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup asum
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 3/93 to return if incx .le. 0.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DTEMP
INTEGER I,M,MP1,NINCX
* ..
* .. Intrinsic Functions ..
INTRINSIC DABS,MOD
* ..
DASUM = 0.0d0
DTEMP = 0.0d0
IF (N.LE.0 .OR. INCX.LE.0) RETURN
IF (INCX.EQ.1) THEN
* code for increment equal to 1
*
*
* clean-up loop
*
M = MOD(N,6)
IF (M.NE.0) THEN
DO I = 1,M
DTEMP = DTEMP + DABS(DX(I))
END DO
IF (N.LT.6) THEN
DASUM = DTEMP
RETURN
END IF
END IF
MP1 = M + 1
DO I = MP1,N,6
DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I+1)) +
$ DABS(DX(I+2)) + DABS(DX(I+3)) +
$ DABS(DX(I+4)) + DABS(DX(I+5))
END DO
ELSE
*
* code for increment not equal to 1
*
NINCX = N*INCX
DO I = 1,NINCX,INCX
DTEMP = DTEMP + DABS(DX(I))
END DO
END IF
DASUM = DTEMP
RETURN
*
* End of DASUM
*
END
*> \brief \b DAXPY
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION DA
* INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*),DY(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DAXPY constant times a vector plus a vector.
*> uses unrolled loops for increments equal to one.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in] DA
*> \verbatim
*> DA is DOUBLE PRECISION
*> On entry, DA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*>
*> \param[in,out] DY
*> \verbatim
*> DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> storage spacing between elements of DY
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup axpy
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION DA
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*),DY(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
INTEGER I,IX,IY,M,MP1
* ..
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
IF (N.LE.0) RETURN
IF (DA.EQ.0.0d0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* code for both increments equal to 1
*
*
* clean-up loop
*
M = MOD(N,4)
IF (M.NE.0) THEN
DO I = 1,M
DY(I) = DY(I) + DA*DX(I)
END DO
END IF
IF (N.LT.4) RETURN
MP1 = M + 1
DO I = MP1,N,4
DY(I) = DY(I) + DA*DX(I)
DY(I+1) = DY(I+1) + DA*DX(I+1)
DY(I+2) = DY(I+2) + DA*DX(I+2)
DY(I+3) = DY(I+3) + DA*DX(I+3)
END DO
ELSE
*
* code for unequal increments or equal increments
* not equal to 1
*
IX = 1
IY = 1
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO I = 1,N
DY(IY) = DY(IY) + DA*DX(IX)
IX = IX + INCX
IY = IY + INCY
END DO
END IF
RETURN
*
* End of DAXPY
*
END
*> \brief \b DCOPY
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
*
* .. Scalar Arguments ..
* INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*),DY(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DCOPY copies a vector, x, to a vector, y.
*> uses unrolled loops for increments equal to 1.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*>
*> \param[out] DY
*> \verbatim
*> DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> storage spacing between elements of DY
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup copy
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*),DY(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
INTEGER I,IX,IY,M,MP1
* ..
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
IF (N.LE.0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* code for both increments equal to 1
*
*
* clean-up loop
*
M = MOD(N,7)
IF (M.NE.0) THEN
DO I = 1,M
DY(I) = DX(I)
END DO
IF (N.LT.7) RETURN
END IF
MP1 = M + 1
DO I = MP1,N,7
DY(I) = DX(I)
DY(I+1) = DX(I+1)
DY(I+2) = DX(I+2)
DY(I+3) = DX(I+3)
DY(I+4) = DX(I+4)
DY(I+5) = DX(I+5)
DY(I+6) = DX(I+6)
END DO
ELSE
*
* code for unequal increments or equal increments
* not equal to 1
*
IX = 1
IY = 1
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO I = 1,N
DY(IY) = DX(IX)
IX = IX + INCX
IY = IY + INCY
END DO
END IF
RETURN
*
* End of DCOPY
*
END
*> \brief \b DDOT
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
*
* .. Scalar Arguments ..
* INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*),DY(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DDOT forms the dot product of two vectors.
*> uses unrolled loops for increments equal to one.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*>
*> \param[in] DY
*> \verbatim
*> DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> storage spacing between elements of DY
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup dot
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*),DY(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DTEMP
INTEGER I,IX,IY,M,MP1
* ..
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
DDOT = 0.0d0
DTEMP = 0.0d0
IF (N.LE.0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* code for both increments equal to 1
*
*
* clean-up loop
*
M = MOD(N,5)
IF (M.NE.0) THEN
DO I = 1,M
DTEMP = DTEMP + DX(I)*DY(I)
END DO
IF (N.LT.5) THEN
DDOT=DTEMP
RETURN
END IF
END IF
MP1 = M + 1
DO I = MP1,N,5
DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) +
$ DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
END DO
ELSE
*
* code for unequal increments or equal increments
* not equal to 1
*
IX = 1
IY = 1
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO I = 1,N
DTEMP = DTEMP + DX(IX)*DY(IY)
IX = IX + INCX
IY = IY + INCY
END DO
END IF
DDOT = DTEMP
RETURN
*
* End of DDOT
*
END
*> \brief \b DGBMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER INCX,INCY,KL,KU,LDA,M,N
* CHARACTER TRANS
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DGBMV performs one of the matrix-vector operations
*>
*> y := alpha*A*x + beta*y, or y := alpha*A**T*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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x + beta*y.
*>
*> TRANS = 'C' or 'c' y := alpha*A**T*x + beta*y.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> On entry, M specifies the number of rows of the matrix A.
*> M must be at least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of columns of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] KL
*> \verbatim
*> KL is INTEGER
*> On entry, KL specifies the number of sub-diagonals of the
*> matrix A. KL must satisfy 0 .le. KL.
*> \endverbatim
*>
*> \param[in] KU
*> \verbatim
*> KU is INTEGER
*> On entry, KU specifies the number of super-diagonals of the
*> matrix A. KU must satisfy 0 .le. KU.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta. When BETA is
*> supplied as zero then Y need not be set on input.
*> \endverbatim
*>
*> \param[in,out] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, 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.
*> If either m or n is zero, then Y not referenced and the function
*> performs a quick return.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup gbmv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,
+ BETA,Y,INCY)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER INCX,INCY,KL,KU,LDA,M,N
CHARACTER TRANS
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. 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,LENX,LENY
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX,MIN
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 1
ELSE IF (M.LT.0) THEN
INFO = 2
ELSE IF (N.LT.0) THEN
INFO = 3
ELSE IF (KL.LT.0) THEN
INFO = 4
ELSE IF (KU.LT.0) THEN
INFO = 5
ELSE IF (LDA.LT. (KL+KU+1)) THEN
INFO = 8
ELSE IF (INCX.EQ.0) THEN
INFO = 10
ELSE IF (INCY.EQ.0) THEN
INFO = 13
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DGBMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
+ ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
*
* Set LENX and LENY, the lengths of the vectors x and y, and set
* up the start points in X and Y.
*
IF (LSAME(TRANS,'N')) THEN
LENX = N
LENY = M
ELSE
LENX = M
LENY = N
END IF
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (LENX-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (LENY-1)*INCY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the band part of A.
*
* First form y := beta*y.
*
IF (BETA.NE.ONE) THEN
IF (INCY.EQ.1) THEN
IF (BETA.EQ.ZERO) THEN
DO 10 I = 1,LENY
Y(I) = ZERO
10 CONTINUE
ELSE
DO 20 I = 1,LENY
Y(I) = BETA*Y(I)
20 CONTINUE
END IF
ELSE
IY = KY
IF (BETA.EQ.ZERO) THEN
DO 30 I = 1,LENY
Y(IY) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40 I = 1,LENY
Y(IY) = BETA*Y(IY)
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF (ALPHA.EQ.ZERO) RETURN
KUP1 = KU + 1
IF (LSAME(TRANS,'N')) THEN
*
* Form y := alpha*A*x + y.
*
JX = KX
IF (INCY.EQ.1) THEN
DO 60 J = 1,N
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
JX = JX + INCX
60 CONTINUE
ELSE
DO 80 J = 1,N
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
JX = JX + INCX
IF (J.GT.KU) KY = KY + INCY
80 CONTINUE
END IF
ELSE
*
* Form y := alpha*A**T*x + y.
*
JY = KY
IF (INCX.EQ.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.GT.KU) KX = KX + INCX
120 CONTINUE
END IF
END IF
*
RETURN
*
* End of DGBMV
*
END
*> \brief \b DGEMM
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER K,LDA,LDB,LDC,M,N
* CHARACTER TRANSA,TRANSB
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DGEMM performs one of the matrix-matrix operations
*>
*> C := alpha*op( A )*op( B ) + beta*C,
*>
*> where op( X ) is one of
*>
*> op( X ) = X or op( X ) = X**T,
*>
*> alpha and beta are scalars, and A, B and C are matrices, with op( A )
*> an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] TRANSA
*> \verbatim
*> TRANSA is CHARACTER*1
*> On entry, TRANSA specifies the form of op( A ) to be used in
*> the matrix multiplication as follows:
*>
*> TRANSA = 'N' or 'n', op( A ) = A.
*>
*> TRANSA = 'T' or 't', op( A ) = A**T.
*>
*> TRANSA = 'C' or 'c', op( A ) = A**T.
*> \endverbatim
*>
*> \param[in] TRANSB
*> \verbatim
*> TRANSB is CHARACTER*1
*> On entry, TRANSB specifies the form of op( B ) to be used in
*> the matrix multiplication as follows:
*>
*> TRANSB = 'N' or 'n', op( B ) = B.
*>
*> TRANSB = 'T' or 't', op( B ) = B**T.
*>
*> TRANSB = 'C' or 'c', op( B ) = B**T.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> On entry, M specifies the number of rows of the matrix
*> op( A ) and of the matrix C. M must be at least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of columns of the matrix
*> op( B ) and the number of columns of the matrix C. N must be
*> at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*> K is INTEGER
*> On entry, K specifies the number of columns of the matrix
*> op( A ) and the number of rows of the matrix op( B ). K must
*> be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
*> k when TRANSA = 'N' or 'n', and is m otherwise.
*> Before entry with TRANSA = 'N' or 'n', the leading m by k
*> part of the array A must contain the matrix A, otherwise
*> the leading k by m part of the array A must contain the
*> matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. When TRANSA = 'N' or 'n' then
*> LDA must be at least max( 1, m ), otherwise LDA must be at
*> least max( 1, k ).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is
*> n when TRANSB = 'N' or 'n', and is k otherwise.
*> Before entry with TRANSB = 'N' or 'n', the leading k by n
*> part of the array B must contain the matrix B, otherwise
*> the leading n by k part of the array B must contain the
*> matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> On entry, LDB specifies the first dimension of B as declared
*> in the calling (sub) program. When TRANSB = 'N' or 'n' then
*> LDB must be at least max( 1, k ), otherwise LDB must be at
*> least max( 1, n ).
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta. When BETA is
*> supplied as zero then C need not be set on input.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*> C is DOUBLE PRECISION array, dimension ( LDC, N )
*> Before entry, the leading m by n part of the array C must
*> contain the matrix C, except when beta is zero, in which
*> case C need not be set on entry.
*> On exit, the array C is overwritten by the m by n matrix
*> ( alpha*op( A )*op( B ) + beta*C ).
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*> LDC is INTEGER
*> On entry, LDC specifies the first dimension of C as declared
*> in the calling (sub) program. LDC must be at least
*> max( 1, m ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup gemm
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 3 Blas routine.
*>
*> -- Written on 8-February-1989.
*> Jack Dongarra, Argonne National Laboratory.
*> Iain Duff, AERE Harwell.
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
*> Sven Hammarling, Numerical Algorithms Group Ltd.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,
+ BETA,C,LDC)
*
* -- Reference BLAS level3 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER K,LDA,LDB,LDC,M,N
CHARACTER TRANSA,TRANSB
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I,INFO,J,L,NROWA,NROWB
LOGICAL NOTA,NOTB
* ..
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
*
* Set NOTA and NOTB as true if A and B respectively are not
* transposed and set NROWA and NROWB as the number of rows of A
* and B respectively.
*
NOTA = LSAME(TRANSA,'N')
NOTB = LSAME(TRANSB,'N')
IF (NOTA) THEN
NROWA = M
ELSE
NROWA = K
END IF
IF (NOTB) THEN
NROWB = K
ELSE
NROWB = N
END IF
*
* Test the input parameters.
*
INFO = 0
IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND.
+ (.NOT.LSAME(TRANSA,'T'))) THEN
INFO = 1
ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND.
+ (.NOT.LSAME(TRANSB,'T'))) THEN
INFO = 2
ELSE IF (M.LT.0) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (K.LT.0) THEN
INFO = 5
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
INFO = 8
ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
INFO = 10
ELSE IF (LDC.LT.MAX(1,M)) THEN
INFO = 13
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DGEMM ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
+ (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
*
* And if alpha.eq.zero.
*
IF (ALPHA.EQ.ZERO) THEN
IF (BETA.EQ.ZERO) THEN
DO 20 J = 1,N
DO 10 I = 1,M
C(I,J) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1,N
DO 30 I = 1,M
C(I,J) = BETA*C(I,J)
30 CONTINUE
40 CONTINUE
END IF
RETURN
END IF
*
* Start the operations.
*
IF (NOTB) THEN
IF (NOTA) THEN
*
* Form C := alpha*A*B + beta*C.
*
DO 90 J = 1,N
IF (BETA.EQ.ZERO) THEN
DO 50 I = 1,M
C(I,J) = ZERO
50 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 60 I = 1,M
C(I,J) = BETA*C(I,J)
60 CONTINUE
END IF
DO 80 L = 1,K
TEMP = ALPHA*B(L,J)
DO 70 I = 1,M
C(I,J) = C(I,J) + TEMP*A(I,L)
70 CONTINUE
80 CONTINUE
90 CONTINUE
ELSE
*
* Form C := alpha*A**T*B + beta*C
*
DO 120 J = 1,N
DO 110 I = 1,M
TEMP = ZERO
DO 100 L = 1,K
TEMP = TEMP + A(L,I)*B(L,J)
100 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP
ELSE
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
END IF
110 CONTINUE
120 CONTINUE
END IF
ELSE
IF (NOTA) THEN
*
* Form C := alpha*A*B**T + beta*C
*
DO 170 J = 1,N
IF (BETA.EQ.ZERO) THEN
DO 130 I = 1,M
C(I,J) = ZERO
130 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 140 I = 1,M
C(I,J) = BETA*C(I,J)
140 CONTINUE
END IF
DO 160 L = 1,K
TEMP = ALPHA*B(J,L)
DO 150 I = 1,M
C(I,J) = C(I,J) + TEMP*A(I,L)
150 CONTINUE
160 CONTINUE
170 CONTINUE
ELSE
*
* Form C := alpha*A**T*B**T + beta*C
*
DO 200 J = 1,N
DO 190 I = 1,M
TEMP = ZERO
DO 180 L = 1,K
TEMP = TEMP + A(L,I)*B(J,L)
180 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP
ELSE
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
END IF
190 CONTINUE
200 CONTINUE
END IF
END IF
*
RETURN
*
* End of DGEMM
*
END
*> \brief \b DGEMMTR
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DGEMMTR(UPLO,TRANSA,TRANSB,N,K,ALPHA,A,LDA,B,LDB,BETA,
* C,LDC)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER K,LDA,LDB,LDC,N
* CHARACTER TRANSA,TRANSB, UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DGEMMTR performs one of the matrix-matrix operations
*>
*> C := alpha*op( A )*op( B ) + beta*C,
*>
*> where op( X ) is one of
*>
*> op( X ) = X or op( X ) = X**T,
*>
*> alpha and beta are scalars, and A, B and C are matrices, with op( A )
*> an n by k matrix, op( B ) a k by n matrix and C an n by n matrix.
*> Thereby, the routine only accesses and updates the upper or lower
*> triangular part of the result matrix C. This behaviour can be used if
*> the resulting matrix C is known to be symmetric.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> On entry, UPLO specifies whether the lower or the upper
*> triangular part of C is access and updated.
*>
*> UPLO = 'L' or 'l', the lower triangular part of C is used.
*>
*> UPLO = 'U' or 'u', the upper triangular part of C is used.
*> \endverbatim
*
*> \param[in] TRANSA
*> \verbatim
*> TRANSA is CHARACTER*1
*> On entry, TRANSA specifies the form of op( A ) to be used in
*> the matrix multiplication as follows:
*>
*> TRANSA = 'N' or 'n', op( A ) = A.
*>
*> TRANSA = 'T' or 't', op( A ) = A**T.
*>
*> TRANSA = 'C' or 'c', op( A ) = A**T.
*> \endverbatim
*>
*> \param[in] TRANSB
*> \verbatim
*> TRANSB is CHARACTER*1
*> On entry, TRANSB specifies the form of op( B ) to be used in
*> the matrix multiplication as follows:
*>
*> TRANSB = 'N' or 'n', op( B ) = B.
*>
*> TRANSB = 'T' or 't', op( B ) = B**T.
*>
*> TRANSB = 'C' or 'c', op( B ) = B**T.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of rows and columns of
*> the matrix C, the number of columns of op(B) and the number
*> of rows of op(A). N must be at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*> K is INTEGER
*> On entry, K specifies the number of columns of the matrix
*> op( A ) and the number of rows of the matrix op( B ). K must
*> be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
*> k when TRANSA = 'N' or 'n', and is n otherwise.
*> Before entry with TRANSA = 'N' or 'n', the leading n by k
*> part of the array A must contain the matrix A, otherwise
*> the leading k by m part of the array A must contain the
*> matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. When TRANSA = 'N' or 'n' then
*> LDA must be at least max( 1, n ), otherwise LDA must be at
*> least max( 1, k ).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is
*> n when TRANSB = 'N' or 'n', and is k otherwise.
*> Before entry with TRANSB = 'N' or 'n', the leading k by n
*> part of the array B must contain the matrix B, otherwise
*> the leading n by k part of the array B must contain the
*> matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> On entry, LDB specifies the first dimension of B as declared
*> in the calling (sub) program. When TRANSB = 'N' or 'n' then
*> LDB must be at least max( 1, k ), otherwise LDB must be at
*> least max( 1, n ).
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta. When BETA is
*> supplied as zero then C need not be set on input.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*> C is DOUBLE PRECISION array, dimension ( LDC, N )
*> Before entry, the leading n by n part of the array C must
*> contain the matrix C, except when beta is zero, in which
*> case C need not be set on entry.
*> On exit, the upper or lower triangular part of the matrix
*> C is overwritten by the n by n matrix
*> ( alpha*op( A )*op( B ) + beta*C ).
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*> LDC is INTEGER
*> On entry, LDC specifies the first dimension of C as declared
*> in the calling (sub) program. LDC must be at least
*> max( 1, n ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Martin Koehler
*
*> \ingroup gemmtr
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 3 Blas routine.
*>
*> -- Written on 19-July-2023.
*> Martin Koehler, MPI Magdeburg
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DGEMMTR(UPLO,TRANSA,TRANSB,N,K,ALPHA,A,LDA,B,LDB,
+ BETA,C,LDC)
IMPLICIT NONE
*
* -- Reference BLAS level3 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER K,LDA,LDB,LDC,N
CHARACTER TRANSA,TRANSB,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I,INFO,J,L,NROWA,NROWB, ISTART, ISTOP
LOGICAL NOTA,NOTB, UPPER
* ..
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
*
* Set NOTA and NOTB as true if A and B respectively are not
* transposed and set NROWA and NROWB as the number of rows of A
* and B respectively.
*
NOTA = LSAME(TRANSA,'N')
NOTB = LSAME(TRANSB,'N')
IF (NOTA) THEN
NROWA = N
ELSE
NROWA = K
END IF
IF (NOTB) THEN
NROWB = K
ELSE
NROWB = N
END IF
UPPER = LSAME(UPLO, 'U')
*
* Test the input parameters.
*
INFO = 0
IF ((.NOT. UPPER) .AND. (.NOT. LSAME(UPLO, 'L'))) THEN
INFO = 1
ELSE IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND.
+ (.NOT.LSAME(TRANSA,'T'))) THEN
INFO = 2
ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND.
+ (.NOT.LSAME(TRANSB,'T'))) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (K.LT.0) THEN
INFO = 5
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
INFO = 8
ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
INFO = 10
ELSE IF (LDC.LT.MAX(1,N)) THEN
INFO = 13
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DGEMMTR',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (N.EQ.0) RETURN
*
* And if alpha.eq.zero.
*
IF (ALPHA.EQ.ZERO) THEN
IF (BETA.EQ.ZERO) THEN
DO 20 J = 1,N
IF (UPPER) THEN
ISTART = 1
ISTOP = J
ELSE
ISTART = J
ISTOP = N
END IF
DO 10 I = ISTART, ISTOP
C(I,J) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1,N
IF (UPPER) THEN
ISTART = 1
ISTOP = J
ELSE
ISTART = J
ISTOP = N
END IF
DO 30 I = ISTART, ISTOP
C(I,J) = BETA*C(I,J)
30 CONTINUE
40 CONTINUE
END IF
RETURN
END IF
*
* Start the operations.
*
IF (NOTB) THEN
IF (NOTA) THEN
*
* Form C := alpha*A*B + beta*C.
*
DO 90 J = 1,N
IF (UPPER) THEN
ISTART = 1
ISTOP = J
ELSE
ISTART = J
ISTOP = N
END IF
IF (BETA.EQ.ZERO) THEN
DO 50 I = ISTART, ISTOP
C(I,J) = ZERO
50 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 60 I = ISTART, ISTOP
C(I,J) = BETA*C(I,J)
60 CONTINUE
END IF
DO 80 L = 1,K
TEMP = ALPHA*B(L,J)
DO 70 I = ISTART, ISTOP
C(I,J) = C(I,J) + TEMP*A(I,L)
70 CONTINUE
80 CONTINUE
90 CONTINUE
ELSE
*
* Form C := alpha*A**T*B + beta*C
*
DO 120 J = 1,N
IF (UPPER) THEN
ISTART = 1
ISTOP = J
ELSE
ISTART = J
ISTOP = N
END IF
DO 110 I = ISTART, ISTOP
TEMP = ZERO
DO 100 L = 1,K
TEMP = TEMP + A(L,I)*B(L,J)
100 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP
ELSE
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
END IF
110 CONTINUE
120 CONTINUE
END IF
ELSE
IF (NOTA) THEN
*
* Form C := alpha*A*B**T + beta*C
*
DO 170 J = 1,N
IF (UPPER) THEN
ISTART = 1
ISTOP = J
ELSE
ISTART = J
ISTOP = N
END IF
IF (BETA.EQ.ZERO) THEN
DO 130 I = ISTART,ISTOP
C(I,J) = ZERO
130 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 140 I = ISTART,ISTOP
C(I,J) = BETA*C(I,J)
140 CONTINUE
END IF
DO 160 L = 1,K
TEMP = ALPHA*B(J,L)
DO 150 I = ISTART,ISTOP
C(I,J) = C(I,J) + TEMP*A(I,L)
150 CONTINUE
160 CONTINUE
170 CONTINUE
ELSE
*
* Form C := alpha*A**T*B**T + beta*C
*
DO 200 J = 1,N
IF (UPPER) THEN
ISTART = 1
ISTOP = J
ELSE
ISTART = J
ISTOP = N
END IF
DO 190 I = ISTART, ISTOP
TEMP = ZERO
DO 180 L = 1,K
TEMP = TEMP + A(L,I)*B(J,L)
180 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP
ELSE
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
END IF
190 CONTINUE
200 CONTINUE
END IF
END IF
*
RETURN
*
* End of SGEMM
*
END
*> \brief \b DGEMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER INCX,INCY,LDA,M,N
* CHARACTER TRANS
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DGEMV performs one of the matrix-vector operations
*>
*> y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y,
*>
*> where alpha and beta are scalars, x and y are vectors and A is an
*> m by n matrix.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x + beta*y.
*>
*> TRANS = 'C' or 'c' y := alpha*A**T*x + beta*y.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> On entry, M specifies the number of rows of the matrix A.
*> M must be at least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of columns of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, N )
*> Before entry, the leading m by n part of the array A must
*> contain the matrix of coefficients.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta. When BETA is
*> supplied as zero then Y need not be set on input.
*> \endverbatim
*>
*> \param[in,out] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, 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.
*> If either m or n is zero, then Y not referenced and the function
*> performs a quick return.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup gemv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER INCX,INCY,LDA,M,N
CHARACTER TRANS
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 1
ELSE IF (M.LT.0) THEN
INFO = 2
ELSE IF (N.LT.0) THEN
INFO = 3
ELSE IF (LDA.LT.MAX(1,M)) THEN
INFO = 6
ELSE IF (INCX.EQ.0) THEN
INFO = 8
ELSE IF (INCY.EQ.0) THEN
INFO = 11
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DGEMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
+ ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
*
* Set LENX and LENY, the lengths of the vectors x and y, and set
* up the start points in X and Y.
*
IF (LSAME(TRANS,'N')) THEN
LENX = N
LENY = M
ELSE
LENX = M
LENY = N
END IF
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (LENX-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (LENY-1)*INCY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
* First form y := beta*y.
*
IF (BETA.NE.ONE) THEN
IF (INCY.EQ.1) THEN
IF (BETA.EQ.ZERO) THEN
DO 10 I = 1,LENY
Y(I) = ZERO
10 CONTINUE
ELSE
DO 20 I = 1,LENY
Y(I) = BETA*Y(I)
20 CONTINUE
END IF
ELSE
IY = KY
IF (BETA.EQ.ZERO) THEN
DO 30 I = 1,LENY
Y(IY) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40 I = 1,LENY
Y(IY) = BETA*Y(IY)
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF (ALPHA.EQ.ZERO) RETURN
IF (LSAME(TRANS,'N')) THEN
*
* Form y := alpha*A*x + y.
*
JX = KX
IF (INCY.EQ.1) THEN
DO 60 J = 1,N
TEMP = ALPHA*X(JX)
DO 50 I = 1,M
Y(I) = Y(I) + TEMP*A(I,J)
50 CONTINUE
JX = JX + INCX
60 CONTINUE
ELSE
DO 80 J = 1,N
TEMP = ALPHA*X(JX)
IY = KY
DO 70 I = 1,M
Y(IY) = Y(IY) + TEMP*A(I,J)
IY = IY + INCY
70 CONTINUE
JX = JX + INCX
80 CONTINUE
END IF
ELSE
*
* Form y := alpha*A**T*x + y.
*
JY = KY
IF (INCX.EQ.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 of DGEMV
*
END
*> \brief \b DGER
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA
* INTEGER INCX,INCY,LDA,M,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DGER performs the rank 1 operation
*>
*> A := alpha*x*y**T + 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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> On entry, M specifies the number of rows of the matrix A.
*> M must be at least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of columns of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( m - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the m
*> element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCY ) ).
*> Before entry, the incremented array Y must contain the n
*> element vector y.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup ger
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER INCX,INCY,LDA,M,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (M.LT.0) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (INCX.EQ.0) THEN
INFO = 5
ELSE IF (INCY.EQ.0) THEN
INFO = 7
ELSE IF (LDA.LT.MAX(1,M)) THEN
INFO = 9
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DGER ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((M.EQ.0) .OR. (N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF (INCY.GT.0) THEN
JY = 1
ELSE
JY = 1 - (N-1)*INCY
END IF
IF (INCX.EQ.1) THEN
DO 20 J = 1,N
IF (Y(JY).NE.ZERO) THEN
TEMP = ALPHA*Y(JY)
DO 10 I = 1,M
A(I,J) = A(I,J) + X(I)*TEMP
10 CONTINUE
END IF
JY = JY + INCY
20 CONTINUE
ELSE
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (M-1)*INCX
END IF
DO 40 J = 1,N
IF (Y(JY).NE.ZERO) THEN
TEMP = ALPHA*Y(JY)
IX = KX
DO 30 I = 1,M
A(I,J) = A(I,J) + X(IX)*TEMP
IX = IX + INCX
30 CONTINUE
END IF
JY = JY + INCY
40 CONTINUE
END IF
*
RETURN
*
* End of DGER
*
END
*> \brief \b DROT
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DROT(N,DX,INCX,DY,INCY,C,S)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION C,S
* INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*),DY(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DROT applies a plane rotation.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in,out] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*>
*> \param[in,out] DY
*> \verbatim
*> DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> storage spacing between elements of DY
*> \endverbatim
*>
*> \param[in] C
*> \verbatim
*> C is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] S
*> \verbatim
*> S is DOUBLE PRECISION
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup rot
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DROT(N,DX,INCX,DY,INCY,C,S)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION C,S
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*),DY(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DTEMP
INTEGER I,IX,IY
* ..
IF (N.LE.0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* code for both increments equal to 1
*
DO I = 1,N
DTEMP = C*DX(I) + S*DY(I)
DY(I) = C*DY(I) - S*DX(I)
DX(I) = DTEMP
END DO
ELSE
*
* code for unequal increments or equal increments not equal
* to 1
*
IX = 1
IY = 1
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO I = 1,N
DTEMP = C*DX(IX) + S*DY(IY)
DY(IY) = C*DY(IY) - S*DX(IX)
DX(IX) = DTEMP
IX = IX + INCX
IY = IY + INCY
END DO
END IF
RETURN
*
* End of DROT
*
END
*> \brief \b DROTM
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DROTM(N,DX,INCX,DY,INCY,DPARAM)
*
* .. Scalar Arguments ..
* INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DPARAM(5),DX(*),DY(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
*>
*> (DX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF DX ARE IN
*> (DY**T)
*>
*> DX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE
*> LX = (-INCX)*N, AND SIMILARLY FOR SY USING LY AND INCY.
*> WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
*>
*> DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0
*>
*> (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0)
*> H=( ) ( ) ( ) ( )
*> (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0).
*> SEE DROTMG FOR A DESCRIPTION OF DATA STORAGE IN DPARAM.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in,out] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*>
*> \param[in,out] DY
*> \verbatim
*> DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> storage spacing between elements of DY
*> \endverbatim
*>
*> \param[in] DPARAM
*> \verbatim
*> DPARAM is DOUBLE PRECISION array, dimension (5)
*> DPARAM(1)=DFLAG
*> DPARAM(2)=DH11
*> DPARAM(3)=DH21
*> DPARAM(4)=DH12
*> DPARAM(5)=DH22
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup rotm
*
* =====================================================================
SUBROUTINE DROTM(N,DX,INCX,DY,INCY,DPARAM)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DPARAM(5),DX(*),DY(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,TWO,W,Z,ZERO
INTEGER I,KX,KY,NSTEPS
* ..
* .. Data statements ..
DATA ZERO,TWO/0.D0,2.D0/
* ..
*
DFLAG = DPARAM(1)
IF (N.LE.0 .OR. (DFLAG+TWO.EQ.ZERO)) RETURN
IF (INCX.EQ.INCY.AND.INCX.GT.0) THEN
*
NSTEPS = N*INCX
IF (DFLAG.LT.ZERO) THEN
DH11 = DPARAM(2)
DH12 = DPARAM(4)
DH21 = DPARAM(3)
DH22 = DPARAM(5)
DO I = 1,NSTEPS,INCX
W = DX(I)
Z = DY(I)
DX(I) = W*DH11 + Z*DH12
DY(I) = W*DH21 + Z*DH22
END DO
ELSE IF (DFLAG.EQ.ZERO) THEN
DH12 = DPARAM(4)
DH21 = DPARAM(3)
DO I = 1,NSTEPS,INCX
W = DX(I)
Z = DY(I)
DX(I) = W + Z*DH12
DY(I) = W*DH21 + Z
END DO
ELSE
DH11 = DPARAM(2)
DH22 = DPARAM(5)
DO I = 1,NSTEPS,INCX
W = DX(I)
Z = DY(I)
DX(I) = W*DH11 + Z
DY(I) = -W + DH22*Z
END DO
END IF
ELSE
KX = 1
KY = 1
IF (INCX.LT.0) KX = 1 + (1-N)*INCX
IF (INCY.LT.0) KY = 1 + (1-N)*INCY
*
IF (DFLAG.LT.ZERO) THEN
DH11 = DPARAM(2)
DH12 = DPARAM(4)
DH21 = DPARAM(3)
DH22 = DPARAM(5)
DO I = 1,N
W = DX(KX)
Z = DY(KY)
DX(KX) = W*DH11 + Z*DH12
DY(KY) = W*DH21 + Z*DH22
KX = KX + INCX
KY = KY + INCY
END DO
ELSE IF (DFLAG.EQ.ZERO) THEN
DH12 = DPARAM(4)
DH21 = DPARAM(3)
DO I = 1,N
W = DX(KX)
Z = DY(KY)
DX(KX) = W + Z*DH12
DY(KY) = W*DH21 + Z
KX = KX + INCX
KY = KY + INCY
END DO
ELSE
DH11 = DPARAM(2)
DH22 = DPARAM(5)
DO I = 1,N
W = DX(KX)
Z = DY(KY)
DX(KX) = W*DH11 + Z
DY(KY) = -W + DH22*Z
KX = KX + INCX
KY = KY + INCY
END DO
END IF
END IF
RETURN
*
* End of DROTM
*
END
*> \brief \b DROTMG
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION DD1,DD2,DX1,DY1
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DPARAM(5)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
*> THE SECOND COMPONENT OF THE 2-VECTOR (DSQRT(DD1)*DX1,DSQRT(DD2)*> DY2)**T.
*> WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
*>
*> DFLAG=-1.D0 DFLAG=0.D0 DFLAG=1.D0 DFLAG=-2.D0
*>
*> (DH11 DH12) (1.D0 DH12) (DH11 1.D0) (1.D0 0.D0)
*> H=( ) ( ) ( ) ( )
*> (DH21 DH22), (DH21 1.D0), (-1.D0 DH22), (0.D0 1.D0).
*> LOCATIONS 2-4 OF DPARAM CONTAIN DH11, DH21, DH12, AND DH22
*> RESPECTIVELY. (VALUES OF 1.D0, -1.D0, OR 0.D0 IMPLIED BY THE
*> VALUE OF DPARAM(1) ARE NOT STORED IN DPARAM.)
*>
*> THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
*> INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
*> OF DD1 AND DD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
*>
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in,out] DD1
*> \verbatim
*> DD1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DD2
*> \verbatim
*> DD2 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in,out] DX1
*> \verbatim
*> DX1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[in] DY1
*> \verbatim
*> DY1 is DOUBLE PRECISION
*> \endverbatim
*>
*> \param[out] DPARAM
*> \verbatim
*> DPARAM is DOUBLE PRECISION array, dimension (5)
*> DPARAM(1)=DFLAG
*> DPARAM(2)=DH11
*> DPARAM(3)=DH21
*> DPARAM(4)=DH12
*> DPARAM(5)=DH22
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup rotmg
*
* =====================================================================
SUBROUTINE DROTMG(DD1,DD2,DX1,DY1,DPARAM)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION DD1,DD2,DX1,DY1
* ..
* .. Array Arguments ..
DOUBLE PRECISION DPARAM(5)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,DP1,DP2,DQ1,DQ2,DTEMP,
$ DU,GAM,GAMSQ,ONE,RGAMSQ,TWO,ZERO
* ..
* .. Intrinsic Functions ..
INTRINSIC DABS
* ..
* .. Data statements ..
*
DATA ZERO,ONE,TWO/0.D0,1.D0,2.D0/
DATA GAM,GAMSQ,RGAMSQ/4096.D0,16777216.D0,5.9604645D-8/
* ..
IF (DD1.LT.ZERO) THEN
* GO ZERO-H-D-AND-DX1..
DFLAG = -ONE
DH11 = ZERO
DH12 = ZERO
DH21 = ZERO
DH22 = ZERO
*
DD1 = ZERO
DD2 = ZERO
DX1 = ZERO
ELSE
* CASE-DD1-NONNEGATIVE
DP2 = DD2*DY1
IF (DP2.EQ.ZERO) THEN
DFLAG = -TWO
DPARAM(1) = DFLAG
RETURN
END IF
* REGULAR-CASE..
DP1 = DD1*DX1
DQ2 = DP2*DY1
DQ1 = DP1*DX1
*
IF (DABS(DQ1).GT.DABS(DQ2)) THEN
DH21 = -DY1/DX1
DH12 = DP2/DP1
*
DU = ONE - DH12*DH21
*
IF (DU.GT.ZERO) THEN
DFLAG = ZERO
DD1 = DD1/DU
DD2 = DD2/DU
DX1 = DX1*DU
ELSE
* This code path if here for safety. We do not expect this
* condition to ever hold except in edge cases with rounding
* errors. See DOI: 10.1145/355841.355847
DFLAG = -ONE
DH11 = ZERO
DH12 = ZERO
DH21 = ZERO
DH22 = ZERO
*
DD1 = ZERO
DD2 = ZERO
DX1 = ZERO
END IF
ELSE
IF (DQ2.LT.ZERO) THEN
* GO ZERO-H-D-AND-DX1..
DFLAG = -ONE
DH11 = ZERO
DH12 = ZERO
DH21 = ZERO
DH22 = ZERO
*
DD1 = ZERO
DD2 = ZERO
DX1 = ZERO
ELSE
DFLAG = ONE
DH11 = DP1/DP2
DH22 = DX1/DY1
DU = ONE + DH11*DH22
DTEMP = DD2/DU
DD2 = DD1/DU
DD1 = DTEMP
DX1 = DY1*DU
END IF
END IF
* PROCEDURE..SCALE-CHECK
IF (DD1.NE.ZERO) THEN
DO WHILE ((DD1.LE.RGAMSQ) .OR. (DD1.GE.GAMSQ))
IF (DFLAG.EQ.ZERO) THEN
DH11 = ONE
DH22 = ONE
DFLAG = -ONE
ELSE
DH21 = -ONE
DH12 = ONE
DFLAG = -ONE
END IF
IF (DD1.LE.RGAMSQ) THEN
DD1 = DD1*GAM**2
DX1 = DX1/GAM
DH11 = DH11/GAM
DH12 = DH12/GAM
ELSE
DD1 = DD1/GAM**2
DX1 = DX1*GAM
DH11 = DH11*GAM
DH12 = DH12*GAM
END IF
ENDDO
END IF
IF (DD2.NE.ZERO) THEN
DO WHILE ( (DABS(DD2).LE.RGAMSQ) .OR. (DABS(DD2).GE.GAMSQ) )
IF (DFLAG.EQ.ZERO) THEN
DH11 = ONE
DH22 = ONE
DFLAG = -ONE
ELSE
DH21 = -ONE
DH12 = ONE
DFLAG = -ONE
END IF
IF (DABS(DD2).LE.RGAMSQ) THEN
DD2 = DD2*GAM**2
DH21 = DH21/GAM
DH22 = DH22/GAM
ELSE
DD2 = DD2/GAM**2
DH21 = DH21*GAM
DH22 = DH22*GAM
END IF
END DO
END IF
END IF
IF (DFLAG.LT.ZERO) THEN
DPARAM(2) = DH11
DPARAM(3) = DH21
DPARAM(4) = DH12
DPARAM(5) = DH22
ELSE IF (DFLAG.EQ.ZERO) THEN
DPARAM(3) = DH21
DPARAM(4) = DH12
ELSE
DPARAM(2) = DH11
DPARAM(5) = DH22
END IF
DPARAM(1) = DFLAG
RETURN
*
* End of DROTMG
*
END
*> \brief \b DSBMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER INCX,INCY,K,LDA,N
* CHARACTER UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*> K is INTEGER
*> On entry, K specifies the number of super-diagonals of the
*> matrix A. K must satisfy 0 .le. K.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. LDA must be at least
*> ( k + 1 ).
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the
*> vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta.
*> \endverbatim
*>
*> \param[in,out] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hbmv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER INCX,INCY,K,LDA,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP1,TEMP2
INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX,MIN
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (K.LT.0) THEN
INFO = 3
ELSE IF (LDA.LT. (K+1)) THEN
INFO = 6
ELSE IF (INCX.EQ.0) THEN
INFO = 8
ELSE IF (INCY.EQ.0) THEN
INFO = 11
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSBMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
*
* Set up the start points in X and Y.
*
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (N-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (N-1)*INCY
END IF
*
* Start the operations. In this version the elements of the array A
* are accessed sequentially with one pass through A.
*
* First form y := beta*y.
*
IF (BETA.NE.ONE) THEN
IF (INCY.EQ.1) THEN
IF (BETA.EQ.ZERO) THEN
DO 10 I = 1,N
Y(I) = ZERO
10 CONTINUE
ELSE
DO 20 I = 1,N
Y(I) = BETA*Y(I)
20 CONTINUE
END IF
ELSE
IY = KY
IF (BETA.EQ.ZERO) THEN
DO 30 I = 1,N
Y(IY) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40 I = 1,N
Y(IY) = BETA*Y(IY)
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF (ALPHA.EQ.ZERO) RETURN
IF (LSAME(UPLO,'U')) THEN
*
* Form y when upper triangle of A is stored.
*
KPLUS1 = K + 1
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 60 J = 1,N
TEMP1 = ALPHA*X(J)
TEMP2 = ZERO
L = KPLUS1 - J
DO 50 I = MAX(1,J-K),J - 1
Y(I) = Y(I) + TEMP1*A(L+I,J)
TEMP2 = TEMP2 + A(L+I,J)*X(I)
50 CONTINUE
Y(J) = Y(J) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
60 CONTINUE
ELSE
JX = KX
JY = KY
DO 80 J = 1,N
TEMP1 = ALPHA*X(JX)
TEMP2 = ZERO
IX = KX
IY = KY
L = KPLUS1 - J
DO 70 I = MAX(1,J-K),J - 1
Y(IY) = Y(IY) + TEMP1*A(L+I,J)
TEMP2 = TEMP2 + A(L+I,J)*X(IX)
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y(JY) = Y(JY) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
IF (J.GT.K) THEN
KX = KX + INCX
KY = KY + INCY
END IF
80 CONTINUE
END IF
ELSE
*
* Form y when lower triangle of A is stored.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 100 J = 1,N
TEMP1 = ALPHA*X(J)
TEMP2 = ZERO
Y(J) = Y(J) + TEMP1*A(1,J)
L = 1 - J
DO 90 I = J + 1,MIN(N,J+K)
Y(I) = Y(I) + TEMP1*A(L+I,J)
TEMP2 = TEMP2 + A(L+I,J)*X(I)
90 CONTINUE
Y(J) = Y(J) + ALPHA*TEMP2
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120 J = 1,N
TEMP1 = ALPHA*X(JX)
TEMP2 = ZERO
Y(JY) = Y(JY) + TEMP1*A(1,J)
L = 1 - J
IX = JX
IY = JY
DO 110 I = J + 1,MIN(N,J+K)
IX = IX + INCX
IY = IY + INCY
Y(IY) = Y(IY) + TEMP1*A(L+I,J)
TEMP2 = TEMP2 + A(L+I,J)*X(IX)
110 CONTINUE
Y(JY) = Y(JY) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
120 CONTINUE
END IF
END IF
*
RETURN
*
* End of DSBMV
*
END
*> \brief \b DSCAL
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSCAL(N,DA,DX,INCX)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION DA
* INTEGER INCX,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSCAL scales a vector by a constant.
*> uses unrolled loops for increment equal to 1.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in] DA
*> \verbatim
*> DA is DOUBLE PRECISION
*> On entry, DA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in,out] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup scal
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 3/93 to return if incx .le. 0.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSCAL(N,DA,DX,INCX)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION DA
INTEGER INCX,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
INTEGER I,M,MP1,NINCX
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER (ONE=1.0D+0)
* ..
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
IF (N.LE.0 .OR. INCX.LE.0 .OR. DA.EQ.ONE) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1
*
*
* clean-up loop
*
M = MOD(N,5)
IF (M.NE.0) THEN
DO I = 1,M
DX(I) = DA*DX(I)
END DO
IF (N.LT.5) RETURN
END IF
MP1 = M + 1
DO I = MP1,N,5
DX(I) = DA*DX(I)
DX(I+1) = DA*DX(I+1)
DX(I+2) = DA*DX(I+2)
DX(I+3) = DA*DX(I+3)
DX(I+4) = DA*DX(I+4)
END DO
ELSE
*
* code for increment not equal to 1
*
NINCX = N*INCX
DO I = 1,NINCX,INCX
DX(I) = DA*DX(I)
END DO
END IF
RETURN
*
* End of DSCAL
*
END
*> \brief \b DSDOT
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* DOUBLE PRECISION FUNCTION DSDOT(N,SX,INCX,SY,INCY)
*
* .. Scalar Arguments ..
* INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
* REAL SX(*),SY(*)
* ..
*
* AUTHORS
* =======
* Lawson, C. L., (JPL), Hanson, R. J., (SNLA),
* Kincaid, D. R., (U. of Texas), Krogh, F. T., (JPL)
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> Compute the inner product of two vectors with extended
*> precision accumulation and result.
*>
*> Returns D.P. dot product accumulated in D.P., for S.P. SX and SY
*> DSDOT = sum for I = 0 to N-1 of SX(LX+I*INCX) * SY(LY+I*INCY),
*> where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
*> defined in a similar way using INCY.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in] SX
*> \verbatim
*> SX is REAL array, dimension(N)
*> single precision vector with N elements
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of SX
*> \endverbatim
*>
*> \param[in] SY
*> \verbatim
*> SY is REAL array, dimension(N)
*> single precision vector with N elements
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> storage spacing between elements of SY
*> \endverbatim
*>
*> \result DSDOT
*> \verbatim
*> DSDOT is DOUBLE PRECISION
*> DSDOT double precision dot product (zero if N.LE.0)
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup dot
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*> \endverbatim
*
*> \par References:
* ================
*>
*> \verbatim
*>
*>
*> C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
*> Krogh, Basic linear algebra subprograms for Fortran
*> usage, Algorithm No. 539, Transactions on Mathematical
*> Software 5, 3 (September 1979), pp. 308-323.
*>
*> REVISION HISTORY (YYMMDD)
*>
*> 791001 DATE WRITTEN
*> 890831 Modified array declarations. (WRB)
*> 890831 REVISION DATE from Version 3.2
*> 891214 Prologue converted to Version 4.0 format. (BAB)
*> 920310 Corrected definition of LX in DESCRIPTION. (WRB)
*> 920501 Reformatted the REFERENCES section. (WRB)
*> 070118 Reformat to LAPACK style (JL)
*> \endverbatim
*>
* =====================================================================
DOUBLE PRECISION FUNCTION DSDOT(N,SX,INCX,SY,INCY)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
REAL SX(*),SY(*)
* ..
*
* Authors:
* ========
* Lawson, C. L., (JPL), Hanson, R. J., (SNLA),
* Kincaid, D. R., (U. of Texas), Krogh, F. T., (JPL)
*
* =====================================================================
*
* .. Local Scalars ..
INTEGER I,KX,KY,NS
* ..
* .. Intrinsic Functions ..
INTRINSIC DBLE
* ..
DSDOT = 0.0D0
IF (N.LE.0) RETURN
IF (INCX.EQ.INCY .AND. INCX.GT.0) THEN
*
* Code for equal, positive, non-unit increments.
*
NS = N*INCX
DO I = 1,NS,INCX
DSDOT = DSDOT + DBLE(SX(I))*DBLE(SY(I))
END DO
ELSE
*
* Code for unequal or nonpositive increments.
*
KX = 1
KY = 1
IF (INCX.LT.0) KX = 1 + (1-N)*INCX
IF (INCY.LT.0) KY = 1 + (1-N)*INCY
DO I = 1,N
DSDOT = DSDOT + DBLE(SX(KX))*DBLE(SY(KY))
KX = KX + INCX
KY = KY + INCY
END DO
END IF
RETURN
*
* End of DSDOT
*
END
*> \brief \b DSPMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER INCX,INCY,N
* CHARACTER UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION AP(*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] AP
*> \verbatim
*> AP is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the n
*> element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta. When BETA is
*> supplied as zero then Y need not be set on input.
*> \endverbatim
*>
*> \param[in,out] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hpmv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER INCX,INCY,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION AP(*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP1,TEMP2
INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (INCX.EQ.0) THEN
INFO = 6
ELSE IF (INCY.EQ.0) THEN
INFO = 9
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSPMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
*
* Set up the start points in X and Y.
*
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (N-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (N-1)*INCY
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
* First form y := beta*y.
*
IF (BETA.NE.ONE) THEN
IF (INCY.EQ.1) THEN
IF (BETA.EQ.ZERO) THEN
DO 10 I = 1,N
Y(I) = ZERO
10 CONTINUE
ELSE
DO 20 I = 1,N
Y(I) = BETA*Y(I)
20 CONTINUE
END IF
ELSE
IY = KY
IF (BETA.EQ.ZERO) THEN
DO 30 I = 1,N
Y(IY) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40 I = 1,N
Y(IY) = BETA*Y(IY)
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF (ALPHA.EQ.ZERO) RETURN
KK = 1
IF (LSAME(UPLO,'U')) THEN
*
* Form y when AP contains the upper triangle.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 60 J = 1,N
TEMP1 = ALPHA*X(J)
TEMP2 = ZERO
K = KK
DO 50 I = 1,J - 1
Y(I) = Y(I) + TEMP1*AP(K)
TEMP2 = TEMP2 + AP(K)*X(I)
K = K + 1
50 CONTINUE
Y(J) = Y(J) + TEMP1*AP(KK+J-1) + ALPHA*TEMP2
KK = KK + J
60 CONTINUE
ELSE
JX = KX
JY = KY
DO 80 J = 1,N
TEMP1 = ALPHA*X(JX)
TEMP2 = ZERO
IX = KX
IY = KY
DO 70 K = KK,KK + J - 2
Y(IY) = Y(IY) + TEMP1*AP(K)
TEMP2 = TEMP2 + AP(K)*X(IX)
IX = IX + INCX
IY = IY + INCY
70 CONTINUE
Y(JY) = Y(JY) + TEMP1*AP(KK+J-1) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
KK = KK + J
80 CONTINUE
END IF
ELSE
*
* Form y when AP contains the lower triangle.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 100 J = 1,N
TEMP1 = ALPHA*X(J)
TEMP2 = ZERO
Y(J) = Y(J) + TEMP1*AP(KK)
K = KK + 1
DO 90 I = J + 1,N
Y(I) = Y(I) + TEMP1*AP(K)
TEMP2 = TEMP2 + AP(K)*X(I)
K = K + 1
90 CONTINUE
Y(J) = Y(J) + ALPHA*TEMP2
KK = KK + (N-J+1)
100 CONTINUE
ELSE
JX = KX
JY = KY
DO 120 J = 1,N
TEMP1 = ALPHA*X(JX)
TEMP2 = ZERO
Y(JY) = Y(JY) + TEMP1*AP(KK)
IX = JX
IY = JY
DO 110 K = KK + 1,KK + N - J
IX = IX + INCX
IY = IY + INCY
Y(IY) = Y(IY) + TEMP1*AP(K)
TEMP2 = TEMP2 + AP(K)*X(IX)
110 CONTINUE
Y(JY) = Y(JY) + ALPHA*TEMP2
JX = JX + INCX
JY = JY + INCY
KK = KK + (N-J+1)
120 CONTINUE
END IF
END IF
*
RETURN
*
* End of DSPMV
*
END
*> \brief \b DSPR2
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA
* INTEGER INCX,INCY,N
* CHARACTER UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION AP(*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSPR2 performs the symmetric rank 2 operation
*>
*> A := alpha*x*y**T + alpha*y*x**T + 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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the n
*> element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCY ) ).
*> Before entry, the incremented array Y must contain the n
*> element vector y.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*>
*> \param[in,out] AP
*> \verbatim
*> AP is DOUBLE PRECISION array, 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.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hpr2
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER INCX,INCY,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION AP(*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (INCX.EQ.0) THEN
INFO = 5
ELSE IF (INCY.EQ.0) THEN
INFO = 7
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSPR2 ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
*
* Set up the start points in X and Y if the increments are not both
* unity.
*
IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (N-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (N-1)*INCY
END IF
JX = KX
JY = KY
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
KK = 1
IF (LSAME(UPLO,'U')) THEN
*
* Form A when upper triangle is stored in AP.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 20 J = 1,N
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.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).NE.ZERO) .OR. (Y(JY).NE.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.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 60 J = 1,N
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.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).NE.ZERO) .OR. (Y(JY).NE.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 of DSPR2
*
END
*> \brief \b DSPR
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA
* INTEGER INCX,N
* CHARACTER UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION AP(*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSPR performs the symmetric rank 1 operation
*>
*> A := alpha*x*x**T + 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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the n
*> element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in,out] AP
*> \verbatim
*> AP is DOUBLE PRECISION array, 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.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hpr
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSPR(UPLO,N,ALPHA,X,INCX,AP)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER INCX,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION AP(*),X(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (INCX.EQ.0) THEN
INFO = 5
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSPR ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
*
* Set the start point in X if the increment is not unity.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of the array AP
* are accessed sequentially with one pass through AP.
*
KK = 1
IF (LSAME(UPLO,'U')) THEN
*
* Form A when upper triangle is stored in AP.
*
IF (INCX.EQ.1) THEN
DO 20 J = 1,N
IF (X(J).NE.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).NE.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.EQ.1) THEN
DO 60 J = 1,N
IF (X(J).NE.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).NE.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 of DSPR
*
END
*> \brief \b DSWAP
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSWAP(N,DX,INCX,DY,INCY)
*
* .. Scalar Arguments ..
* INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*),DY(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSWAP interchanges two vectors.
*> uses unrolled loops for increments equal to 1.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in,out] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*>
*> \param[in,out] DY
*> \verbatim
*> DY is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCY ) )
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> storage spacing between elements of DY
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup swap
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSWAP(N,DX,INCX,DY,INCY)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,INCY,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*),DY(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DTEMP
INTEGER I,IX,IY,M,MP1
* ..
* .. Intrinsic Functions ..
INTRINSIC MOD
* ..
IF (N.LE.0) RETURN
IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN
*
* code for both increments equal to 1
*
*
* clean-up loop
*
M = MOD(N,3)
IF (M.NE.0) THEN
DO I = 1,M
DTEMP = DX(I)
DX(I) = DY(I)
DY(I) = DTEMP
END DO
IF (N.LT.3) RETURN
END IF
MP1 = M + 1
DO I = MP1,N,3
DTEMP = DX(I)
DX(I) = DY(I)
DY(I) = DTEMP
DTEMP = DX(I+1)
DX(I+1) = DY(I+1)
DY(I+1) = DTEMP
DTEMP = DX(I+2)
DX(I+2) = DY(I+2)
DY(I+2) = DTEMP
END DO
ELSE
*
* code for unequal increments or equal increments not equal
* to 1
*
IX = 1
IY = 1
IF (INCX.LT.0) IX = (-N+1)*INCX + 1
IF (INCY.LT.0) IY = (-N+1)*INCY + 1
DO I = 1,N
DTEMP = DX(IX)
DX(IX) = DY(IY)
DY(IY) = DTEMP
IX = IX + INCX
IY = IY + INCY
END DO
END IF
RETURN
*
* End of DSWAP
*
END
*> \brief \b DSYMM
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER LDA,LDB,LDC,M,N
* CHARACTER SIDE,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSYMM performs one of the matrix-matrix operations
*>
*> C := alpha*A*B + beta*C,
*>
*> or
*>
*> C := alpha*B*A + beta*C,
*>
*> where alpha and beta are scalars, A is a symmetric matrix and B and
*> C are m by n matrices.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] SIDE
*> \verbatim
*> SIDE is CHARACTER*1
*> On entry, SIDE specifies whether the symmetric matrix A
*> appears on the left or right in the operation as follows:
*>
*> SIDE = 'L' or 'l' C := alpha*A*B + beta*C,
*>
*> SIDE = 'R' or 'r' C := alpha*B*A + beta*C,
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> On entry, UPLO specifies whether the upper or lower
*> triangular part of the symmetric matrix A is to be
*> referenced as follows:
*>
*> UPLO = 'U' or 'u' Only the upper triangular part of the
*> symmetric matrix is to be referenced.
*>
*> UPLO = 'L' or 'l' Only the lower triangular part of the
*> symmetric matrix is to be referenced.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> On entry, M specifies the number of rows of the matrix C.
*> M must be at least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of columns of the matrix C.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
*> m when SIDE = 'L' or 'l' and is n otherwise.
*> Before entry with SIDE = 'L' or 'l', the m by m part of
*> the array A must contain the symmetric matrix, such that
*> when UPLO = 'U' or 'u', the leading m by m upper triangular
*> part of the array A must contain the upper triangular part
*> of the symmetric matrix and the strictly lower triangular
*> part of A is not referenced, and when UPLO = 'L' or 'l',
*> the leading m by m lower triangular part of the array A
*> must contain the lower triangular part of the symmetric
*> matrix and the strictly upper triangular part of A is not
*> referenced.
*> Before entry with SIDE = 'R' or 'r', the n by n part of
*> the array A must contain the symmetric matrix, such that
*> when UPLO = 'U' or 'u', the leading n by n upper triangular
*> part of the array A must contain the upper triangular part
*> of the symmetric matrix and the strictly lower triangular
*> part of A is not referenced, and when UPLO = 'L' or 'l',
*> the leading n by n lower triangular part of the array A
*> must contain the lower triangular part of the symmetric
*> matrix and the strictly upper triangular part of A is not
*> referenced.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. When SIDE = 'L' or 'l' then
*> LDA must be at least max( 1, m ), otherwise LDA must be at
*> least max( 1, n ).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension ( LDB, N )
*> Before entry, the leading m by n part of the array B must
*> contain the matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> On entry, LDB specifies the first dimension of B as declared
*> in the calling (sub) program. LDB must be at least
*> max( 1, m ).
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta. When BETA is
*> supplied as zero then C need not be set on input.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*> C is DOUBLE PRECISION array, dimension ( LDC, N )
*> Before entry, the leading m by n part of the array C must
*> contain the matrix C, except when beta is zero, in which
*> case C need not be set on entry.
*> On exit, the array C is overwritten by the m by n updated
*> matrix.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*> LDC is INTEGER
*> On entry, LDC specifies the first dimension of C as declared
*> in the calling (sub) program. LDC must be at least
*> max( 1, m ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hemm
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 3 Blas routine.
*>
*> -- Written on 8-February-1989.
*> Jack Dongarra, Argonne National Laboratory.
*> Iain Duff, AERE Harwell.
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
*> Sven Hammarling, Numerical Algorithms Group Ltd.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
* -- Reference BLAS level3 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER LDA,LDB,LDC,M,N
CHARACTER SIDE,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP1,TEMP2
INTEGER I,INFO,J,K,NROWA
LOGICAL UPPER
* ..
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
*
* Set NROWA as the number of rows of A.
*
IF (LSAME(SIDE,'L')) THEN
NROWA = M
ELSE
NROWA = N
END IF
UPPER = LSAME(UPLO,'U')
*
* Test the input parameters.
*
INFO = 0
IF ((.NOT.LSAME(SIDE,'L')) .AND.
+ (.NOT.LSAME(SIDE,'R'))) THEN
INFO = 1
ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
INFO = 2
ELSE IF (M.LT.0) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
INFO = 7
ELSE IF (LDB.LT.MAX(1,M)) THEN
INFO = 9
ELSE IF (LDC.LT.MAX(1,M)) THEN
INFO = 12
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSYMM ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
+ ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
*
* And when alpha.eq.zero.
*
IF (ALPHA.EQ.ZERO) THEN
IF (BETA.EQ.ZERO) THEN
DO 20 J = 1,N
DO 10 I = 1,M
C(I,J) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1,N
DO 30 I = 1,M
C(I,J) = BETA*C(I,J)
30 CONTINUE
40 CONTINUE
END IF
RETURN
END IF
*
* Start the operations.
*
IF (LSAME(SIDE,'L')) THEN
*
* Form C := alpha*A*B + beta*C.
*
IF (UPPER) THEN
DO 70 J = 1,N
DO 60 I = 1,M
TEMP1 = ALPHA*B(I,J)
TEMP2 = ZERO
DO 50 K = 1,I - 1
C(K,J) = C(K,J) + TEMP1*A(K,I)
TEMP2 = TEMP2 + B(K,J)*A(K,I)
50 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
ELSE
C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
+ ALPHA*TEMP2
END IF
60 CONTINUE
70 CONTINUE
ELSE
DO 100 J = 1,N
DO 90 I = M,1,-1
TEMP1 = ALPHA*B(I,J)
TEMP2 = ZERO
DO 80 K = I + 1,M
C(K,J) = C(K,J) + TEMP1*A(K,I)
TEMP2 = TEMP2 + B(K,J)*A(K,I)
80 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
ELSE
C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
+ ALPHA*TEMP2
END IF
90 CONTINUE
100 CONTINUE
END IF
ELSE
*
* Form C := alpha*B*A + beta*C.
*
DO 170 J = 1,N
TEMP1 = ALPHA*A(J,J)
IF (BETA.EQ.ZERO) THEN
DO 110 I = 1,M
C(I,J) = TEMP1*B(I,J)
110 CONTINUE
ELSE
DO 120 I = 1,M
C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
120 CONTINUE
END IF
DO 140 K = 1,J - 1
IF (UPPER) THEN
TEMP1 = ALPHA*A(K,J)
ELSE
TEMP1 = ALPHA*A(J,K)
END IF
DO 130 I = 1,M
C(I,J) = C(I,J) + TEMP1*B(I,K)
130 CONTINUE
140 CONTINUE
DO 160 K = J + 1,N
IF (UPPER) THEN
TEMP1 = ALPHA*A(J,K)
ELSE
TEMP1 = ALPHA*A(K,J)
END IF
DO 150 I = 1,M
C(I,J) = C(I,J) + TEMP1*B(I,K)
150 CONTINUE
160 CONTINUE
170 CONTINUE
END IF
*
RETURN
*
* End of DSYMM
*
END
*> \brief \b DSYMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER INCX,INCY,LDA,N
* CHARACTER UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the n
*> element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta. When BETA is
*> supplied as zero then Y need not be set on input.
*> \endverbatim
*>
*> \param[in,out] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup hemv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER INCX,INCY,LDA,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (LDA.LT.MAX(1,N)) THEN
INFO = 5
ELSE IF (INCX.EQ.0) THEN
INFO = 7
ELSE IF (INCY.EQ.0) THEN
INFO = 10
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSYMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
*
* Set up the start points in X and Y.
*
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (N-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (N-1)*INCY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the triangular part
* of A.
*
* First form y := beta*y.
*
IF (BETA.NE.ONE) THEN
IF (INCY.EQ.1) THEN
IF (BETA.EQ.ZERO) THEN
DO 10 I = 1,N
Y(I) = ZERO
10 CONTINUE
ELSE
DO 20 I = 1,N
Y(I) = BETA*Y(I)
20 CONTINUE
END IF
ELSE
IY = KY
IF (BETA.EQ.ZERO) THEN
DO 30 I = 1,N
Y(IY) = ZERO
IY = IY + INCY
30 CONTINUE
ELSE
DO 40 I = 1,N
Y(IY) = BETA*Y(IY)
IY = IY + INCY
40 CONTINUE
END IF
END IF
END IF
IF (ALPHA.EQ.ZERO) RETURN
IF (LSAME(UPLO,'U')) THEN
*
* Form y when A is stored in upper triangle.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 60 J = 1,N
TEMP1 = ALPHA*X(J)
TEMP2 = ZERO
DO 50 I = 1,J - 1
Y(I) = Y(I) + TEMP1*A(I,J)
TEMP2 = TEMP2 + 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.EQ.1) .AND. (INCY.EQ.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 of DSYMV
*
END
*> \brief \b DSYR2
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSYR2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA
* INTEGER INCX,INCY,LDA,N
* CHARACTER UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSYR2 performs the symmetric rank 2 operation
*>
*> A := alpha*x*y**T + alpha*y*x**T + A,
*>
*> where alpha is a scalar, x and y are n element vectors and A is an n
*> by n symmetric matrix.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the n
*> element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCY ) ).
*> Before entry, the incremented array Y must contain the n
*> element vector y.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*> INCY is INTEGER
*> On entry, INCY specifies the increment for the elements of
*> Y. INCY must not be zero.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup her2
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSYR2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER INCX,INCY,LDA,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*),Y(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (INCX.EQ.0) THEN
INFO = 5
ELSE IF (INCY.EQ.0) THEN
INFO = 7
ELSE IF (LDA.LT.MAX(1,N)) THEN
INFO = 9
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSYR2 ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
*
* Set up the start points in X and Y if the increments are not both
* unity.
*
IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
IF (INCX.GT.0) THEN
KX = 1
ELSE
KX = 1 - (N-1)*INCX
END IF
IF (INCY.GT.0) THEN
KY = 1
ELSE
KY = 1 - (N-1)*INCY
END IF
JX = KX
JY = KY
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the triangular part
* of A.
*
IF (LSAME(UPLO,'U')) THEN
*
* Form A when A is stored in the upper triangle.
*
IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 20 J = 1,N
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.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).NE.ZERO) .OR. (Y(JY).NE.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.EQ.1) .AND. (INCY.EQ.1)) THEN
DO 60 J = 1,N
IF ((X(J).NE.ZERO) .OR. (Y(J).NE.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).NE.ZERO) .OR. (Y(JY).NE.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 of DSYR2
*
END
*> \brief \b DSYR2K
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER K,LDA,LDB,LDC,N
* CHARACTER TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSYR2K performs one of the symmetric rank 2k operations
*>
*> C := alpha*A*B**T + alpha*B*A**T + beta*C,
*>
*> or
*>
*> C := alpha*A**T*B + alpha*B**T*A + beta*C,
*>
*> where alpha and beta are scalars, C is an n by n symmetric matrix
*> and A and B are n by k matrices in the first case and k by n
*> matrices in the second case.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> On entry, UPLO specifies whether the upper or lower
*> triangular part of the array C is to be referenced as
*> follows:
*>
*> UPLO = 'U' or 'u' Only the upper triangular part of C
*> is to be referenced.
*>
*> UPLO = 'L' or 'l' Only the lower triangular part of C
*> is to be referenced.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is CHARACTER*1
*> On entry, TRANS specifies the operation to be performed as
*> follows:
*>
*> TRANS = 'N' or 'n' C := alpha*A*B**T + alpha*B*A**T +
*> beta*C.
*>
*> TRANS = 'T' or 't' C := alpha*A**T*B + alpha*B**T*A +
*> beta*C.
*>
*> TRANS = 'C' or 'c' C := alpha*A**T*B + alpha*B**T*A +
*> beta*C.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix C. N must be
*> at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*> K is INTEGER
*> On entry with TRANS = 'N' or 'n', K specifies the number
*> of columns of the matrices A and B, and on entry with
*> TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
*> of rows of the matrices A and B. K must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
*> k when TRANS = 'N' or 'n', and is n otherwise.
*> Before entry with TRANS = 'N' or 'n', the leading n by k
*> part of the array A must contain the matrix A, otherwise
*> the leading k by n part of the array A must contain the
*> matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. When TRANS = 'N' or 'n'
*> then LDA must be at least max( 1, n ), otherwise LDA must
*> be at least max( 1, k ).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension ( LDB, kb ), where kb is
*> k when TRANS = 'N' or 'n', and is n otherwise.
*> Before entry with TRANS = 'N' or 'n', the leading n by k
*> part of the array B must contain the matrix B, otherwise
*> the leading k by n part of the array B must contain the
*> matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> On entry, LDB specifies the first dimension of B as declared
*> in the calling (sub) program. When TRANS = 'N' or 'n'
*> then LDB must be at least max( 1, n ), otherwise LDB must
*> be at least max( 1, k ).
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*> C is DOUBLE PRECISION array, dimension ( LDC, N )
*> Before entry with UPLO = 'U' or 'u', the leading n by n
*> upper triangular part of the array C must contain the upper
*> triangular part of the symmetric matrix and the strictly
*> lower triangular part of C is not referenced. On exit, the
*> upper triangular part of the array C is overwritten by the
*> upper triangular part of the updated matrix.
*> Before entry with UPLO = 'L' or 'l', the leading n by n
*> lower triangular part of the array C must contain the lower
*> triangular part of the symmetric matrix and the strictly
*> upper triangular part of C is not referenced. On exit, the
*> lower triangular part of the array C is overwritten by the
*> lower triangular part of the updated matrix.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*> LDC is INTEGER
*> On entry, LDC specifies the first dimension of C as declared
*> in the calling (sub) program. LDC must be at least
*> max( 1, n ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup her2k
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 3 Blas routine.
*>
*>
*> -- Written on 8-February-1989.
*> Jack Dongarra, Argonne National Laboratory.
*> Iain Duff, AERE Harwell.
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
*> Sven Hammarling, Numerical Algorithms Group Ltd.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
* -- Reference BLAS level3 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER K,LDA,LDB,LDC,N
CHARACTER TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP1,TEMP2
INTEGER I,INFO,J,L,NROWA
LOGICAL UPPER
* ..
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
*
* Test the input parameters.
*
IF (LSAME(TRANS,'N')) THEN
NROWA = N
ELSE
NROWA = K
END IF
UPPER = LSAME(UPLO,'U')
*
INFO = 0
IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
INFO = 1
ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
+ (.NOT.LSAME(TRANS,'T')) .AND.
+ (.NOT.LSAME(TRANS,'C'))) THEN
INFO = 2
ELSE IF (N.LT.0) THEN
INFO = 3
ELSE IF (K.LT.0) THEN
INFO = 4
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
INFO = 7
ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
INFO = 9
ELSE IF (LDC.LT.MAX(1,N)) THEN
INFO = 12
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSYR2K',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
+ (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
*
* And when alpha.eq.zero.
*
IF (ALPHA.EQ.ZERO) THEN
IF (UPPER) THEN
IF (BETA.EQ.ZERO) THEN
DO 20 J = 1,N
DO 10 I = 1,J
C(I,J) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1,N
DO 30 I = 1,J
C(I,J) = BETA*C(I,J)
30 CONTINUE
40 CONTINUE
END IF
ELSE
IF (BETA.EQ.ZERO) THEN
DO 60 J = 1,N
DO 50 I = J,N
C(I,J) = ZERO
50 CONTINUE
60 CONTINUE
ELSE
DO 80 J = 1,N
DO 70 I = J,N
C(I,J) = BETA*C(I,J)
70 CONTINUE
80 CONTINUE
END IF
END IF
RETURN
END IF
*
* Start the operations.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form C := alpha*A*B**T + alpha*B*A**T + C.
*
IF (UPPER) THEN
DO 130 J = 1,N
IF (BETA.EQ.ZERO) THEN
DO 90 I = 1,J
C(I,J) = ZERO
90 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 100 I = 1,J
C(I,J) = BETA*C(I,J)
100 CONTINUE
END IF
DO 120 L = 1,K
IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
TEMP1 = ALPHA*B(J,L)
TEMP2 = ALPHA*A(J,L)
DO 110 I = 1,J
C(I,J) = C(I,J) + A(I,L)*TEMP1 +
+ B(I,L)*TEMP2
110 CONTINUE
END IF
120 CONTINUE
130 CONTINUE
ELSE
DO 180 J = 1,N
IF (BETA.EQ.ZERO) THEN
DO 140 I = J,N
C(I,J) = ZERO
140 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 150 I = J,N
C(I,J) = BETA*C(I,J)
150 CONTINUE
END IF
DO 170 L = 1,K
IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
TEMP1 = ALPHA*B(J,L)
TEMP2 = ALPHA*A(J,L)
DO 160 I = J,N
C(I,J) = C(I,J) + A(I,L)*TEMP1 +
+ B(I,L)*TEMP2
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
END IF
ELSE
*
* Form C := alpha*A**T*B + alpha*B**T*A + C.
*
IF (UPPER) THEN
DO 210 J = 1,N
DO 200 I = 1,J
TEMP1 = ZERO
TEMP2 = ZERO
DO 190 L = 1,K
TEMP1 = TEMP1 + A(L,I)*B(L,J)
TEMP2 = TEMP2 + B(L,I)*A(L,J)
190 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
ELSE
C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
+ ALPHA*TEMP2
END IF
200 CONTINUE
210 CONTINUE
ELSE
DO 240 J = 1,N
DO 230 I = J,N
TEMP1 = ZERO
TEMP2 = ZERO
DO 220 L = 1,K
TEMP1 = TEMP1 + A(L,I)*B(L,J)
TEMP2 = TEMP2 + B(L,I)*A(L,J)
220 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
ELSE
C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
+ ALPHA*TEMP2
END IF
230 CONTINUE
240 CONTINUE
END IF
END IF
*
RETURN
*
* End of DSYR2K
*
END
*> \brief \b DSYR
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSYR(UPLO,N,ALPHA,X,INCX,A,LDA)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA
* INTEGER INCX,LDA,N
* CHARACTER UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSYR performs the symmetric rank 1 operation
*>
*> A := alpha*x*x**T + A,
*>
*> where alpha is a real scalar, x is an n element vector and A is an
*> n by n symmetric matrix.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*> X is DOUBLE PRECISION array, dimension at least
*> ( 1 + ( n - 1 )*abs( INCX ) ).
*> Before entry, the incremented array X must contain the n
*> element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup her
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSYR(UPLO,N,ALPHA,X,INCX,A,LDA)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER INCX,LDA,N
CHARACTER UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (N.LT.0) THEN
INFO = 2
ELSE IF (INCX.EQ.0) THEN
INFO = 5
ELSE IF (LDA.LT.MAX(1,N)) THEN
INFO = 7
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSYR ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
*
* Set the start point in X if the increment is not unity.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through the triangular part
* of A.
*
IF (LSAME(UPLO,'U')) THEN
*
* Form A when A is stored in upper triangle.
*
IF (INCX.EQ.1) THEN
DO 20 J = 1,N
IF (X(J).NE.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).NE.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.EQ.1) THEN
DO 60 J = 1,N
IF (X(J).NE.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).NE.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 of DSYR
*
END
*> \brief \b DSYRK
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA,BETA
* INTEGER K,LDA,LDC,N
* CHARACTER TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),C(LDC,*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DSYRK performs one of the symmetric rank k operations
*>
*> C := alpha*A*A**T + beta*C,
*>
*> or
*>
*> C := alpha*A**T*A + beta*C,
*>
*> where alpha and beta are scalars, C is an n by n symmetric matrix
*> and A is an n by k matrix in the first case and a k by n matrix
*> in the second case.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> On entry, UPLO specifies whether the upper or lower
*> triangular part of the array C is to be referenced as
*> follows:
*>
*> UPLO = 'U' or 'u' Only the upper triangular part of C
*> is to be referenced.
*>
*> UPLO = 'L' or 'l' Only the lower triangular part of C
*> is to be referenced.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is CHARACTER*1
*> On entry, TRANS specifies the operation to be performed as
*> follows:
*>
*> TRANS = 'N' or 'n' C := alpha*A*A**T + beta*C.
*>
*> TRANS = 'T' or 't' C := alpha*A**T*A + beta*C.
*>
*> TRANS = 'C' or 'c' C := alpha*A**T*A + beta*C.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix C. N must be
*> at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*> K is INTEGER
*> On entry with TRANS = 'N' or 'n', K specifies the number
*> of columns of the matrix A, and on entry with
*> TRANS = 'T' or 't' or 'C' or 'c', K specifies the number
*> of rows of the matrix A. K must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, ka ), where ka is
*> k when TRANS = 'N' or 'n', and is n otherwise.
*> Before entry with TRANS = 'N' or 'n', the leading n by k
*> part of the array A must contain the matrix A, otherwise
*> the leading k by n part of the array A must contain the
*> matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. When TRANS = 'N' or 'n'
*> then LDA must be at least max( 1, n ), otherwise LDA must
*> be at least max( 1, k ).
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*> BETA is DOUBLE PRECISION.
*> On entry, BETA specifies the scalar beta.
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*> C is DOUBLE PRECISION array, dimension ( LDC, N )
*> Before entry with UPLO = 'U' or 'u', the leading n by n
*> upper triangular part of the array C must contain the upper
*> triangular part of the symmetric matrix and the strictly
*> lower triangular part of C is not referenced. On exit, the
*> upper triangular part of the array C is overwritten by the
*> upper triangular part of the updated matrix.
*> Before entry with UPLO = 'L' or 'l', the leading n by n
*> lower triangular part of the array C must contain the lower
*> triangular part of the symmetric matrix and the strictly
*> upper triangular part of C is not referenced. On exit, the
*> lower triangular part of the array C is overwritten by the
*> lower triangular part of the updated matrix.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*> LDC is INTEGER
*> On entry, LDC specifies the first dimension of C as declared
*> in the calling (sub) program. LDC must be at least
*> max( 1, n ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup herk
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 3 Blas routine.
*>
*> -- Written on 8-February-1989.
*> Jack Dongarra, Argonne National Laboratory.
*> Iain Duff, AERE Harwell.
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
*> Sven Hammarling, Numerical Algorithms Group Ltd.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DSYRK(UPLO,TRANS,N,K,ALPHA,A,LDA,BETA,C,LDC)
*
* -- Reference BLAS level3 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA,BETA
INTEGER K,LDA,LDC,N
CHARACTER TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),C(LDC,*)
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I,INFO,J,L,NROWA
LOGICAL UPPER
* ..
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
*
* Test the input parameters.
*
IF (LSAME(TRANS,'N')) THEN
NROWA = N
ELSE
NROWA = K
END IF
UPPER = LSAME(UPLO,'U')
*
INFO = 0
IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
INFO = 1
ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
+ (.NOT.LSAME(TRANS,'T')) .AND.
+ (.NOT.LSAME(TRANS,'C'))) THEN
INFO = 2
ELSE IF (N.LT.0) THEN
INFO = 3
ELSE IF (K.LT.0) THEN
INFO = 4
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
INFO = 7
ELSE IF (LDC.LT.MAX(1,N)) THEN
INFO = 10
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DSYRK ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
+ (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
*
* And when alpha.eq.zero.
*
IF (ALPHA.EQ.ZERO) THEN
IF (UPPER) THEN
IF (BETA.EQ.ZERO) THEN
DO 20 J = 1,N
DO 10 I = 1,J
C(I,J) = ZERO
10 CONTINUE
20 CONTINUE
ELSE
DO 40 J = 1,N
DO 30 I = 1,J
C(I,J) = BETA*C(I,J)
30 CONTINUE
40 CONTINUE
END IF
ELSE
IF (BETA.EQ.ZERO) THEN
DO 60 J = 1,N
DO 50 I = J,N
C(I,J) = ZERO
50 CONTINUE
60 CONTINUE
ELSE
DO 80 J = 1,N
DO 70 I = J,N
C(I,J) = BETA*C(I,J)
70 CONTINUE
80 CONTINUE
END IF
END IF
RETURN
END IF
*
* Start the operations.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form C := alpha*A*A**T + beta*C.
*
IF (UPPER) THEN
DO 130 J = 1,N
IF (BETA.EQ.ZERO) THEN
DO 90 I = 1,J
C(I,J) = ZERO
90 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 100 I = 1,J
C(I,J) = BETA*C(I,J)
100 CONTINUE
END IF
DO 120 L = 1,K
IF (A(J,L).NE.ZERO) THEN
TEMP = ALPHA*A(J,L)
DO 110 I = 1,J
C(I,J) = C(I,J) + TEMP*A(I,L)
110 CONTINUE
END IF
120 CONTINUE
130 CONTINUE
ELSE
DO 180 J = 1,N
IF (BETA.EQ.ZERO) THEN
DO 140 I = J,N
C(I,J) = ZERO
140 CONTINUE
ELSE IF (BETA.NE.ONE) THEN
DO 150 I = J,N
C(I,J) = BETA*C(I,J)
150 CONTINUE
END IF
DO 170 L = 1,K
IF (A(J,L).NE.ZERO) THEN
TEMP = ALPHA*A(J,L)
DO 160 I = J,N
C(I,J) = C(I,J) + TEMP*A(I,L)
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
END IF
ELSE
*
* Form C := alpha*A**T*A + beta*C.
*
IF (UPPER) THEN
DO 210 J = 1,N
DO 200 I = 1,J
TEMP = ZERO
DO 190 L = 1,K
TEMP = TEMP + A(L,I)*A(L,J)
190 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP
ELSE
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
END IF
200 CONTINUE
210 CONTINUE
ELSE
DO 240 J = 1,N
DO 230 I = J,N
TEMP = ZERO
DO 220 L = 1,K
TEMP = TEMP + A(L,I)*A(L,J)
220 CONTINUE
IF (BETA.EQ.ZERO) THEN
C(I,J) = ALPHA*TEMP
ELSE
C(I,J) = ALPHA*TEMP + BETA*C(I,J)
END IF
230 CONTINUE
240 CONTINUE
END IF
END IF
*
RETURN
*
* End of DSYRK
*
END
*> \brief \b DTBMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,K,LDA,N
* CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTBMV performs one of the matrix-vector operations
*>
*> x := A*x, or x := A**T*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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x.
*>
*> TRANS = 'C' or 'c' x := A**T*x.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*> K is 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.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. LDA must be at least
*> ( k + 1 ).
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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
*> transformed vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tbmv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,K,LDA,N
CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D+0)
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L
LOGICAL NOUNIT
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX,MIN
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (.NOT.LSAME(TRANS,'N') .AND.
+ .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 2
ELSE IF (.NOT.LSAME(DIAG,'U') .AND.
+ .NOT.LSAME(DIAG,'N')) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (K.LT.0) THEN
INFO = 5
ELSE IF (LDA.LT. (K+1)) THEN
INFO = 7
ELSE IF (INCX.EQ.0) THEN
INFO = 9
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTBMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (N.EQ.0) RETURN
*
NOUNIT = LSAME(DIAG,'N')
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form x := A*x.
*
IF (LSAME(UPLO,'U')) THEN
KPLUS1 = K + 1
IF (INCX.EQ.1) THEN
DO 20 J = 1,N
IF (X(J).NE.ZERO) THEN
TEMP = X(J)
L = KPLUS1 - J
DO 10 I = MAX(1,J-K),J - 1
X(I) = X(I) + TEMP*A(L+I,J)
10 CONTINUE
IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J)
END IF
20 CONTINUE
ELSE
JX = KX
DO 40 J = 1,N
IF (X(JX).NE.ZERO) THEN
TEMP = X(JX)
IX = KX
L = KPLUS1 - J
DO 30 I = MAX(1,J-K),J - 1
X(IX) = X(IX) + TEMP*A(L+I,J)
IX = IX + INCX
30 CONTINUE
IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J)
END IF
JX = JX + INCX
IF (J.GT.K) KX = KX + INCX
40 CONTINUE
END IF
ELSE
IF (INCX.EQ.1) THEN
DO 60 J = N,1,-1
IF (X(J).NE.ZERO) THEN
TEMP = X(J)
L = 1 - J
DO 50 I = MIN(N,J+K),J + 1,-1
X(I) = X(I) + TEMP*A(L+I,J)
50 CONTINUE
IF (NOUNIT) X(J) = X(J)*A(1,J)
END IF
60 CONTINUE
ELSE
KX = KX + (N-1)*INCX
JX = KX
DO 80 J = N,1,-1
IF (X(JX).NE.ZERO) THEN
TEMP = X(JX)
IX = KX
L = 1 - J
DO 70 I = MIN(N,J+K),J + 1,-1
X(IX) = X(IX) + TEMP*A(L+I,J)
IX = IX - INCX
70 CONTINUE
IF (NOUNIT) X(JX) = X(JX)*A(1,J)
END IF
JX = JX - INCX
IF ((N-J).GE.K) KX = KX - INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := A**T*x.
*
IF (LSAME(UPLO,'U')) THEN
KPLUS1 = K + 1
IF (INCX.EQ.1) THEN
DO 100 J = N,1,-1
TEMP = X(J)
L = KPLUS1 - J
IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
DO 90 I = J - 1,MAX(1,J-K),-1
TEMP = TEMP + A(L+I,J)*X(I)
90 CONTINUE
X(J) = TEMP
100 CONTINUE
ELSE
KX = KX + (N-1)*INCX
JX = KX
DO 120 J = N,1,-1
TEMP = X(JX)
KX = KX - INCX
IX = KX
L = KPLUS1 - J
IF (NOUNIT) TEMP = TEMP*A(KPLUS1,J)
DO 110 I = J - 1,MAX(1,J-K),-1
TEMP = TEMP + A(L+I,J)*X(IX)
IX = IX - INCX
110 CONTINUE
X(JX) = TEMP
JX = JX - INCX
120 CONTINUE
END IF
ELSE
IF (INCX.EQ.1) THEN
DO 140 J = 1,N
TEMP = X(J)
L = 1 - J
IF (NOUNIT) TEMP = TEMP*A(1,J)
DO 130 I = J + 1,MIN(N,J+K)
TEMP = TEMP + A(L+I,J)*X(I)
130 CONTINUE
X(J) = TEMP
140 CONTINUE
ELSE
JX = KX
DO 160 J = 1,N
TEMP = X(JX)
KX = KX + INCX
IX = KX
L = 1 - J
IF (NOUNIT) TEMP = TEMP*A(1,J)
DO 150 I = J + 1,MIN(N,J+K)
TEMP = TEMP + A(L+I,J)*X(IX)
IX = IX + INCX
150 CONTINUE
X(JX) = TEMP
JX = JX + INCX
160 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTBMV
*
END
*> \brief \b DTBSV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,K,LDA,N
* CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTBSV solves one of the systems of equations
*>
*> A*x = b, or A**T*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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x = b.
*>
*> TRANS = 'C' or 'c' A**T*x = b.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*> K is 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.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. LDA must be at least
*> ( k + 1 ).
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tbsv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,K,LDA,N
CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER (ZERO=0.0D+0)
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I,INFO,IX,J,JX,KPLUS1,KX,L
LOGICAL NOUNIT
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX,MIN
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (.NOT.LSAME(TRANS,'N') .AND.
+ .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 2
ELSE IF (.NOT.LSAME(DIAG,'U') .AND.
+ .NOT.LSAME(DIAG,'N')) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (K.LT.0) THEN
INFO = 5
ELSE IF (LDA.LT. (K+1)) THEN
INFO = 7
ELSE IF (INCX.EQ.0) THEN
INFO = 9
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTBSV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (N.EQ.0) RETURN
*
NOUNIT = LSAME(DIAG,'N')
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed by sequentially with one pass through A.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form x := inv( A )*x.
*
IF (LSAME(UPLO,'U')) THEN
KPLUS1 = K + 1
IF (INCX.EQ.1) THEN
DO 20 J = N,1,-1
IF (X(J).NE.ZERO) THEN
L = KPLUS1 - J
IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J)
TEMP = X(J)
DO 10 I = J - 1,MAX(1,J-K),-1
X(I) = X(I) - TEMP*A(L+I,J)
10 CONTINUE
END IF
20 CONTINUE
ELSE
KX = KX + (N-1)*INCX
JX = KX
DO 40 J = N,1,-1
KX = KX - INCX
IF (X(JX).NE.ZERO) THEN
IX = KX
L = KPLUS1 - J
IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J)
TEMP = X(JX)
DO 30 I = J - 1,MAX(1,J-K),-1
X(IX) = X(IX) - TEMP*A(L+I,J)
IX = IX - INCX
30 CONTINUE
END IF
JX = JX - INCX
40 CONTINUE
END IF
ELSE
IF (INCX.EQ.1) THEN
DO 60 J = 1,N
IF (X(J).NE.ZERO) THEN
L = 1 - J
IF (NOUNIT) X(J) = X(J)/A(1,J)
TEMP = X(J)
DO 50 I = J + 1,MIN(N,J+K)
X(I) = X(I) - TEMP*A(L+I,J)
50 CONTINUE
END IF
60 CONTINUE
ELSE
JX = KX
DO 80 J = 1,N
KX = KX + INCX
IF (X(JX).NE.ZERO) THEN
IX = KX
L = 1 - J
IF (NOUNIT) X(JX) = X(JX)/A(1,J)
TEMP = X(JX)
DO 70 I = J + 1,MIN(N,J+K)
X(IX) = X(IX) - TEMP*A(L+I,J)
IX = IX + INCX
70 CONTINUE
END IF
JX = JX + INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := inv( A**T)*x.
*
IF (LSAME(UPLO,'U')) THEN
KPLUS1 = K + 1
IF (INCX.EQ.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) TEMP = TEMP/A(KPLUS1,J)
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) TEMP = TEMP/A(KPLUS1,J)
X(JX) = TEMP
JX = JX + INCX
IF (J.GT.K) KX = KX + INCX
120 CONTINUE
END IF
ELSE
IF (INCX.EQ.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) TEMP = TEMP/A(1,J)
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) TEMP = TEMP/A(1,J)
X(JX) = TEMP
JX = JX - INCX
IF ((N-J).GE.K) KX = KX - INCX
160 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTBSV
*
END
*> \brief \b DTPMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,N
* CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION AP(*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTPMV performs one of the matrix-vector operations
*>
*> x := A*x, or x := A**T*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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x.
*>
*> TRANS = 'C' or 'c' x := A**T*x.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] AP
*> \verbatim
*> AP is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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
*> transformed vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tpmv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,N
CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION AP(*),X(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (.NOT.LSAME(TRANS,'N') .AND.
+ .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 2
ELSE IF (.NOT.LSAME(DIAG,'U') .AND.
+ .NOT.LSAME(DIAG,'N')) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (INCX.EQ.0) THEN
INFO = 7
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTPMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (N.EQ.0) RETURN
*
NOUNIT = LSAME(DIAG,'N')
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of AP are
* accessed sequentially with one pass through AP.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form x:= A*x.
*
IF (LSAME(UPLO,'U')) THEN
KK = 1
IF (INCX.EQ.1) THEN
DO 20 J = 1,N
IF (X(J).NE.ZERO) THEN
TEMP = X(J)
K = KK
DO 10 I = 1,J - 1
X(I) = X(I) + TEMP*AP(K)
K = K + 1
10 CONTINUE
IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
END IF
KK = KK + J
20 CONTINUE
ELSE
JX = KX
DO 40 J = 1,N
IF (X(JX).NE.ZERO) THEN
TEMP = X(JX)
IX = KX
DO 30 K = KK,KK + J - 2
X(IX) = X(IX) + TEMP*AP(K)
IX = IX + INCX
30 CONTINUE
IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
END IF
JX = JX + INCX
KK = KK + J
40 CONTINUE
END IF
ELSE
KK = (N* (N+1))/2
IF (INCX.EQ.1) THEN
DO 60 J = N,1,-1
IF (X(J).NE.ZERO) THEN
TEMP = X(J)
K = KK
DO 50 I = N,J + 1,-1
X(I) = X(I) + TEMP*AP(K)
K = K - 1
50 CONTINUE
IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
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).NE.ZERO) THEN
TEMP = X(JX)
IX = KX
DO 70 K = KK,KK - (N- (J+1)),-1
X(IX) = X(IX) + TEMP*AP(K)
IX = IX - INCX
70 CONTINUE
IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
END IF
JX = JX - INCX
KK = KK - (N-J+1)
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := A**T*x.
*
IF (LSAME(UPLO,'U')) THEN
KK = (N* (N+1))/2
IF (INCX.EQ.1) THEN
DO 100 J = N,1,-1
TEMP = X(J)
IF (NOUNIT) TEMP = TEMP*AP(KK)
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) TEMP = TEMP*AP(KK)
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.EQ.1) THEN
DO 140 J = 1,N
TEMP = X(J)
IF (NOUNIT) TEMP = TEMP*AP(KK)
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) TEMP = TEMP*AP(KK)
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 of DTPMV
*
END
*> \brief \b DTPSV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,N
* CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION AP(*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTPSV solves one of the systems of equations
*>
*> A*x = b, or A**T*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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x = b.
*>
*> TRANS = 'C' or 'c' A**T*x = b.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] AP
*> \verbatim
*> AP is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup tpsv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,N
CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION AP(*),X(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (.NOT.LSAME(TRANS,'N') .AND.
+ .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 2
ELSE IF (.NOT.LSAME(DIAG,'U') .AND.
+ .NOT.LSAME(DIAG,'N')) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (INCX.EQ.0) THEN
INFO = 7
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTPSV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (N.EQ.0) RETURN
*
NOUNIT = LSAME(DIAG,'N')
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of AP are
* accessed sequentially with one pass through AP.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form x := inv( A )*x.
*
IF (LSAME(UPLO,'U')) THEN
KK = (N* (N+1))/2
IF (INCX.EQ.1) THEN
DO 20 J = N,1,-1
IF (X(J).NE.ZERO) THEN
IF (NOUNIT) X(J) = X(J)/AP(KK)
TEMP = X(J)
K = KK - 1
DO 10 I = J - 1,1,-1
X(I) = X(I) - TEMP*AP(K)
K = K - 1
10 CONTINUE
END IF
KK = KK - J
20 CONTINUE
ELSE
JX = KX + (N-1)*INCX
DO 40 J = N,1,-1
IF (X(JX).NE.ZERO) THEN
IF (NOUNIT) X(JX) = X(JX)/AP(KK)
TEMP = X(JX)
IX = JX
DO 30 K = KK - 1,KK - J + 1,-1
IX = IX - INCX
X(IX) = X(IX) - TEMP*AP(K)
30 CONTINUE
END IF
JX = JX - INCX
KK = KK - J
40 CONTINUE
END IF
ELSE
KK = 1
IF (INCX.EQ.1) THEN
DO 60 J = 1,N
IF (X(J).NE.ZERO) THEN
IF (NOUNIT) X(J) = X(J)/AP(KK)
TEMP = X(J)
K = KK + 1
DO 50 I = J + 1,N
X(I) = X(I) - TEMP*AP(K)
K = K + 1
50 CONTINUE
END IF
KK = KK + (N-J+1)
60 CONTINUE
ELSE
JX = KX
DO 80 J = 1,N
IF (X(JX).NE.ZERO) THEN
IF (NOUNIT) X(JX) = X(JX)/AP(KK)
TEMP = X(JX)
IX = JX
DO 70 K = KK + 1,KK + N - J
IX = IX + INCX
X(IX) = X(IX) - TEMP*AP(K)
70 CONTINUE
END IF
JX = JX + INCX
KK = KK + (N-J+1)
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := inv( A**T )*x.
*
IF (LSAME(UPLO,'U')) THEN
KK = 1
IF (INCX.EQ.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) TEMP = TEMP/AP(KK+J-1)
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) TEMP = TEMP/AP(KK+J-1)
X(JX) = TEMP
JX = JX + INCX
KK = KK + J
120 CONTINUE
END IF
ELSE
KK = (N* (N+1))/2
IF (INCX.EQ.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) TEMP = TEMP/AP(KK-N+J)
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) TEMP = TEMP/AP(KK-N+J)
X(JX) = TEMP
JX = JX - INCX
KK = KK - (N-J+1)
160 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTPSV
*
END
*> \brief \b DTRMM
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA
* INTEGER LDA,LDB,M,N
* CHARACTER DIAG,SIDE,TRANSA,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),B(LDB,*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTRMM performs one of the matrix-matrix operations
*>
*> B := alpha*op( A )*B, or B := alpha*B*op( A ),
*>
*> where alpha is a scalar, B is an m by n matrix, A is a unit, or
*> non-unit, upper or lower triangular matrix and op( A ) is one of
*>
*> op( A ) = A or op( A ) = A**T.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] SIDE
*> \verbatim
*> SIDE is CHARACTER*1
*> On entry, SIDE specifies whether op( A ) multiplies B from
*> the left or right as follows:
*>
*> SIDE = 'L' or 'l' B := alpha*op( A )*B.
*>
*> SIDE = 'R' or 'r' B := alpha*B*op( A ).
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> On entry, UPLO specifies whether the matrix A is an upper or
*> lower triangular matrix as follows:
*>
*> UPLO = 'U' or 'u' A is an upper triangular matrix.
*>
*> UPLO = 'L' or 'l' A is a lower triangular matrix.
*> \endverbatim
*>
*> \param[in] TRANSA
*> \verbatim
*> TRANSA is CHARACTER*1
*> On entry, TRANSA specifies the form of op( A ) to be used in
*> the matrix multiplication as follows:
*>
*> TRANSA = 'N' or 'n' op( A ) = A.
*>
*> TRANSA = 'T' or 't' op( A ) = A**T.
*>
*> TRANSA = 'C' or 'c' op( A ) = A**T.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> On entry, M specifies the number of rows of B. M must be at
*> least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of columns of B. N must be
*> at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha. When alpha is
*> zero then A is not referenced and B need not be set before
*> entry.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, k ), where k is m
*> when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
*> Before entry with UPLO = 'U' or 'u', the leading k by k
*> upper triangular part of the array A must contain the upper
*> triangular matrix and the strictly lower triangular part of
*> A is not referenced.
*> Before entry with UPLO = 'L' or 'l', the leading k by k
*> lower triangular part of the array A must contain the lower
*> triangular matrix and the strictly upper triangular part of
*> A is not referenced.
*> Note that when DIAG = 'U' or 'u', the diagonal elements of
*> A are not referenced either, but are assumed to be unity.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. When SIDE = 'L' or 'l' then
*> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
*> then LDA must be at least max( 1, n ).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension ( LDB, N )
*> Before entry, the leading m by n part of the array B must
*> contain the matrix B, and on exit is overwritten by the
*> transformed matrix.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> On entry, LDB specifies the first dimension of B as declared
*> in the calling (sub) program. LDB must be at least
*> max( 1, m ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup trmm
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 3 Blas routine.
*>
*> -- Written on 8-February-1989.
*> Jack Dongarra, Argonne National Laboratory.
*> Iain Duff, AERE Harwell.
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
*> Sven Hammarling, Numerical Algorithms Group Ltd.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*
* -- Reference BLAS level3 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER LDA,LDB,M,N
CHARACTER DIAG,SIDE,TRANSA,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),B(LDB,*)
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I,INFO,J,K,NROWA
LOGICAL LSIDE,NOUNIT,UPPER
* ..
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
*
* Test the input parameters.
*
LSIDE = LSAME(SIDE,'L')
IF (LSIDE) THEN
NROWA = M
ELSE
NROWA = N
END IF
NOUNIT = LSAME(DIAG,'N')
UPPER = LSAME(UPLO,'U')
*
INFO = 0
IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
INFO = 1
ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
INFO = 2
ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
+ (.NOT.LSAME(TRANSA,'T')) .AND.
+ (.NOT.LSAME(TRANSA,'C'))) THEN
INFO = 3
ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND.
+ (.NOT.LSAME(DIAG,'N'))) THEN
INFO = 4
ELSE IF (M.LT.0) THEN
INFO = 5
ELSE IF (N.LT.0) THEN
INFO = 6
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
INFO = 9
ELSE IF (LDB.LT.MAX(1,M)) THEN
INFO = 11
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTRMM ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (M.EQ.0 .OR. N.EQ.0) RETURN
*
* And when alpha.eq.zero.
*
IF (ALPHA.EQ.ZERO) THEN
DO 20 J = 1,N
DO 10 I = 1,M
B(I,J) = ZERO
10 CONTINUE
20 CONTINUE
RETURN
END IF
*
* Start the operations.
*
IF (LSIDE) THEN
IF (LSAME(TRANSA,'N')) THEN
*
* Form B := alpha*A*B.
*
IF (UPPER) THEN
DO 50 J = 1,N
DO 40 K = 1,M
IF (B(K,J).NE.ZERO) THEN
TEMP = ALPHA*B(K,J)
DO 30 I = 1,K - 1
B(I,J) = B(I,J) + TEMP*A(I,K)
30 CONTINUE
IF (NOUNIT) TEMP = TEMP*A(K,K)
B(K,J) = TEMP
END IF
40 CONTINUE
50 CONTINUE
ELSE
DO 80 J = 1,N
DO 70 K = M,1,-1
IF (B(K,J).NE.ZERO) THEN
TEMP = ALPHA*B(K,J)
B(K,J) = TEMP
IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
DO 60 I = K + 1,M
B(I,J) = B(I,J) + TEMP*A(I,K)
60 CONTINUE
END IF
70 CONTINUE
80 CONTINUE
END IF
ELSE
*
* Form B := alpha*A**T*B.
*
IF (UPPER) THEN
DO 110 J = 1,N
DO 100 I = M,1,-1
TEMP = B(I,J)
IF (NOUNIT) TEMP = TEMP*A(I,I)
DO 90 K = 1,I - 1
TEMP = TEMP + A(K,I)*B(K,J)
90 CONTINUE
B(I,J) = ALPHA*TEMP
100 CONTINUE
110 CONTINUE
ELSE
DO 140 J = 1,N
DO 130 I = 1,M
TEMP = B(I,J)
IF (NOUNIT) TEMP = TEMP*A(I,I)
DO 120 K = I + 1,M
TEMP = TEMP + A(K,I)*B(K,J)
120 CONTINUE
B(I,J) = ALPHA*TEMP
130 CONTINUE
140 CONTINUE
END IF
END IF
ELSE
IF (LSAME(TRANSA,'N')) THEN
*
* Form B := alpha*B*A.
*
IF (UPPER) THEN
DO 180 J = N,1,-1
TEMP = ALPHA
IF (NOUNIT) TEMP = TEMP*A(J,J)
DO 150 I = 1,M
B(I,J) = TEMP*B(I,J)
150 CONTINUE
DO 170 K = 1,J - 1
IF (A(K,J).NE.ZERO) THEN
TEMP = ALPHA*A(K,J)
DO 160 I = 1,M
B(I,J) = B(I,J) + TEMP*B(I,K)
160 CONTINUE
END IF
170 CONTINUE
180 CONTINUE
ELSE
DO 220 J = 1,N
TEMP = ALPHA
IF (NOUNIT) TEMP = TEMP*A(J,J)
DO 190 I = 1,M
B(I,J) = TEMP*B(I,J)
190 CONTINUE
DO 210 K = J + 1,N
IF (A(K,J).NE.ZERO) THEN
TEMP = ALPHA*A(K,J)
DO 200 I = 1,M
B(I,J) = B(I,J) + TEMP*B(I,K)
200 CONTINUE
END IF
210 CONTINUE
220 CONTINUE
END IF
ELSE
*
* Form B := alpha*B*A**T.
*
IF (UPPER) THEN
DO 260 K = 1,N
DO 240 J = 1,K - 1
IF (A(J,K).NE.ZERO) THEN
TEMP = ALPHA*A(J,K)
DO 230 I = 1,M
B(I,J) = B(I,J) + TEMP*B(I,K)
230 CONTINUE
END IF
240 CONTINUE
TEMP = ALPHA
IF (NOUNIT) TEMP = TEMP*A(K,K)
IF (TEMP.NE.ONE) THEN
DO 250 I = 1,M
B(I,K) = TEMP*B(I,K)
250 CONTINUE
END IF
260 CONTINUE
ELSE
DO 300 K = N,1,-1
DO 280 J = K + 1,N
IF (A(J,K).NE.ZERO) THEN
TEMP = ALPHA*A(J,K)
DO 270 I = 1,M
B(I,J) = B(I,J) + TEMP*B(I,K)
270 CONTINUE
END IF
280 CONTINUE
TEMP = ALPHA
IF (NOUNIT) TEMP = TEMP*A(K,K)
IF (TEMP.NE.ONE) THEN
DO 290 I = 1,M
B(I,K) = TEMP*B(I,K)
290 CONTINUE
END IF
300 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTRMM
*
END
*> \brief \b DTRMV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,LDA,N
* CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTRMV performs one of the matrix-vector operations
*>
*> x := A*x, or x := A**T*x,
*>
*> where x is an n element vector and A is an n by n unit, or non-unit,
*> upper or lower triangular matrix.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x.
*>
*> TRANS = 'C' or 'c' x := A**T*x.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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
*> transformed vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup trmv
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 2 Blas routine.
*> The vector and matrix arguments are not referenced when N = 0, or M = 0
*>
*> -- 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.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,LDA,N
CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (.NOT.LSAME(TRANS,'N') .AND.
+ .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 2
ELSE IF (.NOT.LSAME(DIAG,'U') .AND.
+ .NOT.LSAME(DIAG,'N')) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (LDA.LT.MAX(1,N)) THEN
INFO = 6
ELSE IF (INCX.EQ.0) THEN
INFO = 8
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTRMV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (N.EQ.0) RETURN
*
NOUNIT = LSAME(DIAG,'N')
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form x := A*x.
*
IF (LSAME(UPLO,'U')) THEN
IF (INCX.EQ.1) THEN
DO 20 J = 1,N
IF (X(J).NE.ZERO) THEN
TEMP = X(J)
DO 10 I = 1,J - 1
X(I) = X(I) + TEMP*A(I,J)
10 CONTINUE
IF (NOUNIT) X(J) = X(J)*A(J,J)
END IF
20 CONTINUE
ELSE
JX = KX
DO 40 J = 1,N
IF (X(JX).NE.ZERO) THEN
TEMP = X(JX)
IX = KX
DO 30 I = 1,J - 1
X(IX) = X(IX) + TEMP*A(I,J)
IX = IX + INCX
30 CONTINUE
IF (NOUNIT) X(JX) = X(JX)*A(J,J)
END IF
JX = JX + INCX
40 CONTINUE
END IF
ELSE
IF (INCX.EQ.1) THEN
DO 60 J = N,1,-1
IF (X(J).NE.ZERO) THEN
TEMP = X(J)
DO 50 I = N,J + 1,-1
X(I) = X(I) + TEMP*A(I,J)
50 CONTINUE
IF (NOUNIT) X(J) = X(J)*A(J,J)
END IF
60 CONTINUE
ELSE
KX = KX + (N-1)*INCX
JX = KX
DO 80 J = N,1,-1
IF (X(JX).NE.ZERO) THEN
TEMP = X(JX)
IX = KX
DO 70 I = N,J + 1,-1
X(IX) = X(IX) + TEMP*A(I,J)
IX = IX - INCX
70 CONTINUE
IF (NOUNIT) X(JX) = X(JX)*A(J,J)
END IF
JX = JX - INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := A**T*x.
*
IF (LSAME(UPLO,'U')) THEN
IF (INCX.EQ.1) THEN
DO 100 J = N,1,-1
TEMP = X(J)
IF (NOUNIT) TEMP = TEMP*A(J,J)
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) TEMP = TEMP*A(J,J)
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.EQ.1) THEN
DO 140 J = 1,N
TEMP = X(J)
IF (NOUNIT) TEMP = TEMP*A(J,J)
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) TEMP = TEMP*A(J,J)
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 of DTRMV
*
END
*> \brief \b DTRSM
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*
* .. Scalar Arguments ..
* DOUBLE PRECISION ALPHA
* INTEGER LDA,LDB,M,N
* CHARACTER DIAG,SIDE,TRANSA,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),B(LDB,*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTRSM solves one of the matrix equations
*>
*> op( A )*X = alpha*B, or X*op( A ) = alpha*B,
*>
*> where alpha is a scalar, X and B are m by n matrices, A is a unit, or
*> non-unit, upper or lower triangular matrix and op( A ) is one of
*>
*> op( A ) = A or op( A ) = A**T.
*>
*> The matrix X is overwritten on B.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] SIDE
*> \verbatim
*> SIDE is CHARACTER*1
*> On entry, SIDE specifies whether op( A ) appears on the left
*> or right of X as follows:
*>
*> SIDE = 'L' or 'l' op( A )*X = alpha*B.
*>
*> SIDE = 'R' or 'r' X*op( A ) = alpha*B.
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*> UPLO is CHARACTER*1
*> On entry, UPLO specifies whether the matrix A is an upper or
*> lower triangular matrix as follows:
*>
*> UPLO = 'U' or 'u' A is an upper triangular matrix.
*>
*> UPLO = 'L' or 'l' A is a lower triangular matrix.
*> \endverbatim
*>
*> \param[in] TRANSA
*> \verbatim
*> TRANSA is CHARACTER*1
*> On entry, TRANSA specifies the form of op( A ) to be used in
*> the matrix multiplication as follows:
*>
*> TRANSA = 'N' or 'n' op( A ) = A.
*>
*> TRANSA = 'T' or 't' op( A ) = A**T.
*>
*> TRANSA = 'C' or 'c' op( A ) = A**T.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*> M is INTEGER
*> On entry, M specifies the number of rows of B. M must be at
*> least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the number of columns of B. N must be
*> at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*> ALPHA is DOUBLE PRECISION.
*> On entry, ALPHA specifies the scalar alpha. When alpha is
*> zero then A is not referenced and B need not be set before
*> entry.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension ( LDA, k ),
*> where k is m when SIDE = 'L' or 'l'
*> and k is n when SIDE = 'R' or 'r'.
*> Before entry with UPLO = 'U' or 'u', the leading k by k
*> upper triangular part of the array A must contain the upper
*> triangular matrix and the strictly lower triangular part of
*> A is not referenced.
*> Before entry with UPLO = 'L' or 'l', the leading k by k
*> lower triangular part of the array A must contain the lower
*> triangular matrix and the strictly upper triangular part of
*> A is not referenced.
*> Note that when DIAG = 'U' or 'u', the diagonal elements of
*> A are not referenced either, but are assumed to be unity.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> On entry, LDA specifies the first dimension of A as declared
*> in the calling (sub) program. When SIDE = 'L' or 'l' then
*> LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
*> then LDA must be at least max( 1, n ).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension ( LDB, N )
*> Before entry, the leading m by n part of the array B must
*> contain the right-hand side matrix B, and on exit is
*> overwritten by the solution matrix X.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> On entry, LDB specifies the first dimension of B as declared
*> in the calling (sub) program. LDB must be at least
*> max( 1, m ).
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup trsm
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> Level 3 Blas routine.
*>
*>
*> -- Written on 8-February-1989.
*> Jack Dongarra, Argonne National Laboratory.
*> Iain Duff, AERE Harwell.
*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
*> Sven Hammarling, Numerical Algorithms Group Ltd.
*> \endverbatim
*>
* =====================================================================
SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*
* -- Reference BLAS level3 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
DOUBLE PRECISION ALPHA
INTEGER LDA,LDB,M,N
CHARACTER DIAG,SIDE,TRANSA,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),B(LDB,*)
* ..
*
* =====================================================================
*
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Local Scalars ..
DOUBLE PRECISION TEMP
INTEGER I,INFO,J,K,NROWA
LOGICAL LSIDE,NOUNIT,UPPER
* ..
* .. Parameters ..
DOUBLE PRECISION ONE,ZERO
PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
* ..
*
* Test the input parameters.
*
LSIDE = LSAME(SIDE,'L')
IF (LSIDE) THEN
NROWA = M
ELSE
NROWA = N
END IF
NOUNIT = LSAME(DIAG,'N')
UPPER = LSAME(UPLO,'U')
*
INFO = 0
IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
INFO = 1
ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
INFO = 2
ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
+ (.NOT.LSAME(TRANSA,'T')) .AND.
+ (.NOT.LSAME(TRANSA,'C'))) THEN
INFO = 3
ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND.
+ (.NOT.LSAME(DIAG,'N'))) THEN
INFO = 4
ELSE IF (M.LT.0) THEN
INFO = 5
ELSE IF (N.LT.0) THEN
INFO = 6
ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
INFO = 9
ELSE IF (LDB.LT.MAX(1,M)) THEN
INFO = 11
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTRSM ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (M.EQ.0 .OR. N.EQ.0) RETURN
*
* And when alpha.eq.zero.
*
IF (ALPHA.EQ.ZERO) THEN
DO 20 J = 1,N
DO 10 I = 1,M
B(I,J) = ZERO
10 CONTINUE
20 CONTINUE
RETURN
END IF
*
* Start the operations.
*
IF (LSIDE) THEN
IF (LSAME(TRANSA,'N')) THEN
*
* Form B := alpha*inv( A )*B.
*
IF (UPPER) THEN
DO 60 J = 1,N
IF (ALPHA.NE.ONE) THEN
DO 30 I = 1,M
B(I,J) = ALPHA*B(I,J)
30 CONTINUE
END IF
DO 50 K = M,1,-1
IF (B(K,J).NE.ZERO) THEN
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
DO 40 I = 1,K - 1
B(I,J) = B(I,J) - B(K,J)*A(I,K)
40 CONTINUE
END IF
50 CONTINUE
60 CONTINUE
ELSE
DO 100 J = 1,N
IF (ALPHA.NE.ONE) THEN
DO 70 I = 1,M
B(I,J) = ALPHA*B(I,J)
70 CONTINUE
END IF
DO 90 K = 1,M
IF (B(K,J).NE.ZERO) THEN
IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
DO 80 I = K + 1,M
B(I,J) = B(I,J) - B(K,J)*A(I,K)
80 CONTINUE
END IF
90 CONTINUE
100 CONTINUE
END IF
ELSE
*
* Form B := alpha*inv( A**T )*B.
*
IF (UPPER) THEN
DO 130 J = 1,N
DO 120 I = 1,M
TEMP = ALPHA*B(I,J)
DO 110 K = 1,I - 1
TEMP = TEMP - A(K,I)*B(K,J)
110 CONTINUE
IF (NOUNIT) TEMP = TEMP/A(I,I)
B(I,J) = TEMP
120 CONTINUE
130 CONTINUE
ELSE
DO 160 J = 1,N
DO 150 I = M,1,-1
TEMP = ALPHA*B(I,J)
DO 140 K = I + 1,M
TEMP = TEMP - A(K,I)*B(K,J)
140 CONTINUE
IF (NOUNIT) TEMP = TEMP/A(I,I)
B(I,J) = TEMP
150 CONTINUE
160 CONTINUE
END IF
END IF
ELSE
IF (LSAME(TRANSA,'N')) THEN
*
* Form B := alpha*B*inv( A ).
*
IF (UPPER) THEN
DO 210 J = 1,N
IF (ALPHA.NE.ONE) THEN
DO 170 I = 1,M
B(I,J) = ALPHA*B(I,J)
170 CONTINUE
END IF
DO 190 K = 1,J - 1
IF (A(K,J).NE.ZERO) THEN
DO 180 I = 1,M
B(I,J) = B(I,J) - A(K,J)*B(I,K)
180 CONTINUE
END IF
190 CONTINUE
IF (NOUNIT) THEN
TEMP = ONE/A(J,J)
DO 200 I = 1,M
B(I,J) = TEMP*B(I,J)
200 CONTINUE
END IF
210 CONTINUE
ELSE
DO 260 J = N,1,-1
IF (ALPHA.NE.ONE) THEN
DO 220 I = 1,M
B(I,J) = ALPHA*B(I,J)
220 CONTINUE
END IF
DO 240 K = J + 1,N
IF (A(K,J).NE.ZERO) THEN
DO 230 I = 1,M
B(I,J) = B(I,J) - A(K,J)*B(I,K)
230 CONTINUE
END IF
240 CONTINUE
IF (NOUNIT) THEN
TEMP = ONE/A(J,J)
DO 250 I = 1,M
B(I,J) = TEMP*B(I,J)
250 CONTINUE
END IF
260 CONTINUE
END IF
ELSE
*
* Form B := alpha*B*inv( A**T ).
*
IF (UPPER) THEN
DO 310 K = N,1,-1
IF (NOUNIT) THEN
TEMP = ONE/A(K,K)
DO 270 I = 1,M
B(I,K) = TEMP*B(I,K)
270 CONTINUE
END IF
DO 290 J = 1,K - 1
IF (A(J,K).NE.ZERO) THEN
TEMP = A(J,K)
DO 280 I = 1,M
B(I,J) = B(I,J) - TEMP*B(I,K)
280 CONTINUE
END IF
290 CONTINUE
IF (ALPHA.NE.ONE) THEN
DO 300 I = 1,M
B(I,K) = ALPHA*B(I,K)
300 CONTINUE
END IF
310 CONTINUE
ELSE
DO 360 K = 1,N
IF (NOUNIT) THEN
TEMP = ONE/A(K,K)
DO 320 I = 1,M
B(I,K) = TEMP*B(I,K)
320 CONTINUE
END IF
DO 340 J = K + 1,N
IF (A(J,K).NE.ZERO) THEN
TEMP = A(J,K)
DO 330 I = 1,M
B(I,J) = B(I,J) - TEMP*B(I,K)
330 CONTINUE
END IF
340 CONTINUE
IF (ALPHA.NE.ONE) THEN
DO 350 I = 1,M
B(I,K) = ALPHA*B(I,K)
350 CONTINUE
END IF
360 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTRSM
*
END
*> \brief \b DTRSV
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,LDA,N
* CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
* DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> DTRSV solves one of the systems of equations
*>
*> A*x = b, or A**T*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.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] UPLO
*> \verbatim
*> UPLO is 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.
*> \endverbatim
*>
*> \param[in] TRANS
*> \verbatim
*> TRANS is 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**T*x = b.
*>
*> TRANS = 'C' or 'c' A**T*x = b.
*> \endverbatim
*>
*> \param[in] DIAG
*> \verbatim
*> DIAG is 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.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> On entry, N specifies the order of the matrix A.
*> N must be at least zero.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is 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 ).
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*> X is DOUBLE PRECISION array, 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.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> On entry, INCX specifies the increment for the elements of
*> X. INCX must not be zero.
*>
*> Level 2 Blas routine.
*>
*> -- Written on 22-October-1986.
*> Jack Dongarra, Argonne National Lab.
*> Jeremy Du Croz, Nag Central Office.
*> Sven Hammarling, Nag Central Office.
*> Richard Hanson, Sandia National Labs.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup trsv
*
* =====================================================================
SUBROUTINE DTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
*
* -- Reference BLAS level2 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,LDA,N
CHARACTER DIAG,TRANS,UPLO
* ..
* .. Array Arguments ..
DOUBLE PRECISION A(LDA,*),X(*)
* ..
*
* =====================================================================
*
* .. 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
* ..
*
* Test the input parameters.
*
INFO = 0
IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
INFO = 1
ELSE IF (.NOT.LSAME(TRANS,'N') .AND.
+ .NOT.LSAME(TRANS,'T') .AND.
+ .NOT.LSAME(TRANS,'C')) THEN
INFO = 2
ELSE IF (.NOT.LSAME(DIAG,'U') .AND.
+ .NOT.LSAME(DIAG,'N')) THEN
INFO = 3
ELSE IF (N.LT.0) THEN
INFO = 4
ELSE IF (LDA.LT.MAX(1,N)) THEN
INFO = 6
ELSE IF (INCX.EQ.0) THEN
INFO = 8
END IF
IF (INFO.NE.0) THEN
CALL XERBLA('DTRSV ',INFO)
RETURN
END IF
*
* Quick return if possible.
*
IF (N.EQ.0) RETURN
*
NOUNIT = LSAME(DIAG,'N')
*
* Set up the start point in X if the increment is not unity. This
* will be ( N - 1 )*INCX too small for descending loops.
*
IF (INCX.LE.0) THEN
KX = 1 - (N-1)*INCX
ELSE IF (INCX.NE.1) THEN
KX = 1
END IF
*
* Start the operations. In this version the elements of A are
* accessed sequentially with one pass through A.
*
IF (LSAME(TRANS,'N')) THEN
*
* Form x := inv( A )*x.
*
IF (LSAME(UPLO,'U')) THEN
IF (INCX.EQ.1) THEN
DO 20 J = N,1,-1
IF (X(J).NE.ZERO) THEN
IF (NOUNIT) X(J) = X(J)/A(J,J)
TEMP = X(J)
DO 10 I = J - 1,1,-1
X(I) = X(I) - TEMP*A(I,J)
10 CONTINUE
END IF
20 CONTINUE
ELSE
JX = KX + (N-1)*INCX
DO 40 J = N,1,-1
IF (X(JX).NE.ZERO) THEN
IF (NOUNIT) X(JX) = X(JX)/A(J,J)
TEMP = X(JX)
IX = JX
DO 30 I = J - 1,1,-1
IX = IX - INCX
X(IX) = X(IX) - TEMP*A(I,J)
30 CONTINUE
END IF
JX = JX - INCX
40 CONTINUE
END IF
ELSE
IF (INCX.EQ.1) THEN
DO 60 J = 1,N
IF (X(J).NE.ZERO) THEN
IF (NOUNIT) X(J) = X(J)/A(J,J)
TEMP = X(J)
DO 50 I = J + 1,N
X(I) = X(I) - TEMP*A(I,J)
50 CONTINUE
END IF
60 CONTINUE
ELSE
JX = KX
DO 80 J = 1,N
IF (X(JX).NE.ZERO) THEN
IF (NOUNIT) X(JX) = X(JX)/A(J,J)
TEMP = X(JX)
IX = JX
DO 70 I = J + 1,N
IX = IX + INCX
X(IX) = X(IX) - TEMP*A(I,J)
70 CONTINUE
END IF
JX = JX + INCX
80 CONTINUE
END IF
END IF
ELSE
*
* Form x := inv( A**T )*x.
*
IF (LSAME(UPLO,'U')) THEN
IF (INCX.EQ.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) TEMP = TEMP/A(J,J)
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) TEMP = TEMP/A(J,J)
X(JX) = TEMP
JX = JX + INCX
120 CONTINUE
END IF
ELSE
IF (INCX.EQ.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) TEMP = TEMP/A(J,J)
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) TEMP = TEMP/A(J,J)
X(JX) = TEMP
JX = JX - INCX
160 CONTINUE
END IF
END IF
END IF
*
RETURN
*
* End of DTRSV
*
END
*> \brief \b IDAMAX
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* INTEGER FUNCTION IDAMAX(N,DX,INCX)
*
* .. Scalar Arguments ..
* INTEGER INCX,N
* ..
* .. Array Arguments ..
* DOUBLE PRECISION DX(*)
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> IDAMAX finds the index of the first element having maximum absolute value.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> number of elements in input vector(s)
*> \endverbatim
*>
*> \param[in] DX
*> \verbatim
*> DX is DOUBLE PRECISION array, dimension ( 1 + ( N - 1 )*abs( INCX ) )
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*> INCX is INTEGER
*> storage spacing between elements of DX
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup iamax
*
*> \par Further Details:
* =====================
*>
*> \verbatim
*>
*> jack dongarra, linpack, 3/11/78.
*> modified 3/93 to return if incx .le. 0.
*> modified 12/3/93, array(1) declarations changed to array(*)
*> \endverbatim
*>
* =====================================================================
INTEGER FUNCTION IDAMAX(N,DX,INCX)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INCX,N
* ..
* .. Array Arguments ..
DOUBLE PRECISION DX(*)
* ..
*
* =====================================================================
*
* .. Local Scalars ..
DOUBLE PRECISION DMAX
INTEGER I,IX
* ..
* .. Intrinsic Functions ..
INTRINSIC DABS
* ..
IDAMAX = 0
IF (N.LT.1 .OR. INCX.LE.0) RETURN
IDAMAX = 1
IF (N.EQ.1) RETURN
IF (INCX.EQ.1) THEN
*
* code for increment equal to 1
*
DMAX = DABS(DX(1))
DO I = 2,N
IF (DABS(DX(I)).GT.DMAX) THEN
IDAMAX = I
DMAX = DABS(DX(I))
END IF
END DO
ELSE
*
* code for increment not equal to 1
*
IX = 1
DMAX = DABS(DX(1))
IX = IX + INCX
DO I = 2,N
IF (DABS(DX(IX)).GT.DMAX) THEN
IDAMAX = I
DMAX = DABS(DX(IX))
END IF
IX = IX + INCX
END DO
END IF
RETURN
*
* End of IDAMAX
*
END
*> \brief \b LSAME
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
* Definition:
* ===========
*
* LOGICAL FUNCTION LSAME(CA,CB)
*
* .. Scalar Arguments ..
* CHARACTER CA,CB
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*> LSAME returns .TRUE. if CA is the same letter as CB regardless of
*> case.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] CA
*> \verbatim
*> CA is CHARACTER*1
*> \endverbatim
*>
*> \param[in] CB
*> \verbatim
*> CB is CHARACTER*1
*> CA and CB specify the single characters to be compared.
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup lsame
*
* =====================================================================
LOGICAL FUNCTION LSAME(CA,CB)
*
* -- Reference BLAS level1 routine --
* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER CA,CB
* ..
*
* =====================================================================
*
* .. Intrinsic Functions ..
INTRINSIC ICHAR
* ..
* .. Local Scalars ..
INTEGER INTA,INTB,ZCODE
* ..
*
* Test if the characters are equal
*
LSAME = CA .EQ. CB
IF (LSAME) RETURN
*
* Now test for equivalence if both characters are alphabetic.
*
ZCODE = ICHAR('Z')
*
* Use 'Z' rather than 'A' so that ASCII can be detected on Prime
* machines, on which ICHAR returns a value with bit 8 set.
* ICHAR('A') on Prime machines returns 193 which is the same as
* ICHAR('A') on an EBCDIC machine.
*
INTA = ICHAR(CA)
INTB = ICHAR(CB)
*
IF (ZCODE.EQ.90 .OR. ZCODE.EQ.122) THEN
*
* ASCII is assumed - ZCODE is the ASCII code of either lower or
* upper case 'Z'.
*
IF (INTA.GE.97 .AND. INTA.LE.122) INTA = INTA - 32
IF (INTB.GE.97 .AND. INTB.LE.122) INTB = INTB - 32
*
ELSE IF (ZCODE.EQ.233 .OR. ZCODE.EQ.169) THEN
*
* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or
* upper case 'Z'.
*
IF (INTA.GE.129 .AND. INTA.LE.137 .OR.
+ INTA.GE.145 .AND. INTA.LE.153 .OR.
+ INTA.GE.162 .AND. INTA.LE.169) INTA = INTA + 64
IF (INTB.GE.129 .AND. INTB.LE.137 .OR.
+ INTB.GE.145 .AND. INTB.LE.153 .OR.
+ INTB.GE.162 .AND. INTB.LE.169) INTB = INTB + 64
*
ELSE IF (ZCODE.EQ.218 .OR. ZCODE.EQ.250) THEN
*
* ASCII is assumed, on Prime machines - ZCODE is the ASCII code
* plus 128 of either lower or upper case 'Z'.
*
IF (INTA.GE.225 .AND. INTA.LE.250) INTA = INTA - 32
IF (INTB.GE.225 .AND. INTB.LE.250) INTB = INTB - 32
END IF
LSAME = INTA .EQ. INTB
*
* RETURN
*
* End of LSAME
*
END