diff --git a/src/64to32blas/GNUmakefile b/src/64to32blas/GNUmakefile
index c470c5f3b1..1bad1835d2 100644
--- a/src/64to32blas/GNUmakefile
+++ b/src/64to32blas/GNUmakefile
@@ -10,6 +10,7 @@
ypotri.o ypotrf.o ysygv.o ygeev.o ygeevx.o \
ygels.o ygelss.o ygtsv.o\
ypftrf.o ypftri.o ytfsm.o \
+ ygetrf.o ygetri.o \
ifily.o\
jscal.o jaxpy.o jgemm.o jheev.o jcopy.o jdotc.o \
ijamax.o jgesvd.o jgeev.o jgetrs.o jgetrf.o \
diff --git a/src/64to32blas/ygetrf.F b/src/64to32blas/ygetrf.F
new file mode 100644
index 0000000000..d642077a1b
--- /dev/null
+++ b/src/64to32blas/ygetrf.F
@@ -0,0 +1,20 @@
+ SUBROUTINE ygetrf( M, N, A, LDA, ipiv, INFO )
+*
+* $id$
+*
+ implicit none
+#include "y64.fh"
+ integer ipiv(*)
+ integer info, lda, m, n
+ double precision a( lda, * )
+c
+ INTGR4 m4,n4,info4,lda4
+c
+ n4=n
+ m4=n
+ lda4=lda
+c
+ call dgetrf( m4, n4, a, lda4, ipiv, info4 )
+ info=info4
+ return
+ end
diff --git a/src/64to32blas/ygetri.F b/src/64to32blas/ygetri.F
new file mode 100644
index 0000000000..cc0392dc97
--- /dev/null
+++ b/src/64to32blas/ygetri.F
@@ -0,0 +1,21 @@
+ SUBROUTINE ygetri( N, A, LDA, ipiv, WORK, LWORK, INFO )
+*
+* $id$
+*
+ implicit none
+#include "y64.fh"
+ integer ipiv(*)
+ integer info, lda, n, lwork
+ double precision a( lda, * )
+ double precision work(*)
+c
+ INTGR4 n4,info4,lda4,lwork4
+c
+ n4=n
+ lda4=lda
+ lwork4=lwork
+c
+ call dgetri( n4, a, lda4, ipiv, work, lwork4, info4 )
+ info=info4
+ return
+ end
diff --git a/src/lapack/GNUmakefile b/src/lapack/GNUmakefile
index a9aca3fb4e..776ea9d556 100644
--- a/src/lapack/GNUmakefile
+++ b/src/lapack/GNUmakefile
@@ -84,7 +84,7 @@ OBJ_OPTIMIZE = dbdsqr.o dgebd2.o dgebrd.o dgelq2.o \
dlanhs.o dlarfx.o dlahrd.o dlanv2.o \
dlarra.o dlarrb.o dlarrc.o dlarrd.o dlarrf.o dlarrk.o dlarrv.o\
dlar1v.o dlaneg.o dgelss.o dormtr.o dormql.o dorm2l.o \
- dsfrk.o dsytrf.o \
+ dsfrk.o dsytrf.o dgetri.o dtrsm.o dtrmv.o dtrmm.o\
dpftrf.o dpftri.o dtfsm.o dtftri.o zlacp2.o \
zsysv.o zsytrf.o zsytf2.o zsytrs.o zsytrs2.o zsyr.o \
zsyconv.o zlasyf.o
diff --git a/src/lapack/double/dgetri.f b/src/lapack/double/dgetri.f
new file mode 100644
index 0000000000..92ef90c186
--- /dev/null
+++ b/src/lapack/double/dgetri.f
@@ -0,0 +1,258 @@
+*> \brief \b DGETRI
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DGETRI + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
+*
+* .. Scalar Arguments ..
+* INTEGER INFO, LDA, LWORK, N
+* ..
+* .. Array Arguments ..
+* INTEGER IPIV( * )
+* DOUBLE PRECISION A( LDA, * ), WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DGETRI computes the inverse of a matrix using the LU factorization
+*> computed by DGETRF.
+*>
+*> This method inverts U and then computes inv(A) by solving the system
+*> inv(A)*L = inv(U) for inv(A).
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is DOUBLE PRECISION array, dimension (LDA,N)
+*> On entry, the factors L and U from the factorization
+*> A = P*L*U as computed by DGETRF.
+*> On exit, if INFO = 0, the inverse of the original matrix A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*> IPIV is INTEGER array, dimension (N)
+*> The pivot indices from DGETRF; for 1<=i<=N, row i of the
+*> matrix was interchanged with row IPIV(i).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
+*> On exit, if INFO=0, then WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The dimension of the array WORK. LWORK >= max(1,N).
+*> For optimal performance LWORK >= N*NB, where NB is
+*> the optimal blocksize returned by ILAENV.
+*>
+*> If LWORK = -1, then a workspace query is assumed; the routine
+*> only calculates the optimal size of the WORK array, returns
+*> this value as the first entry of the WORK array, and no error
+*> message related to LWORK is issued by XERBLA.
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
+*> singular and its inverse could not be computed.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \ingroup doubleGEcomputational
+*
+* =====================================================================
+ SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
+*
+* -- LAPACK computational routine --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+* .. Scalar Arguments ..
+ INTEGER INFO, LDA, LWORK, N
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ DOUBLE PRECISION A( LDA, * ), WORK( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL LQUERY
+ INTEGER I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB,
+ $ NBMIN, NN
+* ..
+* .. External Functions ..
+ INTEGER ILAENV
+ EXTERNAL ILAENV
+* ..
+* .. External Subroutines ..
+ EXTERNAL DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 )
+ LWKOPT = N*NB
+ WORK( 1 ) = LWKOPT
+ LQUERY = ( LWORK.EQ.-1 )
+ IF( N.LT.0 ) THEN
+ INFO = -1
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -3
+ ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
+ INFO = -6
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DGETRI', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Form inv(U). If INFO > 0 from DTRTRI, then U is singular,
+* and the inverse is not computed.
+*
+ CALL DTRTRI( 'Upper', 'Non-unit', N, A, LDA, INFO )
+ IF( INFO.GT.0 )
+ $ RETURN
+*
+ NBMIN = 2
+ LDWORK = N
+ IF( NB.GT.1 .AND. NB.LT.N ) THEN
+ IWS = MAX( LDWORK*NB, 1 )
+ IF( LWORK.LT.IWS ) THEN
+ NB = LWORK / LDWORK
+ NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) )
+ END IF
+ ELSE
+ IWS = N
+ END IF
+*
+* Solve the equation inv(A)*L = inv(U) for inv(A).
+*
+ IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
+*
+* Use unblocked code.
+*
+ DO 20 J = N, 1, -1
+*
+* Copy current column of L to WORK and replace with zeros.
+*
+ DO 10 I = J + 1, N
+ WORK( I ) = A( I, J )
+ A( I, J ) = ZERO
+ 10 CONTINUE
+*
+* Compute current column of inv(A).
+*
+ IF( J.LT.N )
+ $ CALL DGEMV( 'No transpose', N, N-J, -ONE, A( 1, J+1 ),
+ $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
+ 20 CONTINUE
+ ELSE
+*
+* Use blocked code.
+*
+ NN = ( ( N-1 ) / NB )*NB + 1
+ DO 50 J = NN, 1, -NB
+ JB = MIN( NB, N-J+1 )
+*
+* Copy current block column of L to WORK and replace with
+* zeros.
+*
+ DO 40 JJ = J, J + JB - 1
+ DO 30 I = JJ + 1, N
+ WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
+ A( I, JJ ) = ZERO
+ 30 CONTINUE
+ 40 CONTINUE
+*
+* Compute current block column of inv(A).
+*
+ IF( J+JB.LE.N )
+ $ CALL DGEMM( 'No transpose', 'No transpose', N, JB,
+ $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
+ $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
+ CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', N, JB,
+ $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
+ 50 CONTINUE
+ END IF
+*
+* Apply column interchanges.
+*
+ DO 60 J = N - 1, 1, -1
+ JP = IPIV( J )
+ IF( JP.NE.J )
+ $ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
+ 60 CONTINUE
+*
+ WORK( 1 ) = IWS
+ RETURN
+*
+* End of DGETRI
+*
+ END
diff --git a/src/lapack/double/dtrmm.f b/src/lapack/double/dtrmm.f
new file mode 100644
index 0000000000..b2cc0a1fa8
--- /dev/null
+++ b/src/lapack/double/dtrmm.f
@@ -0,0 +1,412 @@
+*> \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 double_blas_level3
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> Level 3 Blas routine.
+*>
+*> -- Written on 8-February-1989.
+*> Jack Dongarra, Argonne National Laboratory.
+*> Iain Duff, AERE Harwell.
+*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
+*> Sven Hammarling, Numerical Algorithms Group Ltd.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+*
+* -- Reference BLAS level3 routine --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+* .. Scalar Arguments ..
+ DOUBLE PRECISION ALPHA
+ INTEGER LDA,LDB,M,N
+ CHARACTER DIAG,SIDE,TRANSA,UPLO
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*),B(LDB,*)
+* ..
+*
+* =====================================================================
+*
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX
+* ..
+* .. Local Scalars ..
+ DOUBLE PRECISION TEMP
+ INTEGER I,INFO,J,K,NROWA
+ LOGICAL LSIDE,NOUNIT,UPPER
+* ..
+* .. Parameters ..
+ DOUBLE PRECISION ONE,ZERO
+ PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
+* ..
+*
+* Test the input parameters.
+*
+ LSIDE = LSAME(SIDE,'L')
+ IF (LSIDE) THEN
+ NROWA = M
+ ELSE
+ NROWA = N
+ END IF
+ NOUNIT = LSAME(DIAG,'N')
+ UPPER = LSAME(UPLO,'U')
+*
+ INFO = 0
+ IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
+ INFO = 1
+ ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
+ INFO = 2
+ ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
+ + (.NOT.LSAME(TRANSA,'T')) .AND.
+ + (.NOT.LSAME(TRANSA,'C'))) THEN
+ INFO = 3
+ ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
+ INFO = 4
+ ELSE IF (M.LT.0) THEN
+ INFO = 5
+ ELSE IF (N.LT.0) THEN
+ INFO = 6
+ ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
+ INFO = 9
+ ELSE IF (LDB.LT.MAX(1,M)) THEN
+ INFO = 11
+ END IF
+ IF (INFO.NE.0) THEN
+ CALL XERBLA('DTRMM ',INFO)
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF (M.EQ.0 .OR. N.EQ.0) RETURN
+*
+* And when alpha.eq.zero.
+*
+ IF (ALPHA.EQ.ZERO) THEN
+ DO 20 J = 1,N
+ DO 10 I = 1,M
+ B(I,J) = ZERO
+ 10 CONTINUE
+ 20 CONTINUE
+ RETURN
+ END IF
+*
+* Start the operations.
+*
+ IF (LSIDE) THEN
+ IF (LSAME(TRANSA,'N')) THEN
+*
+* Form B := alpha*A*B.
+*
+ IF (UPPER) THEN
+ DO 50 J = 1,N
+ DO 40 K = 1,M
+ IF (B(K,J).NE.ZERO) THEN
+ TEMP = ALPHA*B(K,J)
+ DO 30 I = 1,K - 1
+ B(I,J) = B(I,J) + TEMP*A(I,K)
+ 30 CONTINUE
+ IF (NOUNIT) TEMP = TEMP*A(K,K)
+ B(K,J) = TEMP
+ END IF
+ 40 CONTINUE
+ 50 CONTINUE
+ ELSE
+ DO 80 J = 1,N
+ DO 70 K = M,1,-1
+ IF (B(K,J).NE.ZERO) THEN
+ TEMP = ALPHA*B(K,J)
+ B(K,J) = TEMP
+ IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
+ DO 60 I = K + 1,M
+ B(I,J) = B(I,J) + TEMP*A(I,K)
+ 60 CONTINUE
+ END IF
+ 70 CONTINUE
+ 80 CONTINUE
+ END IF
+ ELSE
+*
+* Form B := alpha*A**T*B.
+*
+ IF (UPPER) THEN
+ DO 110 J = 1,N
+ DO 100 I = M,1,-1
+ TEMP = B(I,J)
+ IF (NOUNIT) TEMP = TEMP*A(I,I)
+ DO 90 K = 1,I - 1
+ TEMP = TEMP + A(K,I)*B(K,J)
+ 90 CONTINUE
+ B(I,J) = ALPHA*TEMP
+ 100 CONTINUE
+ 110 CONTINUE
+ ELSE
+ DO 140 J = 1,N
+ DO 130 I = 1,M
+ TEMP = B(I,J)
+ IF (NOUNIT) TEMP = TEMP*A(I,I)
+ DO 120 K = I + 1,M
+ TEMP = TEMP + A(K,I)*B(K,J)
+ 120 CONTINUE
+ B(I,J) = ALPHA*TEMP
+ 130 CONTINUE
+ 140 CONTINUE
+ END IF
+ END IF
+ ELSE
+ IF (LSAME(TRANSA,'N')) THEN
+*
+* Form B := alpha*B*A.
+*
+ IF (UPPER) THEN
+ DO 180 J = N,1,-1
+ TEMP = ALPHA
+ IF (NOUNIT) TEMP = TEMP*A(J,J)
+ DO 150 I = 1,M
+ B(I,J) = TEMP*B(I,J)
+ 150 CONTINUE
+ DO 170 K = 1,J - 1
+ IF (A(K,J).NE.ZERO) THEN
+ TEMP = ALPHA*A(K,J)
+ DO 160 I = 1,M
+ B(I,J) = B(I,J) + TEMP*B(I,K)
+ 160 CONTINUE
+ END IF
+ 170 CONTINUE
+ 180 CONTINUE
+ ELSE
+ DO 220 J = 1,N
+ TEMP = ALPHA
+ IF (NOUNIT) TEMP = TEMP*A(J,J)
+ DO 190 I = 1,M
+ B(I,J) = TEMP*B(I,J)
+ 190 CONTINUE
+ DO 210 K = J + 1,N
+ IF (A(K,J).NE.ZERO) THEN
+ TEMP = ALPHA*A(K,J)
+ DO 200 I = 1,M
+ B(I,J) = B(I,J) + TEMP*B(I,K)
+ 200 CONTINUE
+ END IF
+ 210 CONTINUE
+ 220 CONTINUE
+ END IF
+ ELSE
+*
+* Form B := alpha*B*A**T.
+*
+ IF (UPPER) THEN
+ DO 260 K = 1,N
+ DO 240 J = 1,K - 1
+ IF (A(J,K).NE.ZERO) THEN
+ TEMP = ALPHA*A(J,K)
+ DO 230 I = 1,M
+ B(I,J) = B(I,J) + TEMP*B(I,K)
+ 230 CONTINUE
+ END IF
+ 240 CONTINUE
+ TEMP = ALPHA
+ IF (NOUNIT) TEMP = TEMP*A(K,K)
+ IF (TEMP.NE.ONE) THEN
+ DO 250 I = 1,M
+ B(I,K) = TEMP*B(I,K)
+ 250 CONTINUE
+ END IF
+ 260 CONTINUE
+ ELSE
+ DO 300 K = N,1,-1
+ DO 280 J = K + 1,N
+ IF (A(J,K).NE.ZERO) THEN
+ TEMP = ALPHA*A(J,K)
+ DO 270 I = 1,M
+ B(I,J) = B(I,J) + TEMP*B(I,K)
+ 270 CONTINUE
+ END IF
+ 280 CONTINUE
+ TEMP = ALPHA
+ IF (NOUNIT) TEMP = TEMP*A(K,K)
+ IF (TEMP.NE.ONE) THEN
+ DO 290 I = 1,M
+ B(I,K) = TEMP*B(I,K)
+ 290 CONTINUE
+ END IF
+ 300 CONTINUE
+ END IF
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of DTRMM
+*
+ END
diff --git a/src/lapack/double/dtrmv.f b/src/lapack/double/dtrmv.f
new file mode 100644
index 0000000000..e8af8e6136
--- /dev/null
+++ b/src/lapack/double/dtrmv.f
@@ -0,0 +1,339 @@
+*> \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 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 DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+*
+* -- Reference BLAS level2 routine --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+* .. Scalar Arguments ..
+ INTEGER INCX,LDA,N
+ CHARACTER DIAG,TRANS,UPLO
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*),X(*)
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO
+ PARAMETER (ZERO=0.0D+0)
+* ..
+* .. Local Scalars ..
+ DOUBLE PRECISION TEMP
+ INTEGER I,INFO,IX,J,JX,KX
+ LOGICAL NOUNIT
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX
+* ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
+ INFO = 1
+ ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
+ + .NOT.LSAME(TRANS,'C')) THEN
+ INFO = 2
+ ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
+ INFO = 3
+ ELSE IF (N.LT.0) THEN
+ INFO = 4
+ ELSE IF (LDA.LT.MAX(1,N)) THEN
+ INFO = 6
+ ELSE IF (INCX.EQ.0) THEN
+ INFO = 8
+ END IF
+ IF (INFO.NE.0) THEN
+ CALL XERBLA('DTRMV ',INFO)
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF (N.EQ.0) RETURN
+*
+ NOUNIT = LSAME(DIAG,'N')
+*
+* Set up the start point in X if the increment is not unity. This
+* will be ( N - 1 )*INCX too small for descending loops.
+*
+ IF (INCX.LE.0) THEN
+ KX = 1 - (N-1)*INCX
+ ELSE IF (INCX.NE.1) THEN
+ KX = 1
+ END IF
+*
+* Start the operations. In this version the elements of A are
+* accessed sequentially with one pass through A.
+*
+ IF (LSAME(TRANS,'N')) THEN
+*
+* Form x := A*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ IF (INCX.EQ.1) THEN
+ DO 20 J = 1,N
+ IF (X(J).NE.ZERO) THEN
+ TEMP = X(J)
+ DO 10 I = 1,J - 1
+ X(I) = X(I) + TEMP*A(I,J)
+ 10 CONTINUE
+ IF (NOUNIT) X(J) = X(J)*A(J,J)
+ END IF
+ 20 CONTINUE
+ ELSE
+ JX = KX
+ DO 40 J = 1,N
+ IF (X(JX).NE.ZERO) THEN
+ TEMP = X(JX)
+ IX = KX
+ DO 30 I = 1,J - 1
+ X(IX) = X(IX) + TEMP*A(I,J)
+ IX = IX + INCX
+ 30 CONTINUE
+ IF (NOUNIT) X(JX) = X(JX)*A(J,J)
+ END IF
+ JX = JX + INCX
+ 40 CONTINUE
+ END IF
+ ELSE
+ IF (INCX.EQ.1) THEN
+ DO 60 J = N,1,-1
+ IF (X(J).NE.ZERO) THEN
+ TEMP = X(J)
+ DO 50 I = N,J + 1,-1
+ X(I) = X(I) + TEMP*A(I,J)
+ 50 CONTINUE
+ IF (NOUNIT) X(J) = X(J)*A(J,J)
+ END IF
+ 60 CONTINUE
+ ELSE
+ KX = KX + (N-1)*INCX
+ JX = KX
+ DO 80 J = N,1,-1
+ IF (X(JX).NE.ZERO) THEN
+ TEMP = X(JX)
+ IX = KX
+ DO 70 I = N,J + 1,-1
+ X(IX) = X(IX) + TEMP*A(I,J)
+ IX = IX - INCX
+ 70 CONTINUE
+ IF (NOUNIT) X(JX) = X(JX)*A(J,J)
+ END IF
+ JX = JX - INCX
+ 80 CONTINUE
+ END IF
+ END IF
+ ELSE
+*
+* Form x := A**T*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ IF (INCX.EQ.1) THEN
+ DO 100 J = N,1,-1
+ TEMP = X(J)
+ IF (NOUNIT) TEMP = TEMP*A(J,J)
+ DO 90 I = J - 1,1,-1
+ TEMP = TEMP + A(I,J)*X(I)
+ 90 CONTINUE
+ X(J) = TEMP
+ 100 CONTINUE
+ ELSE
+ JX = KX + (N-1)*INCX
+ DO 120 J = N,1,-1
+ TEMP = X(JX)
+ IX = JX
+ IF (NOUNIT) TEMP = TEMP*A(J,J)
+ DO 110 I = J - 1,1,-1
+ IX = IX - INCX
+ TEMP = TEMP + A(I,J)*X(IX)
+ 110 CONTINUE
+ X(JX) = TEMP
+ JX = JX - INCX
+ 120 CONTINUE
+ END IF
+ ELSE
+ IF (INCX.EQ.1) THEN
+ DO 140 J = 1,N
+ TEMP = X(J)
+ IF (NOUNIT) TEMP = TEMP*A(J,J)
+ DO 130 I = J + 1,N
+ TEMP = TEMP + A(I,J)*X(I)
+ 130 CONTINUE
+ X(J) = TEMP
+ 140 CONTINUE
+ ELSE
+ JX = KX
+ DO 160 J = 1,N
+ TEMP = X(JX)
+ IX = JX
+ IF (NOUNIT) TEMP = TEMP*A(J,J)
+ DO 150 I = J + 1,N
+ IX = IX + INCX
+ TEMP = TEMP + A(I,J)*X(IX)
+ 150 CONTINUE
+ X(JX) = TEMP
+ JX = JX + INCX
+ 160 CONTINUE
+ END IF
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of DTRMV
+*
+ END
diff --git a/src/lapack/double/dtrsm.f b/src/lapack/double/dtrsm.f
new file mode 100644
index 0000000000..fa8080bc92
--- /dev/null
+++ b/src/lapack/double/dtrsm.f
@@ -0,0 +1,440 @@
+*> \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 double_blas_level3
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> Level 3 Blas routine.
+*>
+*>
+*> -- Written on 8-February-1989.
+*> Jack Dongarra, Argonne National Laboratory.
+*> Iain Duff, AERE Harwell.
+*> Jeremy Du Croz, Numerical Algorithms Group Ltd.
+*> Sven Hammarling, Numerical Algorithms Group Ltd.
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+*
+* -- Reference BLAS level3 routine --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+* .. Scalar Arguments ..
+ DOUBLE PRECISION ALPHA
+ INTEGER LDA,LDB,M,N
+ CHARACTER DIAG,SIDE,TRANSA,UPLO
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*),B(LDB,*)
+* ..
+*
+* =====================================================================
+*
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX
+* ..
+* .. Local Scalars ..
+ DOUBLE PRECISION TEMP
+ INTEGER I,INFO,J,K,NROWA
+ LOGICAL LSIDE,NOUNIT,UPPER
+* ..
+* .. Parameters ..
+ DOUBLE PRECISION ONE,ZERO
+ PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
+* ..
+*
+* Test the input parameters.
+*
+ LSIDE = LSAME(SIDE,'L')
+ IF (LSIDE) THEN
+ NROWA = M
+ ELSE
+ NROWA = N
+ END IF
+ NOUNIT = LSAME(DIAG,'N')
+ UPPER = LSAME(UPLO,'U')
+*
+ INFO = 0
+ IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
+ INFO = 1
+ ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
+ INFO = 2
+ ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
+ + (.NOT.LSAME(TRANSA,'T')) .AND.
+ + (.NOT.LSAME(TRANSA,'C'))) THEN
+ INFO = 3
+ ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
+ INFO = 4
+ ELSE IF (M.LT.0) THEN
+ INFO = 5
+ ELSE IF (N.LT.0) THEN
+ INFO = 6
+ ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
+ INFO = 9
+ ELSE IF (LDB.LT.MAX(1,M)) THEN
+ INFO = 11
+ END IF
+ IF (INFO.NE.0) THEN
+ CALL XERBLA('DTRSM ',INFO)
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF (M.EQ.0 .OR. N.EQ.0) RETURN
+*
+* And when alpha.eq.zero.
+*
+ IF (ALPHA.EQ.ZERO) THEN
+ DO 20 J = 1,N
+ DO 10 I = 1,M
+ B(I,J) = ZERO
+ 10 CONTINUE
+ 20 CONTINUE
+ RETURN
+ END IF
+*
+* Start the operations.
+*
+ IF (LSIDE) THEN
+ IF (LSAME(TRANSA,'N')) THEN
+*
+* Form B := alpha*inv( A )*B.
+*
+ IF (UPPER) THEN
+ DO 60 J = 1,N
+ IF (ALPHA.NE.ONE) THEN
+ DO 30 I = 1,M
+ B(I,J) = ALPHA*B(I,J)
+ 30 CONTINUE
+ END IF
+ DO 50 K = M,1,-1
+ IF (B(K,J).NE.ZERO) THEN
+ IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
+ DO 40 I = 1,K - 1
+ B(I,J) = B(I,J) - B(K,J)*A(I,K)
+ 40 CONTINUE
+ END IF
+ 50 CONTINUE
+ 60 CONTINUE
+ ELSE
+ DO 100 J = 1,N
+ IF (ALPHA.NE.ONE) THEN
+ DO 70 I = 1,M
+ B(I,J) = ALPHA*B(I,J)
+ 70 CONTINUE
+ END IF
+ DO 90 K = 1,M
+ IF (B(K,J).NE.ZERO) THEN
+ IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
+ DO 80 I = K + 1,M
+ B(I,J) = B(I,J) - B(K,J)*A(I,K)
+ 80 CONTINUE
+ END IF
+ 90 CONTINUE
+ 100 CONTINUE
+ END IF
+ ELSE
+*
+* Form B := alpha*inv( A**T )*B.
+*
+ IF (UPPER) THEN
+ DO 130 J = 1,N
+ DO 120 I = 1,M
+ TEMP = ALPHA*B(I,J)
+ DO 110 K = 1,I - 1
+ TEMP = TEMP - A(K,I)*B(K,J)
+ 110 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/A(I,I)
+ B(I,J) = TEMP
+ 120 CONTINUE
+ 130 CONTINUE
+ ELSE
+ DO 160 J = 1,N
+ DO 150 I = M,1,-1
+ TEMP = ALPHA*B(I,J)
+ DO 140 K = I + 1,M
+ TEMP = TEMP - A(K,I)*B(K,J)
+ 140 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/A(I,I)
+ B(I,J) = TEMP
+ 150 CONTINUE
+ 160 CONTINUE
+ END IF
+ END IF
+ ELSE
+ IF (LSAME(TRANSA,'N')) THEN
+*
+* Form B := alpha*B*inv( A ).
+*
+ IF (UPPER) THEN
+ DO 210 J = 1,N
+ IF (ALPHA.NE.ONE) THEN
+ DO 170 I = 1,M
+ B(I,J) = ALPHA*B(I,J)
+ 170 CONTINUE
+ END IF
+ DO 190 K = 1,J - 1
+ IF (A(K,J).NE.ZERO) THEN
+ DO 180 I = 1,M
+ B(I,J) = B(I,J) - A(K,J)*B(I,K)
+ 180 CONTINUE
+ END IF
+ 190 CONTINUE
+ IF (NOUNIT) THEN
+ TEMP = ONE/A(J,J)
+ DO 200 I = 1,M
+ B(I,J) = TEMP*B(I,J)
+ 200 CONTINUE
+ END IF
+ 210 CONTINUE
+ ELSE
+ DO 260 J = N,1,-1
+ IF (ALPHA.NE.ONE) THEN
+ DO 220 I = 1,M
+ B(I,J) = ALPHA*B(I,J)
+ 220 CONTINUE
+ END IF
+ DO 240 K = J + 1,N
+ IF (A(K,J).NE.ZERO) THEN
+ DO 230 I = 1,M
+ B(I,J) = B(I,J) - A(K,J)*B(I,K)
+ 230 CONTINUE
+ END IF
+ 240 CONTINUE
+ IF (NOUNIT) THEN
+ TEMP = ONE/A(J,J)
+ DO 250 I = 1,M
+ B(I,J) = TEMP*B(I,J)
+ 250 CONTINUE
+ END IF
+ 260 CONTINUE
+ END IF
+ ELSE
+*
+* Form B := alpha*B*inv( A**T ).
+*
+ IF (UPPER) THEN
+ DO 310 K = N,1,-1
+ IF (NOUNIT) THEN
+ TEMP = ONE/A(K,K)
+ DO 270 I = 1,M
+ B(I,K) = TEMP*B(I,K)
+ 270 CONTINUE
+ END IF
+ DO 290 J = 1,K - 1
+ IF (A(J,K).NE.ZERO) THEN
+ TEMP = A(J,K)
+ DO 280 I = 1,M
+ B(I,J) = B(I,J) - TEMP*B(I,K)
+ 280 CONTINUE
+ END IF
+ 290 CONTINUE
+ IF (ALPHA.NE.ONE) THEN
+ DO 300 I = 1,M
+ B(I,K) = ALPHA*B(I,K)
+ 300 CONTINUE
+ END IF
+ 310 CONTINUE
+ ELSE
+ DO 360 K = 1,N
+ IF (NOUNIT) THEN
+ TEMP = ONE/A(K,K)
+ DO 320 I = 1,M
+ B(I,K) = TEMP*B(I,K)
+ 320 CONTINUE
+ END IF
+ DO 340 J = K + 1,N
+ IF (A(J,K).NE.ZERO) THEN
+ TEMP = A(J,K)
+ DO 330 I = 1,M
+ B(I,J) = B(I,J) - TEMP*B(I,K)
+ 330 CONTINUE
+ END IF
+ 340 CONTINUE
+ IF (ALPHA.NE.ONE) THEN
+ DO 350 I = 1,M
+ B(I,K) = ALPHA*B(I,K)
+ 350 CONTINUE
+ END IF
+ 360 CONTINUE
+ END IF
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of DTRSM
+*
+ END