*> \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 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 c IF (Y(JY).NE.ZERO) THEN TEMP = ALPHA*Y(JY) DO 10 I = 1,M A(I,J) = A(I,J) + X(I)*TEMP 10 CONTINUE c END IF JY = JY + INCY 20 CONTINUE ELSE IF (INCX.GT.0) THEN KX = 1 ELSE KX = 1 - (M-1)*INCX END IF DO 40 J = 1,N c IF (Y(JY).NE.ZERO) THEN TEMP = ALPHA*Y(JY) IX = KX DO 30 I = 1,M A(I,J) = A(I,J) + X(IX)*TEMP IX = IX + INCX 30 CONTINUE c END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of 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 c 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 c END IF KK = KK + J 20 CONTINUE ELSE DO 40 J = 1,N c IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*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 c END IF JX = JX + INCX JY = JY + INCY KK = KK + J 40 CONTINUE END IF ELSE * * Form A when lower triangle is stored in AP. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60 J = 1,N c IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN TEMP1 = ALPHA*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 c END IF KK = KK + N - J + 1 60 CONTINUE ELSE DO 80 J = 1,N c IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*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 c 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 c 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 c END IF KK = KK + J 20 CONTINUE ELSE JX = KX DO 40 J = 1,N c IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IX = KX DO 30 K = KK,KK + J - 1 AP(K) = AP(K) + X(IX)*TEMP IX = IX + INCX 30 CONTINUE c END IF JX = JX + INCX KK = KK + J 40 CONTINUE END IF ELSE * * Form A when lower triangle is stored in AP. * IF (INCX.EQ.1) THEN DO 60 J = 1,N c IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) K = KK DO 50 I = J,N AP(K) = AP(K) + X(I)*TEMP K = K + 1 50 CONTINUE c END IF KK = KK + N - J + 1 60 CONTINUE ELSE JX = KX DO 80 J = 1,N c IF (X(JX).NE.ZERO) THEN 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 c 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 double_blas_level2 * *> \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 double_blas_level2 * *> \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 c 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 c END IF 20 CONTINUE ELSE DO 40 J = 1,N c IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*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 c END IF JX = JX + INCX JY = JY + INCY 40 CONTINUE END IF ELSE * * Form A when A is stored in the lower triangle. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60 J = 1,N c IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN TEMP1 = ALPHA*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 c END IF 60 CONTINUE ELSE DO 80 J = 1,N c IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*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 c 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 c 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 c END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1,N IF (BETA.EQ.ZERO) THEN DO 140 I = J,N C(I,J) = ZERO 140 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 150 I = J,N C(I,J) = BETA*C(I,J) 150 CONTINUE END IF DO 170 L = 1,K c IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN TEMP1 = ALPHA*B(J,L) TEMP2 = ALPHA*A(J,L) DO 160 I = J,N C(I,J) = C(I,J) + A(I,L)*TEMP1 + + B(I,L)*TEMP2 160 CONTINUE c END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*A**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 c 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 c END IF 20 CONTINUE ELSE JX = KX DO 40 J = 1,N c 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 c END IF JX = JX + INCX 40 CONTINUE END IF ELSE * * Form A when A is stored in lower triangle. * IF (INCX.EQ.1) THEN DO 60 J = 1,N c IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 50 I = J,N A(I,J) = A(I,J) + X(I)*TEMP 50 CONTINUE c END IF 60 CONTINUE ELSE JX = KX DO 80 J = 1,N c 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 c 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 c IF (A(J,L).NE.ZERO) THEN TEMP = ALPHA*A(J,L) DO 110 I = 1,J C(I,J) = C(I,J) + TEMP*A(I,L) 110 CONTINUE c END IF 120 CONTINUE 130 CONTINUE ELSE DO 180 J = 1,N IF (BETA.EQ.ZERO) THEN DO 140 I = J,N C(I,J) = ZERO 140 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 150 I = J,N C(I,J) = BETA*C(I,J) 150 CONTINUE END IF DO 170 L = 1,K c IF (A(J,L).NE.ZERO) THEN TEMP = ALPHA*A(J,L) DO 160 I = J,N C(I,J) = C(I,J) + TEMP*A(I,L) 160 CONTINUE c END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*A**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 c IF (X(J).NE.ZERO) THEN TEMP = X(J) L = KPLUS1 - J DO 10 I = MAX(1,J-K),J - 1 X(I) = X(I) + TEMP*A(L+I,J) 10 CONTINUE IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J) c END IF 20 CONTINUE ELSE JX = KX DO 40 J = 1,N c IF (X(JX).NE.ZERO) THEN TEMP = X(JX) IX = KX L = KPLUS1 - J DO 30 I = MAX(1,J-K),J - 1 X(IX) = X(IX) + TEMP*A(L+I,J) IX = IX + INCX 30 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J) c END IF JX = JX + INCX IF (J.GT.K) KX = KX + INCX 40 CONTINUE END IF ELSE IF (INCX.EQ.1) THEN DO 60 J = N,1,-1 c IF (X(J).NE.ZERO) THEN TEMP = X(J) L = 1 - J DO 50 I = MIN(N,J+K),J + 1,-1 X(I) = X(I) + TEMP*A(L+I,J) 50 CONTINUE IF (NOUNIT) X(J) = X(J)*A(1,J) c END IF 60 CONTINUE ELSE KX = KX + (N-1)*INCX JX = KX DO 80 J = N,1,-1 c IF (X(JX).NE.ZERO) THEN TEMP = X(JX) IX = KX L = 1 - J DO 70 I = MIN(N,J+K),J + 1,-1 X(IX) = X(IX) + TEMP*A(L+I,J) IX = IX - INCX 70 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(1,J) c END IF JX = JX - INCX IF ((N-J).GE.K) KX = KX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A**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 c IF (X(J).NE.ZERO) THEN L = KPLUS1 - J IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J) TEMP = X(J) DO 10 I = J - 1,MAX(1,J-K),-1 X(I) = X(I) - TEMP*A(L+I,J) 10 CONTINUE c END IF 20 CONTINUE ELSE KX = KX + (N-1)*INCX JX = KX DO 40 J = N,1,-1 KX = KX - INCX c IF (X(JX).NE.ZERO) THEN IX = KX L = KPLUS1 - J IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J) TEMP = X(JX) DO 30 I = J - 1,MAX(1,J-K),-1 X(IX) = X(IX) - TEMP*A(L+I,J) IX = IX - INCX 30 CONTINUE c END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF (INCX.EQ.1) THEN DO 60 J = 1,N c IF (X(J).NE.ZERO) THEN L = 1 - J IF (NOUNIT) X(J) = X(J)/A(1,J) TEMP = X(J) DO 50 I = J + 1,MIN(N,J+K) X(I) = X(I) - TEMP*A(L+I,J) 50 CONTINUE c END IF 60 CONTINUE ELSE JX = KX DO 80 J = 1,N KX = KX + INCX c IF (X(JX).NE.ZERO) THEN IX = KX L = 1 - J IF (NOUNIT) X(JX) = X(JX)/A(1,J) TEMP = X(JX) DO 70 I = J + 1,MIN(N,J+K) X(IX) = X(IX) - TEMP*A(L+I,J) IX = IX + INCX 70 CONTINUE c END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form x := inv( A**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 c IF (X(J).NE.ZERO) THEN TEMP = X(J) K = KK DO 10 I = 1,J - 1 X(I) = X(I) + TEMP*AP(K) K = K + 1 10 CONTINUE IF (NOUNIT) X(J) = X(J)*AP(KK+J-1) c END IF KK = KK + J 20 CONTINUE ELSE JX = KX DO 40 J = 1,N c IF (X(JX).NE.ZERO) THEN TEMP = X(JX) IX = KX DO 30 K = KK,KK + J - 2 X(IX) = X(IX) + TEMP*AP(K) IX = IX + INCX 30 CONTINUE IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1) c END IF JX = JX + INCX KK = KK + J 40 CONTINUE END IF ELSE KK = (N* (N+1))/2 IF (INCX.EQ.1) THEN DO 60 J = N,1,-1 c IF (X(J).NE.ZERO) THEN TEMP = X(J) K = KK DO 50 I = N,J + 1,-1 X(I) = X(I) + TEMP*AP(K) K = K - 1 50 CONTINUE IF (NOUNIT) X(J) = X(J)*AP(KK-N+J) c END IF KK = KK - (N-J+1) 60 CONTINUE ELSE KX = KX + (N-1)*INCX JX = KX DO 80 J = N,1,-1 c IF (X(JX).NE.ZERO) THEN TEMP = X(JX) IX = KX DO 70 K = KK,KK - (N- (J+1)),-1 X(IX) = X(IX) + TEMP*AP(K) IX = IX - INCX 70 CONTINUE IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J) c END IF JX = JX - INCX KK = KK - (N-J+1) 80 CONTINUE END IF END IF ELSE * * Form x := A**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 c IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/AP(KK) TEMP = X(J) K = KK - 1 DO 10 I = J - 1,1,-1 X(I) = X(I) - TEMP*AP(K) K = K - 1 10 CONTINUE c END IF KK = KK - J 20 CONTINUE ELSE JX = KX + (N-1)*INCX DO 40 J = N,1,-1 c IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/AP(KK) TEMP = X(JX) IX = JX DO 30 K = KK - 1,KK - J + 1,-1 IX = IX - INCX X(IX) = X(IX) - TEMP*AP(K) 30 CONTINUE c END IF JX = JX - INCX KK = KK - J 40 CONTINUE END IF ELSE KK = 1 IF (INCX.EQ.1) THEN DO 60 J = 1,N c IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/AP(KK) TEMP = X(J) K = KK + 1 DO 50 I = J + 1,N X(I) = X(I) - TEMP*AP(K) K = K + 1 50 CONTINUE c END IF KK = KK + (N-J+1) 60 CONTINUE ELSE JX = KX DO 80 J = 1,N c IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/AP(KK) TEMP = X(JX) IX = JX DO 70 K = KK + 1,KK + N - J IX = IX + INCX X(IX) = X(IX) - TEMP*AP(K) 70 CONTINUE c END IF JX = JX + INCX KK = KK + (N-J+1) 80 CONTINUE END IF END IF ELSE * * Form x := inv( A**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 c IF (B(K,J).NE.ZERO) THEN TEMP = ALPHA*B(K,J) DO 30 I = 1,K - 1 B(I,J) = B(I,J) + TEMP*A(I,K) 30 CONTINUE IF (NOUNIT) TEMP = TEMP*A(K,K) B(K,J) = TEMP c END IF 40 CONTINUE 50 CONTINUE ELSE DO 80 J = 1,N DO 70 K = M,1,-1 c IF (B(K,J).NE.ZERO) THEN TEMP = ALPHA*B(K,J) B(K,J) = TEMP IF (NOUNIT) B(K,J) = B(K,J)*A(K,K) DO 60 I = K + 1,M B(I,J) = B(I,J) + TEMP*A(I,K) 60 CONTINUE c END IF 70 CONTINUE 80 CONTINUE END IF ELSE * * Form B := alpha*A**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 c 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 c 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 c 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 c 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 c 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 c 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 c 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 c 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 c IF (X(J).NE.ZERO) THEN TEMP = X(J) DO 10 I = 1,J - 1 X(I) = X(I) + TEMP*A(I,J) 10 CONTINUE IF (NOUNIT) X(J) = X(J)*A(J,J) c END IF 20 CONTINUE ELSE JX = KX DO 40 J = 1,N c IF (X(JX).NE.ZERO) THEN TEMP = X(JX) IX = KX DO 30 I = 1,J - 1 X(IX) = X(IX) + TEMP*A(I,J) IX = IX + INCX 30 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(J,J) c END IF JX = JX + INCX 40 CONTINUE END IF ELSE IF (INCX.EQ.1) THEN DO 60 J = N,1,-1 c IF (X(J).NE.ZERO) THEN TEMP = X(J) DO 50 I = N,J + 1,-1 X(I) = X(I) + TEMP*A(I,J) 50 CONTINUE IF (NOUNIT) X(J) = X(J)*A(J,J) c END IF 60 CONTINUE ELSE KX = KX + (N-1)*INCX JX = KX DO 80 J = N,1,-1 c IF (X(JX).NE.ZERO) THEN TEMP = X(JX) IX = KX DO 70 I = N,J + 1,-1 X(IX) = X(IX) + TEMP*A(I,J) IX = IX - INCX 70 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(J,J) c END IF JX = JX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A**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 c IF (B(K,J).NE.ZERO) THEN IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) DO 40 I = 1,K - 1 B(I,J) = B(I,J) - B(K,J)*A(I,K) 40 CONTINUE c END IF 50 CONTINUE 60 CONTINUE ELSE DO 100 J = 1,N IF (ALPHA.NE.ONE) THEN DO 70 I = 1,M B(I,J) = ALPHA*B(I,J) 70 CONTINUE END IF DO 90 K = 1,M c IF (B(K,J).NE.ZERO) THEN IF (NOUNIT) B(K,J) = B(K,J)/A(K,K) DO 80 I = K + 1,M B(I,J) = B(I,J) - B(K,J)*A(I,K) 80 CONTINUE c END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A**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 c 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 c 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 c 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 c 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 c 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 c 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 c 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 c 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 c IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/A(J,J) TEMP = X(J) DO 10 I = J - 1,1,-1 X(I) = X(I) - TEMP*A(I,J) 10 CONTINUE c END IF 20 CONTINUE ELSE JX = KX + (N-1)*INCX DO 40 J = N,1,-1 c IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/A(J,J) TEMP = X(JX) IX = JX DO 30 I = J - 1,1,-1 IX = IX - INCX X(IX) = X(IX) - TEMP*A(I,J) 30 CONTINUE c END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF (INCX.EQ.1) THEN DO 60 J = 1,N c IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/A(J,J) TEMP = X(J) DO 50 I = J + 1,N X(I) = X(I) - TEMP*A(I,J) 50 CONTINUE c END IF 60 CONTINUE ELSE JX = KX DO 80 J = 1,N c IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/A(J,J) TEMP = X(JX) IX = JX DO 70 I = J + 1,N IX = IX + INCX X(IX) = X(IX) - TEMP*A(I,J) 70 CONTINUE c END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form x := inv( A**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