From 4c5ae335d43d7af0738233de12a9ce0718b6a7db Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 31 May 2013 15:36:13 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9990
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
lib/linalg/disnan.f | 80 +++++++++
lib/linalg/dlaisnan.f | 91 ++++++++++
lib/linalg/zdotc.f | 101 +++++++++++
lib/linalg/zdscal.f | 94 ++++++++++
lib/linalg/zhpr.f | 279 ++++++++++++++++++++++++++++++
lib/linalg/zpptrf.f | 241 ++++++++++++++++++++++++++
lib/linalg/zpptri.f | 190 ++++++++++++++++++++
lib/linalg/zscal.f | 91 ++++++++++
lib/linalg/ztpmv.f | 388 +++++++++++++++++++++++++++++++++++++++++
lib/linalg/ztpsv.f | 390 ++++++++++++++++++++++++++++++++++++++++++
lib/linalg/ztptri.f | 242 ++++++++++++++++++++++++++
11 files changed, 2187 insertions(+)
create mode 100644 lib/linalg/disnan.f
create mode 100644 lib/linalg/dlaisnan.f
create mode 100644 lib/linalg/zdotc.f
create mode 100644 lib/linalg/zdscal.f
create mode 100644 lib/linalg/zhpr.f
create mode 100644 lib/linalg/zpptrf.f
create mode 100644 lib/linalg/zpptri.f
create mode 100644 lib/linalg/zscal.f
create mode 100644 lib/linalg/ztpmv.f
create mode 100644 lib/linalg/ztpsv.f
create mode 100644 lib/linalg/ztptri.f
diff --git a/lib/linalg/disnan.f b/lib/linalg/disnan.f
new file mode 100644
index 0000000000..355b827955
--- /dev/null
+++ b/lib/linalg/disnan.f
@@ -0,0 +1,80 @@
+*> \brief \b DISNAN tests input for NaN.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DISNAN + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/disnan.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/disnan.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/disnan.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* LOGICAL FUNCTION DISNAN( DIN )
+*
+* .. Scalar Arguments ..
+* DOUBLE PRECISION DIN
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> DISNAN returns .TRUE. if its argument is NaN, and .FALSE.
+*> otherwise. To be replaced by the Fortran 2003 intrinsic in the
+*> future.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] DIN
+*> \verbatim
+*> DIN is DOUBLE PRECISION
+*> Input to test for NaN.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date September 2012
+*
+*> \ingroup auxOTHERauxiliary
+*
+* =====================================================================
+ LOGICAL FUNCTION DISNAN( DIN )
+*
+* -- LAPACK auxiliary routine (version 3.4.2) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* September 2012
+*
+* .. Scalar Arguments ..
+ DOUBLE PRECISION DIN
+* ..
+*
+* =====================================================================
+*
+* .. External Functions ..
+ LOGICAL DLAISNAN
+ EXTERNAL DLAISNAN
+* ..
+* .. Executable Statements ..
+ DISNAN = DLAISNAN(DIN,DIN)
+ RETURN
+ END
diff --git a/lib/linalg/dlaisnan.f b/lib/linalg/dlaisnan.f
new file mode 100644
index 0000000000..58595c5c33
--- /dev/null
+++ b/lib/linalg/dlaisnan.f
@@ -0,0 +1,91 @@
+*> \brief \b DLAISNAN tests input for NaN by comparing two arguments for inequality.
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download DLAISNAN + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaisnan.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaisnan.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaisnan.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* LOGICAL FUNCTION DLAISNAN( DIN1, DIN2 )
+*
+* .. Scalar Arguments ..
+* DOUBLE PRECISION DIN1, DIN2
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> This routine is not for general use. It exists solely to avoid
+*> over-optimization in DISNAN.
+*>
+*> DLAISNAN checks for NaNs by comparing its two arguments for
+*> inequality. NaN is the only floating-point value where NaN != NaN
+*> returns .TRUE. To check for NaNs, pass the same variable as both
+*> arguments.
+*>
+*> A compiler must assume that the two arguments are
+*> not the same variable, and the test will not be optimized away.
+*> Interprocedural or whole-program optimization may delete this
+*> test. The ISNAN functions will be replaced by the correct
+*> Fortran 03 intrinsic once the intrinsic is widely available.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] DIN1
+*> \verbatim
+*> DIN1 is DOUBLE PRECISION
+*> \endverbatim
+*>
+*> \param[in] DIN2
+*> \verbatim
+*> DIN2 is DOUBLE PRECISION
+*> Two numbers to compare for inequality.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date September 2012
+*
+*> \ingroup auxOTHERauxiliary
+*
+* =====================================================================
+ LOGICAL FUNCTION DLAISNAN( DIN1, DIN2 )
+*
+* -- LAPACK auxiliary routine (version 3.4.2) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* September 2012
+*
+* .. Scalar Arguments ..
+ DOUBLE PRECISION DIN1, DIN2
+* ..
+*
+* =====================================================================
+*
+* .. Executable Statements ..
+ DLAISNAN = (DIN1.NE.DIN2)
+ RETURN
+ END
diff --git a/lib/linalg/zdotc.f b/lib/linalg/zdotc.f
new file mode 100644
index 0000000000..660648bbe1
--- /dev/null
+++ b/lib/linalg/zdotc.f
@@ -0,0 +1,101 @@
+*> \brief \b ZDOTC
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* COMPLEX*16 FUNCTION ZDOTC(N,ZX,INCX,ZY,INCY)
+*
+* .. Scalar Arguments ..
+* INTEGER INCX,INCY,N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 ZX(*),ZY(*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZDOTC forms the dot product of a vector.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complex16_blas_level1
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> jack dongarra, 3/11/78.
+*> modified 12/3/93, array(1) declarations changed to array(*)
+*> \endverbatim
+*>
+* =====================================================================
+ COMPLEX*16 FUNCTION ZDOTC(N,ZX,INCX,ZY,INCY)
+*
+* -- Reference BLAS level1 routine (version 3.4.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ INTEGER INCX,INCY,N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 ZX(*),ZY(*)
+* ..
+*
+* =====================================================================
+*
+* .. Local Scalars ..
+ COMPLEX*16 ZTEMP
+ INTEGER I,IX,IY
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DCONJG
+* ..
+ ZTEMP = (0.0d0,0.0d0)
+ ZDOTC = (0.0d0,0.0d0)
+ 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
+ ZTEMP = ZTEMP + DCONJG(ZX(I))*ZY(I)
+ 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
+ ZTEMP = ZTEMP + DCONJG(ZX(IX))*ZY(IY)
+ IX = IX + INCX
+ IY = IY + INCY
+ END DO
+ END IF
+ ZDOTC = ZTEMP
+ RETURN
+ END
diff --git a/lib/linalg/zdscal.f b/lib/linalg/zdscal.f
new file mode 100644
index 0000000000..57a9490237
--- /dev/null
+++ b/lib/linalg/zdscal.f
@@ -0,0 +1,94 @@
+*> \brief \b ZDSCAL
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZDSCAL(N,DA,ZX,INCX)
+*
+* .. Scalar Arguments ..
+* DOUBLE PRECISION DA
+* INTEGER INCX,N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 ZX(*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZDSCAL scales a vector by a constant.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complex16_blas_level1
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> jack dongarra, 3/11/78.
+*> modified 3/93 to return if incx .le. 0.
+*> modified 12/3/93, array(1) declarations changed to array(*)
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE ZDSCAL(N,DA,ZX,INCX)
+*
+* -- Reference BLAS level1 routine (version 3.4.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ DOUBLE PRECISION DA
+ INTEGER INCX,N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 ZX(*)
+* ..
+*
+* =====================================================================
+*
+* .. Local Scalars ..
+ INTEGER I,NINCX
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DCMPLX
+* ..
+ IF (N.LE.0 .OR. INCX.LE.0) RETURN
+ IF (INCX.EQ.1) THEN
+*
+* code for increment equal to 1
+*
+ DO I = 1,N
+ ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
+ END DO
+ ELSE
+*
+* code for increment not equal to 1
+*
+ NINCX = N*INCX
+ DO I = 1,NINCX,INCX
+ ZX(I) = DCMPLX(DA,0.0d0)*ZX(I)
+ END DO
+ END IF
+ RETURN
+ END
diff --git a/lib/linalg/zhpr.f b/lib/linalg/zhpr.f
new file mode 100644
index 0000000000..42e61196ba
--- /dev/null
+++ b/lib/linalg/zhpr.f
@@ -0,0 +1,279 @@
+*> \brief \b ZHPR
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZHPR(UPLO,N,ALPHA,X,INCX,AP)
+*
+* .. Scalar Arguments ..
+* DOUBLE PRECISION ALPHA
+* INTEGER INCX,N
+* CHARACTER UPLO
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 AP(*),X(*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZHPR performs the hermitian rank 1 operation
+*>
+*> A := alpha*x*x**H + A,
+*>
+*> where alpha is a real scalar, x is an n element vector and A is an
+*> n by n hermitian 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 COMPLEX*16 array of dimension at least
+*> ( 1 + ( n - 1 )*abs( INCX ) ).
+*> Before entry, the incremented array X must contain the n
+*> element vector x.
+*> \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 COMPLEX*16 array of DIMENSION at least
+*> ( ( n*( n + 1 ) )/2 ).
+*> Before entry with UPLO = 'U' or 'u', the array AP must
+*> contain the upper triangular part of the hermitian matrix
+*> packed sequentially, column by column, so that AP( 1 )
+*> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
+*> and a( 2, 2 ) respectively, and so on. 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 hermitian matrix
+*> packed sequentially, column by column, so that AP( 1 )
+*> contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
+*> and a( 3, 1 ) respectively, and so on. On exit, the array
+*> AP is overwritten by the lower triangular part of the
+*> updated matrix.
+*> Note that the imaginary parts of the diagonal elements need
+*> not be set, they are assumed to be zero, and on exit they
+*> are set to zero.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complex16_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 ZHPR(UPLO,N,ALPHA,X,INCX,AP)
+*
+* -- Reference BLAS level2 routine (version 3.4.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ DOUBLE PRECISION ALPHA
+ INTEGER INCX,N
+ CHARACTER UPLO
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 AP(*),X(*)
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 ZERO
+ PARAMETER (ZERO= (0.0D+0,0.0D+0))
+* ..
+* .. Local Scalars ..
+ COMPLEX*16 TEMP
+ INTEGER I,INFO,IX,J,JX,K,KK,KX
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DBLE,DCONJG
+* ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
+ INFO = 1
+ ELSE IF (N.LT.0) THEN
+ INFO = 2
+ ELSE IF (INCX.EQ.0) THEN
+ INFO = 5
+ END IF
+ IF (INFO.NE.0) THEN
+ CALL XERBLA('ZHPR ',INFO)
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF ((N.EQ.0) .OR. (ALPHA.EQ.DBLE(ZERO))) RETURN
+*
+* Set the start point in X if the increment is not unity.
+*
+ IF (INCX.LE.0) THEN
+ KX = 1 - (N-1)*INCX
+ ELSE IF (INCX.NE.1) THEN
+ KX = 1
+ END IF
+*
+* Start the operations. In this version the elements of the array AP
+* are accessed sequentially with one pass through AP.
+*
+ KK = 1
+ IF (LSAME(UPLO,'U')) THEN
+*
+* Form A when upper triangle is stored in AP.
+*
+ IF (INCX.EQ.1) THEN
+ DO 20 J = 1,N
+ IF (X(J).NE.ZERO) THEN
+ TEMP = ALPHA*DCONJG(X(J))
+ K = KK
+ DO 10 I = 1,J - 1
+ AP(K) = AP(K) + X(I)*TEMP
+ K = K + 1
+ 10 CONTINUE
+ AP(KK+J-1) = DBLE(AP(KK+J-1)) + DBLE(X(J)*TEMP)
+ ELSE
+ AP(KK+J-1) = DBLE(AP(KK+J-1))
+ END IF
+ KK = KK + J
+ 20 CONTINUE
+ ELSE
+ JX = KX
+ DO 40 J = 1,N
+ IF (X(JX).NE.ZERO) THEN
+ TEMP = ALPHA*DCONJG(X(JX))
+ IX = KX
+ DO 30 K = KK,KK + J - 2
+ AP(K) = AP(K) + X(IX)*TEMP
+ IX = IX + INCX
+ 30 CONTINUE
+ AP(KK+J-1) = DBLE(AP(KK+J-1)) + DBLE(X(JX)*TEMP)
+ ELSE
+ AP(KK+J-1) = DBLE(AP(KK+J-1))
+ END IF
+ JX = JX + INCX
+ KK = KK + J
+ 40 CONTINUE
+ END IF
+ ELSE
+*
+* Form A when lower triangle is stored in AP.
+*
+ IF (INCX.EQ.1) THEN
+ DO 60 J = 1,N
+ IF (X(J).NE.ZERO) THEN
+ TEMP = ALPHA*DCONJG(X(J))
+ AP(KK) = DBLE(AP(KK)) + DBLE(TEMP*X(J))
+ K = KK + 1
+ DO 50 I = J + 1,N
+ AP(K) = AP(K) + X(I)*TEMP
+ K = K + 1
+ 50 CONTINUE
+ ELSE
+ AP(KK) = DBLE(AP(KK))
+ END IF
+ KK = KK + N - J + 1
+ 60 CONTINUE
+ ELSE
+ JX = KX
+ DO 80 J = 1,N
+ IF (X(JX).NE.ZERO) THEN
+ TEMP = ALPHA*DCONJG(X(JX))
+ AP(KK) = DBLE(AP(KK)) + DBLE(TEMP*X(JX))
+ IX = JX
+ DO 70 K = KK + 1,KK + N - J
+ IX = IX + INCX
+ AP(K) = AP(K) + X(IX)*TEMP
+ 70 CONTINUE
+ ELSE
+ AP(KK) = DBLE(AP(KK))
+ END IF
+ JX = JX + INCX
+ KK = KK + N - J + 1
+ 80 CONTINUE
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of ZHPR .
+*
+ END
diff --git a/lib/linalg/zpptrf.f b/lib/linalg/zpptrf.f
new file mode 100644
index 0000000000..c34aff332a
--- /dev/null
+++ b/lib/linalg/zpptrf.f
@@ -0,0 +1,241 @@
+*> \brief \b ZPPTRF
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZPPTRF + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpptrf.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpptrf.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpptrf.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZPPTRF( UPLO, N, AP, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER INFO, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 AP( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZPPTRF computes the Cholesky factorization of a complex Hermitian
+*> positive definite matrix A stored in packed format.
+*>
+*> The factorization has the form
+*> A = U**H * U, if UPLO = 'U', or
+*> A = L * L**H, if UPLO = 'L',
+*> where U is an upper triangular matrix and L is lower triangular.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> = 'U': Upper triangle of A is stored;
+*> = 'L': Lower triangle of A is stored.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] AP
+*> \verbatim
+*> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
+*> On entry, the upper or lower triangle of the Hermitian matrix
+*> A, packed columnwise in a linear array. The j-th column of A
+*> is stored in the array AP as follows:
+*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
+*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
+*> See below for further details.
+*>
+*> On exit, if INFO = 0, the triangular factor U or L from the
+*> Cholesky factorization A = U**H*U or A = L*L**H, in the same
+*> storage format as A.
+*> \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, the leading minor of order i is not
+*> positive definite, and the factorization could not be
+*> completed.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complex16OTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> The packed storage scheme is illustrated by the following example
+*> when N = 4, UPLO = 'U':
+*>
+*> Two-dimensional storage of the Hermitian matrix A:
+*>
+*> a11 a12 a13 a14
+*> a22 a23 a24
+*> a33 a34 (aij = conjg(aji))
+*> a44
+*>
+*> Packed storage of the upper triangle of A:
+*>
+*> AP = [ a11, a12, a22, a13, a23, a33, a14, a24, a34, a44 ]
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE ZPPTRF( UPLO, N, AP, INFO )
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 AP( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER J, JC, JJ
+ DOUBLE PRECISION AJJ
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ COMPLEX*16 ZDOTC
+ EXTERNAL LSAME, ZDOTC
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPSV
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DBLE, SQRT
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZPPTRF', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+ IF( UPPER ) THEN
+*
+* Compute the Cholesky factorization A = U**H * U.
+*
+ JJ = 0
+ DO 10 J = 1, N
+ JC = JJ + 1
+ JJ = JJ + J
+*
+* Compute elements 1:J-1 of column J.
+*
+ IF( J.GT.1 )
+ $ CALL ZTPSV( 'Upper', 'Conjugate transpose', 'Non-unit',
+ $ J-1, AP, AP( JC ), 1 )
+*
+* Compute U(J,J) and test for non-positive-definiteness.
+*
+ AJJ = DBLE( AP( JJ ) ) - ZDOTC( J-1, AP( JC ), 1, AP( JC ),
+ $ 1 )
+ IF( AJJ.LE.ZERO ) THEN
+ AP( JJ ) = AJJ
+ GO TO 30
+ END IF
+ AP( JJ ) = SQRT( AJJ )
+ 10 CONTINUE
+ ELSE
+*
+* Compute the Cholesky factorization A = L * L**H.
+*
+ JJ = 1
+ DO 20 J = 1, N
+*
+* Compute L(J,J) and test for non-positive-definiteness.
+*
+ AJJ = DBLE( AP( JJ ) )
+ IF( AJJ.LE.ZERO ) THEN
+ AP( JJ ) = AJJ
+ GO TO 30
+ END IF
+ AJJ = SQRT( AJJ )
+ AP( JJ ) = AJJ
+*
+* Compute elements J+1:N of column J and update the trailing
+* submatrix.
+*
+ IF( J.LT.N ) THEN
+ CALL ZDSCAL( N-J, ONE / AJJ, AP( JJ+1 ), 1 )
+ CALL ZHPR( 'Lower', N-J, -ONE, AP( JJ+1 ), 1,
+ $ AP( JJ+N-J+1 ) )
+ JJ = JJ + N - J + 1
+ END IF
+ 20 CONTINUE
+ END IF
+ GO TO 40
+*
+ 30 CONTINUE
+ INFO = J
+*
+ 40 CONTINUE
+ RETURN
+*
+* End of ZPPTRF
+*
+ END
diff --git a/lib/linalg/zpptri.f b/lib/linalg/zpptri.f
new file mode 100644
index 0000000000..0946797450
--- /dev/null
+++ b/lib/linalg/zpptri.f
@@ -0,0 +1,190 @@
+*> \brief \b ZPPTRI
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZPPTRI + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpptri.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpptri.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpptri.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZPPTRI( UPLO, N, AP, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER INFO, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 AP( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZPPTRI computes the inverse of a complex Hermitian positive definite
+*> matrix A using the Cholesky factorization A = U**H*U or A = L*L**H
+*> computed by ZPPTRF.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> = 'U': Upper triangular factor is stored in AP;
+*> = 'L': Lower triangular factor is stored in AP.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] AP
+*> \verbatim
+*> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
+*> On entry, the triangular factor U or L from the Cholesky
+*> factorization A = U**H*U or A = L*L**H, packed columnwise as
+*> a linear array. The j-th column of U or L is stored in the
+*> array AP as follows:
+*> if UPLO = 'U', AP(i + (j-1)*j/2) = U(i,j) for 1<=i<=j;
+*> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = L(i,j) for j<=i<=n.
+*>
+*> On exit, the upper or lower triangle of the (Hermitian)
+*> inverse of A, overwriting the input factor U or L.
+*> \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, the (i,i) element of the factor U or L is
+*> zero, and the inverse could not be computed.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complex16OTHERcomputational
+*
+* =====================================================================
+ SUBROUTINE ZPPTRI( UPLO, N, AP, INFO )
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 AP( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ONE
+ PARAMETER ( ONE = 1.0D+0 )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER J, JC, JJ, JJN
+ DOUBLE PRECISION AJJ
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ COMPLEX*16 ZDOTC
+ EXTERNAL LSAME, ZDOTC
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZDSCAL, ZHPR, ZTPMV, ZTPTRI
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DBLE
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZPPTRI', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Invert the triangular Cholesky factor U or L.
+*
+ CALL ZTPTRI( UPLO, 'Non-unit', N, AP, INFO )
+ IF( INFO.GT.0 )
+ $ RETURN
+ IF( UPPER ) THEN
+*
+* Compute the product inv(U) * inv(U)**H.
+*
+ JJ = 0
+ DO 10 J = 1, N
+ JC = JJ + 1
+ JJ = JJ + J
+ IF( J.GT.1 )
+ $ CALL ZHPR( 'Upper', J-1, ONE, AP( JC ), 1, AP )
+ AJJ = AP( JJ )
+ CALL ZDSCAL( J, AJJ, AP( JC ), 1 )
+ 10 CONTINUE
+*
+ ELSE
+*
+* Compute the product inv(L)**H * inv(L).
+*
+ JJ = 1
+ DO 20 J = 1, N
+ JJN = JJ + N - J + 1
+ AP( JJ ) = DBLE( ZDOTC( N-J+1, AP( JJ ), 1, AP( JJ ), 1 ) )
+ IF( J.LT.N )
+ $ CALL ZTPMV( 'Lower', 'Conjugate transpose', 'Non-unit',
+ $ N-J, AP( JJN ), AP( JJ+1 ), 1 )
+ JJ = JJN
+ 20 CONTINUE
+ END IF
+*
+ RETURN
+*
+* End of ZPPTRI
+*
+ END
diff --git a/lib/linalg/zscal.f b/lib/linalg/zscal.f
new file mode 100644
index 0000000000..ad28a10a9b
--- /dev/null
+++ b/lib/linalg/zscal.f
@@ -0,0 +1,91 @@
+*> \brief \b ZSCAL
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZSCAL(N,ZA,ZX,INCX)
+*
+* .. Scalar Arguments ..
+* COMPLEX*16 ZA
+* INTEGER INCX,N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 ZX(*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZSCAL scales a vector by a constant.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complex16_blas_level1
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> jack dongarra, 3/11/78.
+*> modified 3/93 to return if incx .le. 0.
+*> modified 12/3/93, array(1) declarations changed to array(*)
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE ZSCAL(N,ZA,ZX,INCX)
+*
+* -- Reference BLAS level1 routine (version 3.4.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ COMPLEX*16 ZA
+ INTEGER INCX,N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 ZX(*)
+* ..
+*
+* =====================================================================
+*
+* .. Local Scalars ..
+ INTEGER I,NINCX
+* ..
+ IF (N.LE.0 .OR. INCX.LE.0) RETURN
+ IF (INCX.EQ.1) THEN
+*
+* code for increment equal to 1
+*
+ DO I = 1,N
+ ZX(I) = ZA*ZX(I)
+ END DO
+ ELSE
+*
+* code for increment not equal to 1
+*
+ NINCX = N*INCX
+ DO I = 1,NINCX,INCX
+ ZX(I) = ZA*ZX(I)
+ END DO
+ END IF
+ RETURN
+ END
diff --git a/lib/linalg/ztpmv.f b/lib/linalg/ztpmv.f
new file mode 100644
index 0000000000..e277ec1a6e
--- /dev/null
+++ b/lib/linalg/ztpmv.f
@@ -0,0 +1,388 @@
+*> \brief \b ZTPMV
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+*
+* .. Scalar Arguments ..
+* INTEGER INCX,N
+* CHARACTER DIAG,TRANS,UPLO
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 AP(*),X(*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZTPMV performs one of the matrix-vector operations
+*>
+*> x := A*x, or x := A**T*x, or x := A**H*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**H*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 COMPLEX*16 array of DIMENSION at least
+*> ( ( n*( n + 1 ) )/2 ).
+*> Before entry with UPLO = 'U' or 'u', the array AP must
+*> contain the upper triangular 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] X
+*> \verbatim
+*> X is (input/output) COMPLEX*16 array of dimension at least
+*> ( 1 + ( n - 1 )*abs( INCX ) ).
+*> Before entry, the incremented array X must contain the n
+*> element vector x. On exit, X is overwritten with the
+*> tranformed vector x.
+*> \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.
+*
+*> \date November 2011
+*
+*> \ingroup complex16_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 ZTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+*
+* -- Reference BLAS level2 routine (version 3.4.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ INTEGER INCX,N
+ CHARACTER DIAG,TRANS,UPLO
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 AP(*),X(*)
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 ZERO
+ PARAMETER (ZERO= (0.0D+0,0.0D+0))
+* ..
+* .. Local Scalars ..
+ COMPLEX*16 TEMP
+ INTEGER I,INFO,IX,J,JX,K,KK,KX
+ LOGICAL NOCONJ,NOUNIT
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DCONJG
+* ..
+*
+* 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('ZTPMV ',INFO)
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF (N.EQ.0) RETURN
+*
+ NOCONJ = LSAME(TRANS,'T')
+ NOUNIT = LSAME(DIAG,'N')
+*
+* Set up the start point in X if the increment is not unity. This
+* will be ( N - 1 )*INCX too small for descending loops.
+*
+ IF (INCX.LE.0) THEN
+ KX = 1 - (N-1)*INCX
+ ELSE IF (INCX.NE.1) THEN
+ KX = 1
+ END IF
+*
+* Start the operations. In this version the elements of AP are
+* accessed sequentially with one pass through AP.
+*
+ IF (LSAME(TRANS,'N')) THEN
+*
+* Form x:= A*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ KK = 1
+ IF (INCX.EQ.1) THEN
+ DO 20 J = 1,N
+ IF (X(J).NE.ZERO) THEN
+ TEMP = X(J)
+ K = KK
+ DO 10 I = 1,J - 1
+ X(I) = X(I) + TEMP*AP(K)
+ K = K + 1
+ 10 CONTINUE
+ IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
+ END IF
+ KK = KK + J
+ 20 CONTINUE
+ ELSE
+ JX = KX
+ DO 40 J = 1,N
+ IF (X(JX).NE.ZERO) THEN
+ TEMP = X(JX)
+ IX = KX
+ DO 30 K = KK,KK + J - 2
+ X(IX) = X(IX) + TEMP*AP(K)
+ IX = IX + INCX
+ 30 CONTINUE
+ IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
+ END IF
+ JX = JX + INCX
+ KK = KK + J
+ 40 CONTINUE
+ END IF
+ ELSE
+ KK = (N* (N+1))/2
+ IF (INCX.EQ.1) THEN
+ DO 60 J = N,1,-1
+ IF (X(J).NE.ZERO) THEN
+ TEMP = X(J)
+ K = KK
+ DO 50 I = N,J + 1,-1
+ X(I) = X(I) + TEMP*AP(K)
+ K = K - 1
+ 50 CONTINUE
+ IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
+ END IF
+ KK = KK - (N-J+1)
+ 60 CONTINUE
+ ELSE
+ KX = KX + (N-1)*INCX
+ JX = KX
+ DO 80 J = N,1,-1
+ IF (X(JX).NE.ZERO) THEN
+ TEMP = X(JX)
+ IX = KX
+ DO 70 K = KK,KK - (N- (J+1)),-1
+ X(IX) = X(IX) + TEMP*AP(K)
+ IX = IX - INCX
+ 70 CONTINUE
+ IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
+ END IF
+ JX = JX - INCX
+ KK = KK - (N-J+1)
+ 80 CONTINUE
+ END IF
+ END IF
+ ELSE
+*
+* Form x := A**T*x or x := A**H*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ KK = (N* (N+1))/2
+ IF (INCX.EQ.1) THEN
+ DO 110 J = N,1,-1
+ TEMP = X(J)
+ K = KK - 1
+ IF (NOCONJ) THEN
+ IF (NOUNIT) TEMP = TEMP*AP(KK)
+ DO 90 I = J - 1,1,-1
+ TEMP = TEMP + AP(K)*X(I)
+ K = K - 1
+ 90 CONTINUE
+ ELSE
+ IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
+ DO 100 I = J - 1,1,-1
+ TEMP = TEMP + DCONJG(AP(K))*X(I)
+ K = K - 1
+ 100 CONTINUE
+ END IF
+ X(J) = TEMP
+ KK = KK - J
+ 110 CONTINUE
+ ELSE
+ JX = KX + (N-1)*INCX
+ DO 140 J = N,1,-1
+ TEMP = X(JX)
+ IX = JX
+ IF (NOCONJ) THEN
+ IF (NOUNIT) TEMP = TEMP*AP(KK)
+ DO 120 K = KK - 1,KK - J + 1,-1
+ IX = IX - INCX
+ TEMP = TEMP + AP(K)*X(IX)
+ 120 CONTINUE
+ ELSE
+ IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
+ DO 130 K = KK - 1,KK - J + 1,-1
+ IX = IX - INCX
+ TEMP = TEMP + DCONJG(AP(K))*X(IX)
+ 130 CONTINUE
+ END IF
+ X(JX) = TEMP
+ JX = JX - INCX
+ KK = KK - J
+ 140 CONTINUE
+ END IF
+ ELSE
+ KK = 1
+ IF (INCX.EQ.1) THEN
+ DO 170 J = 1,N
+ TEMP = X(J)
+ K = KK + 1
+ IF (NOCONJ) THEN
+ IF (NOUNIT) TEMP = TEMP*AP(KK)
+ DO 150 I = J + 1,N
+ TEMP = TEMP + AP(K)*X(I)
+ K = K + 1
+ 150 CONTINUE
+ ELSE
+ IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
+ DO 160 I = J + 1,N
+ TEMP = TEMP + DCONJG(AP(K))*X(I)
+ K = K + 1
+ 160 CONTINUE
+ END IF
+ X(J) = TEMP
+ KK = KK + (N-J+1)
+ 170 CONTINUE
+ ELSE
+ JX = KX
+ DO 200 J = 1,N
+ TEMP = X(JX)
+ IX = JX
+ IF (NOCONJ) THEN
+ IF (NOUNIT) TEMP = TEMP*AP(KK)
+ DO 180 K = KK + 1,KK + N - J
+ IX = IX + INCX
+ TEMP = TEMP + AP(K)*X(IX)
+ 180 CONTINUE
+ ELSE
+ IF (NOUNIT) TEMP = TEMP*DCONJG(AP(KK))
+ DO 190 K = KK + 1,KK + N - J
+ IX = IX + INCX
+ TEMP = TEMP + DCONJG(AP(K))*X(IX)
+ 190 CONTINUE
+ END IF
+ X(JX) = TEMP
+ JX = JX + INCX
+ KK = KK + (N-J+1)
+ 200 CONTINUE
+ END IF
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of ZTPMV .
+*
+ END
diff --git a/lib/linalg/ztpsv.f b/lib/linalg/ztpsv.f
new file mode 100644
index 0000000000..0e75f9facf
--- /dev/null
+++ b/lib/linalg/ztpsv.f
@@ -0,0 +1,390 @@
+*> \brief \b ZTPSV
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+*
+* .. Scalar Arguments ..
+* INTEGER INCX,N
+* CHARACTER DIAG,TRANS,UPLO
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 AP(*),X(*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZTPSV solves one of the systems of equations
+*>
+*> A*x = b, or A**T*x = b, or A**H*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**H*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 COMPLEX*16 array of DIMENSION at least
+*> ( ( n*( n + 1 ) )/2 ).
+*> Before entry with UPLO = 'U' or 'u', the array AP must
+*> contain the upper triangular 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 COMPLEX*16 array of dimension at least
+*> ( 1 + ( n - 1 )*abs( INCX ) ).
+*> Before entry, the incremented array X must contain the n
+*> element right-hand side vector b. On exit, X is overwritten
+*> with the solution vector x.
+*> \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.
+*
+*> \date November 2011
+*
+*> \ingroup complex16_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 ZTPSV(UPLO,TRANS,DIAG,N,AP,X,INCX)
+*
+* -- Reference BLAS level2 routine (version 3.4.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ INTEGER INCX,N
+ CHARACTER DIAG,TRANS,UPLO
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 AP(*),X(*)
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 ZERO
+ PARAMETER (ZERO= (0.0D+0,0.0D+0))
+* ..
+* .. Local Scalars ..
+ COMPLEX*16 TEMP
+ INTEGER I,INFO,IX,J,JX,K,KK,KX
+ LOGICAL NOCONJ,NOUNIT
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DCONJG
+* ..
+*
+* 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('ZTPSV ',INFO)
+ RETURN
+ END IF
+*
+* Quick return if possible.
+*
+ IF (N.EQ.0) RETURN
+*
+ NOCONJ = LSAME(TRANS,'T')
+ NOUNIT = LSAME(DIAG,'N')
+*
+* Set up the start point in X if the increment is not unity. This
+* will be ( N - 1 )*INCX too small for descending loops.
+*
+ IF (INCX.LE.0) THEN
+ KX = 1 - (N-1)*INCX
+ ELSE IF (INCX.NE.1) THEN
+ KX = 1
+ END IF
+*
+* Start the operations. In this version the elements of AP are
+* accessed sequentially with one pass through AP.
+*
+ IF (LSAME(TRANS,'N')) THEN
+*
+* Form x := inv( A )*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ KK = (N* (N+1))/2
+ IF (INCX.EQ.1) THEN
+ DO 20 J = N,1,-1
+ IF (X(J).NE.ZERO) THEN
+ IF (NOUNIT) X(J) = X(J)/AP(KK)
+ TEMP = X(J)
+ K = KK - 1
+ DO 10 I = J - 1,1,-1
+ X(I) = X(I) - TEMP*AP(K)
+ K = K - 1
+ 10 CONTINUE
+ END IF
+ KK = KK - J
+ 20 CONTINUE
+ ELSE
+ JX = KX + (N-1)*INCX
+ DO 40 J = N,1,-1
+ IF (X(JX).NE.ZERO) THEN
+ IF (NOUNIT) X(JX) = X(JX)/AP(KK)
+ TEMP = X(JX)
+ IX = JX
+ DO 30 K = KK - 1,KK - J + 1,-1
+ IX = IX - INCX
+ X(IX) = X(IX) - TEMP*AP(K)
+ 30 CONTINUE
+ END IF
+ JX = JX - INCX
+ KK = KK - J
+ 40 CONTINUE
+ END IF
+ ELSE
+ KK = 1
+ IF (INCX.EQ.1) THEN
+ DO 60 J = 1,N
+ IF (X(J).NE.ZERO) THEN
+ IF (NOUNIT) X(J) = X(J)/AP(KK)
+ TEMP = X(J)
+ K = KK + 1
+ DO 50 I = J + 1,N
+ X(I) = X(I) - TEMP*AP(K)
+ K = K + 1
+ 50 CONTINUE
+ END IF
+ KK = KK + (N-J+1)
+ 60 CONTINUE
+ ELSE
+ JX = KX
+ DO 80 J = 1,N
+ IF (X(JX).NE.ZERO) THEN
+ IF (NOUNIT) X(JX) = X(JX)/AP(KK)
+ TEMP = X(JX)
+ IX = JX
+ DO 70 K = KK + 1,KK + N - J
+ IX = IX + INCX
+ X(IX) = X(IX) - TEMP*AP(K)
+ 70 CONTINUE
+ END IF
+ JX = JX + INCX
+ KK = KK + (N-J+1)
+ 80 CONTINUE
+ END IF
+ END IF
+ ELSE
+*
+* Form x := inv( A**T )*x or x := inv( A**H )*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ KK = 1
+ IF (INCX.EQ.1) THEN
+ DO 110 J = 1,N
+ TEMP = X(J)
+ K = KK
+ IF (NOCONJ) THEN
+ 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)
+ ELSE
+ DO 100 I = 1,J - 1
+ TEMP = TEMP - DCONJG(AP(K))*X(I)
+ K = K + 1
+ 100 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK+J-1))
+ END IF
+ X(J) = TEMP
+ KK = KK + J
+ 110 CONTINUE
+ ELSE
+ JX = KX
+ DO 140 J = 1,N
+ TEMP = X(JX)
+ IX = KX
+ IF (NOCONJ) THEN
+ DO 120 K = KK,KK + J - 2
+ TEMP = TEMP - AP(K)*X(IX)
+ IX = IX + INCX
+ 120 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/AP(KK+J-1)
+ ELSE
+ DO 130 K = KK,KK + J - 2
+ TEMP = TEMP - DCONJG(AP(K))*X(IX)
+ IX = IX + INCX
+ 130 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK+J-1))
+ END IF
+ X(JX) = TEMP
+ JX = JX + INCX
+ KK = KK + J
+ 140 CONTINUE
+ END IF
+ ELSE
+ KK = (N* (N+1))/2
+ IF (INCX.EQ.1) THEN
+ DO 170 J = N,1,-1
+ TEMP = X(J)
+ K = KK
+ IF (NOCONJ) THEN
+ DO 150 I = N,J + 1,-1
+ TEMP = TEMP - AP(K)*X(I)
+ K = K - 1
+ 150 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
+ ELSE
+ DO 160 I = N,J + 1,-1
+ TEMP = TEMP - DCONJG(AP(K))*X(I)
+ K = K - 1
+ 160 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK-N+J))
+ END IF
+ X(J) = TEMP
+ KK = KK - (N-J+1)
+ 170 CONTINUE
+ ELSE
+ KX = KX + (N-1)*INCX
+ JX = KX
+ DO 200 J = N,1,-1
+ TEMP = X(JX)
+ IX = KX
+ IF (NOCONJ) THEN
+ DO 180 K = KK,KK - (N- (J+1)),-1
+ TEMP = TEMP - AP(K)*X(IX)
+ IX = IX - INCX
+ 180 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/AP(KK-N+J)
+ ELSE
+ DO 190 K = KK,KK - (N- (J+1)),-1
+ TEMP = TEMP - DCONJG(AP(K))*X(IX)
+ IX = IX - INCX
+ 190 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(AP(KK-N+J))
+ END IF
+ X(JX) = TEMP
+ JX = JX - INCX
+ KK = KK - (N-J+1)
+ 200 CONTINUE
+ END IF
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of ZTPSV .
+*
+ END
diff --git a/lib/linalg/ztptri.f b/lib/linalg/ztptri.f
new file mode 100644
index 0000000000..187c9ccac1
--- /dev/null
+++ b/lib/linalg/ztptri.f
@@ -0,0 +1,242 @@
+*> \brief \b ZTPTRI
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZTPTRI + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztptri.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztptri.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztptri.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZTPTRI( UPLO, DIAG, N, AP, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER DIAG, UPLO
+* INTEGER INFO, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 AP( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZTPTRI computes the inverse of a complex upper or lower triangular
+*> matrix A stored in packed format.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> = 'U': A is upper triangular;
+*> = 'L': A is lower triangular.
+*> \endverbatim
+*>
+*> \param[in] DIAG
+*> \verbatim
+*> DIAG is CHARACTER*1
+*> = 'N': A is non-unit triangular;
+*> = 'U': A is unit triangular.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] AP
+*> \verbatim
+*> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
+*> On entry, the upper or lower triangular matrix A, stored
+*> columnwise in a linear array. The j-th column of A is stored
+*> in the array AP as follows:
+*> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
+*> if UPLO = 'L', AP(i + (j-1)*((2*n-j)/2) = A(i,j) for j<=i<=n.
+*> See below for further details.
+*> On exit, the (triangular) inverse of the original matrix, in
+*> the same packed storage format.
+*> \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, A(i,i) is exactly zero. The triangular
+*> matrix is singular and its inverse can not be computed.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2011
+*
+*> \ingroup complex16OTHERcomputational
+*
+*> \par Further Details:
+* =====================
+*>
+*> \verbatim
+*>
+*> A triangular matrix A can be transferred to packed storage using one
+*> of the following program segments:
+*>
+*> UPLO = 'U': UPLO = 'L':
+*>
+*> JC = 1 JC = 1
+*> DO 2 J = 1, N DO 2 J = 1, N
+*> DO 1 I = 1, J DO 1 I = J, N
+*> AP(JC+I-1) = A(I,J) AP(JC+I-J) = A(I,J)
+*> 1 CONTINUE 1 CONTINUE
+*> JC = JC + J JC = JC + N - J + 1
+*> 2 CONTINUE 2 CONTINUE
+*> \endverbatim
+*>
+* =====================================================================
+ SUBROUTINE ZTPTRI( UPLO, DIAG, N, AP, INFO )
+*
+* -- LAPACK computational routine (version 3.4.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2011
+*
+* .. Scalar Arguments ..
+ CHARACTER DIAG, UPLO
+ INTEGER INFO, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 AP( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 ONE, ZERO
+ PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ),
+ $ ZERO = ( 0.0D+0, 0.0D+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL NOUNIT, UPPER
+ INTEGER J, JC, JCLAST, JJ
+ COMPLEX*16 AJJ
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZSCAL, ZTPMV
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ NOUNIT = LSAME( DIAG, 'N' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -3
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZTPTRI', -INFO )
+ RETURN
+ END IF
+*
+* Check for singularity if non-unit.
+*
+ IF( NOUNIT ) THEN
+ IF( UPPER ) THEN
+ JJ = 0
+ DO 10 INFO = 1, N
+ JJ = JJ + INFO
+ IF( AP( JJ ).EQ.ZERO )
+ $ RETURN
+ 10 CONTINUE
+ ELSE
+ JJ = 1
+ DO 20 INFO = 1, N
+ IF( AP( JJ ).EQ.ZERO )
+ $ RETURN
+ JJ = JJ + N - INFO + 1
+ 20 CONTINUE
+ END IF
+ INFO = 0
+ END IF
+*
+ IF( UPPER ) THEN
+*
+* Compute inverse of upper triangular matrix.
+*
+ JC = 1
+ DO 30 J = 1, N
+ IF( NOUNIT ) THEN
+ AP( JC+J-1 ) = ONE / AP( JC+J-1 )
+ AJJ = -AP( JC+J-1 )
+ ELSE
+ AJJ = -ONE
+ END IF
+*
+* Compute elements 1:j-1 of j-th column.
+*
+ CALL ZTPMV( 'Upper', 'No transpose', DIAG, J-1, AP,
+ $ AP( JC ), 1 )
+ CALL ZSCAL( J-1, AJJ, AP( JC ), 1 )
+ JC = JC + J
+ 30 CONTINUE
+*
+ ELSE
+*
+* Compute inverse of lower triangular matrix.
+*
+ JC = N*( N+1 ) / 2
+ DO 40 J = N, 1, -1
+ IF( NOUNIT ) THEN
+ AP( JC ) = ONE / AP( JC )
+ AJJ = -AP( JC )
+ ELSE
+ AJJ = -ONE
+ END IF
+ IF( J.LT.N ) THEN
+*
+* Compute elements j+1:n of j-th column.
+*
+ CALL ZTPMV( 'Lower', 'No transpose', DIAG, N-J,
+ $ AP( JCLAST ), AP( JC+1 ), 1 )
+ CALL ZSCAL( N-J, AJJ, AP( JC+1 ), 1 )
+ END IF
+ JCLAST = JC
+ JC = JC - N + J - 2
+ 40 CONTINUE
+ END IF
+*
+ RETURN
+*
+* End of ZTPTRI
+*
+ END
--
GitLab