diff --git a/fortran/zhegs2.f b/fortran/zhegs2.f
new file mode 100644
index 0000000..aec5263
--- /dev/null
+++ b/fortran/zhegs2.f
@@ -0,0 +1,297 @@
+*> \brief \b ZHEGS2 reduces a Hermitian definite generalized eigenproblem to standard form, using the factorization results obtained from cpotrf (unblocked algorithm).
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZHEGS2 + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER INFO, ITYPE, LDA, LDB, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A( LDA, * ), B( LDB, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZHEGS2 reduces a complex Hermitian-definite generalized
+*> eigenproblem to standard form.
+*>
+*> If ITYPE = 1, the problem is A*x = lambda*B*x,
+*> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
+*>
+*> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
+*> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H *A*L.
+*>
+*> B must have been previously factorized as U**H *U or L*L**H by ZPOTRF.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] ITYPE
+*> \verbatim
+*> ITYPE is INTEGER
+*> = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
+*> = 2 or 3: compute U*A*U**H or L**H *A*L.
+*> \endverbatim
+*>
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> Specifies whether the upper or lower triangular part of the
+*> Hermitian matrix A is stored, and how B has been factorized.
+*> = 'U': Upper triangular
+*> = 'L': Lower triangular
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A and B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N)
+*> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
+*> n by n upper triangular part of A contains the upper
+*> triangular part of the matrix A, and the strictly lower
+*> triangular part of A is not referenced. If UPLO = 'L', the
+*> leading n by n lower triangular part of A contains the lower
+*> triangular part of the matrix A, and the strictly upper
+*> triangular part of A is not referenced.
+*>
+*> On exit, if INFO = 0, the transformed matrix, stored in the
+*> same format as A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is COMPLEX*16 array, dimension (LDB,N)
+*> The triangular factor from the Cholesky factorization of B,
+*> as returned by ZPOTRF.
+*> B is modified by the routine but restored on exit.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit.
+*> < 0: if INFO = -i, the i-th argument had an illegal value.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date December 2016
+*
+*> \ingroup complex16HEcomputational
+*
+* =====================================================================
+ SUBROUTINE ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
+*
+* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* December 2016
+*
+* .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, ITYPE, LDA, LDB, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), B( LDB, * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ONE, HALF
+ PARAMETER ( ONE = 1.0D+0, HALF = 0.5D+0 )
+ COMPLEX*16 CONE
+ PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER K
+ DOUBLE PRECISION AKK, BKK
+ COMPLEX*16 CT
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZAXPY, ZDSCAL, ZHER2, ZLACGV, ZTRMV,
+ $ ZTRSV
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
+ INFO = -1
+ ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -5
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZHEGS2', -INFO )
+ RETURN
+ END IF
+*
+ IF( ITYPE.EQ.1 ) THEN
+ IF( UPPER ) THEN
+*
+* Compute inv(U**H)*A*inv(U)
+*
+ DO 10 K = 1, N
+*
+* Update the upper triangle of A(k:n,k:n)
+*
+ AKK = A( K, K )
+ BKK = B( K, K )
+ AKK = AKK / BKK**2
+ A( K, K ) = AKK
+ IF( K.LT.N ) THEN
+ CALL ZDSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA )
+ CT = -HALF*AKK
+ CALL ZLACGV( N-K, A( K, K+1 ), LDA )
+ CALL ZLACGV( N-K, B( K, K+1 ), LDB )
+ CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
+ $ LDA )
+ CALL ZHER2( UPLO, N-K, -CONE, A( K, K+1 ), LDA,
+ $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA )
+ CALL ZAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ),
+ $ LDA )
+ CALL ZLACGV( N-K, B( K, K+1 ), LDB )
+ CALL ZTRSV( UPLO, 'Conjugate transpose', 'Non-unit',
+ $ N-K, B( K+1, K+1 ), LDB, A( K, K+1 ),
+ $ LDA )
+ CALL ZLACGV( N-K, A( K, K+1 ), LDA )
+ END IF
+ 10 CONTINUE
+ ELSE
+*
+* Compute inv(L)*A*inv(L**H)
+*
+ DO 20 K = 1, N
+*
+* Update the lower triangle of A(k:n,k:n)
+*
+ AKK = A( K, K )
+ BKK = B( K, K )
+ AKK = AKK / BKK**2
+ A( K, K ) = AKK
+ IF( K.LT.N ) THEN
+ CALL ZDSCAL( N-K, ONE / BKK, A( K+1, K ), 1 )
+ CT = -HALF*AKK
+ CALL ZAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
+ CALL ZHER2( UPLO, N-K, -CONE, A( K+1, K ), 1,
+ $ B( K+1, K ), 1, A( K+1, K+1 ), LDA )
+ CALL ZAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 )
+ CALL ZTRSV( UPLO, 'No transpose', 'Non-unit', N-K,
+ $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 )
+ END IF
+ 20 CONTINUE
+ END IF
+ ELSE
+ IF( UPPER ) THEN
+*
+* Compute U*A*U**H
+*
+ DO 30 K = 1, N
+*
+* Update the upper triangle of A(1:k,1:k)
+*
+ AKK = A( K, K )
+ BKK = B( K, K )
+ CALL ZTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B,
+ $ LDB, A( 1, K ), 1 )
+ CT = HALF*AKK
+ CALL ZAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
+ CALL ZHER2( UPLO, K-1, CONE, A( 1, K ), 1, B( 1, K ), 1,
+ $ A, LDA )
+ CALL ZAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 )
+ CALL ZDSCAL( K-1, BKK, A( 1, K ), 1 )
+ A( K, K ) = AKK*BKK**2
+ 30 CONTINUE
+ ELSE
+*
+* Compute L**H *A*L
+*
+ DO 40 K = 1, N
+*
+* Update the lower triangle of A(1:k,1:k)
+*
+ AKK = A( K, K )
+ BKK = B( K, K )
+ CALL ZLACGV( K-1, A( K, 1 ), LDA )
+ CALL ZTRMV( UPLO, 'Conjugate transpose', 'Non-unit', K-1,
+ $ B, LDB, A( K, 1 ), LDA )
+ CT = HALF*AKK
+ CALL ZLACGV( K-1, B( K, 1 ), LDB )
+ CALL ZAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
+ CALL ZHER2( UPLO, K-1, CONE, A( K, 1 ), LDA, B( K, 1 ),
+ $ LDB, A, LDA )
+ CALL ZAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA )
+ CALL ZLACGV( K-1, B( K, 1 ), LDB )
+ CALL ZDSCAL( K-1, BKK, A( K, 1 ), LDA )
+ CALL ZLACGV( K-1, A( K, 1 ), LDA )
+ A( K, K ) = AKK*BKK**2
+ 40 CONTINUE
+ END IF
+ END IF
+ RETURN
+*
+* End of ZHEGS2
+*
+ END
diff --git a/fortran/zhegst.f b/fortran/zhegst.f
new file mode 100644
index 0000000..ec9fa16
--- /dev/null
+++ b/fortran/zhegst.f
@@ -0,0 +1,332 @@
+*> \brief \b ZHEGST
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZHEGST + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER INFO, ITYPE, LDA, LDB, N
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A( LDA, * ), B( LDB, * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZHEGST reduces a complex Hermitian-definite generalized
+*> eigenproblem to standard form.
+*>
+*> If ITYPE = 1, the problem is A*x = lambda*B*x,
+*> and A is overwritten by inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H)
+*>
+*> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
+*> B*A*x = lambda*x, and A is overwritten by U*A*U**H or L**H*A*L.
+*>
+*> B must have been previously factorized as U**H*U or L*L**H by ZPOTRF.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] ITYPE
+*> \verbatim
+*> ITYPE is INTEGER
+*> = 1: compute inv(U**H)*A*inv(U) or inv(L)*A*inv(L**H);
+*> = 2 or 3: compute U*A*U**H or L**H*A*L.
+*> \endverbatim
+*>
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> = 'U': Upper triangle of A is stored and B is factored as
+*> U**H*U;
+*> = 'L': Lower triangle of A is stored and B is factored as
+*> L*L**H.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A and B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA,N)
+*> On entry, the Hermitian matrix A. If UPLO = 'U', the leading
+*> N-by-N upper triangular part of A contains the upper
+*> triangular part of the matrix A, and the strictly lower
+*> triangular part of A is not referenced. If UPLO = 'L', the
+*> leading N-by-N lower triangular part of A contains the lower
+*> triangular part of the matrix A, and the strictly upper
+*> triangular part of A is not referenced.
+*>
+*> On exit, if INFO = 0, the transformed matrix, stored in the
+*> same format as A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is COMPLEX*16 array, dimension (LDB,N)
+*> The triangular factor from the Cholesky factorization of B,
+*> as returned by ZPOTRF.
+*> B is modified by the routine but restored on exit.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date December 2016
+*
+*> \ingroup complex16HEcomputational
+*
+* =====================================================================
+ SUBROUTINE ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
+*
+* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* December 2016
+*
+* .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, ITYPE, LDA, LDB, N
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), B( LDB, * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ DOUBLE PRECISION ONE
+ PARAMETER ( ONE = 1.0D+0 )
+ COMPLEX*16 CONE, HALF
+ PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
+ $ HALF = ( 0.5D+0, 0.0D+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER K, KB, NB
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZHEGS2, ZHEMM, ZHER2K, ZTRMM, ZTRSM
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER ILAENV
+ EXTERNAL LSAME, ILAENV
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
+ INFO = -1
+ ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -5
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZHEGST', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Determine the block size for this environment.
+*
+ NB = ILAENV( 1, 'ZHEGST', UPLO, N, -1, -1, -1 )
+*
+ IF( NB.LE.1 .OR. NB.GE.N ) THEN
+*
+* Use unblocked code
+*
+ CALL ZHEGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
+ ELSE
+*
+* Use blocked code
+*
+ IF( ITYPE.EQ.1 ) THEN
+ IF( UPPER ) THEN
+*
+* Compute inv(U**H)*A*inv(U)
+*
+ DO 10 K = 1, N, NB
+ KB = MIN( N-K+1, NB )
+*
+* Update the upper triangle of A(k:n,k:n)
+*
+ CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
+ $ B( K, K ), LDB, INFO )
+ IF( K+KB.LE.N ) THEN
+ CALL ZTRSM( 'L', UPLO, 'C',
+ $ 'N', KB, N-K-KB+1, CONE,
+ $ B( K, K ), LDB, A( K, K+KB ), LDA )
+ CALL ZHEMM( 'L', UPLO, KB, N-K-KB+1, -HALF,
+ $ A( K, K ), LDA, B( K, K+KB ), LDB,
+ $ CONE, A( K, K+KB ), LDA )
+ CALL ZHER2K( UPLO, 'C', N-K-KB+1,
+ $ KB, -CONE, A( K, K+KB ), LDA,
+ $ B( K, K+KB ), LDB, ONE,
+ $ A( K+KB, K+KB ), LDA )
+ CALL ZHEMM( 'L', UPLO, KB, N-K-KB+1, -HALF,
+ $ A( K, K ), LDA, B( K, K+KB ), LDB,
+ $ CONE, A( K, K+KB ), LDA )
+ CALL ZTRSM( 'R', UPLO, 'N',
+ $ 'N', KB, N-K-KB+1, CONE,
+ $ B( K+KB, K+KB ), LDB, A( K, K+KB ),
+ $ LDA )
+ END IF
+ 10 CONTINUE
+ ELSE
+*
+* Compute inv(L)*A*inv(L**H)
+*
+ DO 20 K = 1, N, NB
+ KB = MIN( N-K+1, NB )
+*
+* Update the lower triangle of A(k:n,k:n)
+*
+ CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
+ $ B( K, K ), LDB, INFO )
+ IF( K+KB.LE.N ) THEN
+ CALL ZTRSM( 'R', UPLO, 'C',
+ $ 'N', N-K-KB+1, KB, CONE,
+ $ B( K, K ), LDB, A( K+KB, K ), LDA )
+ CALL ZHEMM( 'R', UPLO, N-K-KB+1, KB, -HALF,
+ $ A( K, K ), LDA, B( K+KB, K ), LDB,
+ $ CONE, A( K+KB, K ), LDA )
+ CALL ZHER2K( UPLO, 'N', N-K-KB+1, KB,
+ $ -CONE, A( K+KB, K ), LDA,
+ $ B( K+KB, K ), LDB, ONE,
+ $ A( K+KB, K+KB ), LDA )
+ CALL ZHEMM( 'R', UPLO, N-K-KB+1, KB, -HALF,
+ $ A( K, K ), LDA, B( K+KB, K ), LDB,
+ $ CONE, A( K+KB, K ), LDA )
+ CALL ZTRSM( 'L', UPLO, 'N',
+ $ 'N', N-K-KB+1, KB, CONE,
+ $ B( K+KB, K+KB ), LDB, A( K+KB, K ),
+ $ LDA )
+ END IF
+ 20 CONTINUE
+ END IF
+ ELSE
+ IF( UPPER ) THEN
+*
+* Compute U*A*U**H
+*
+ DO 30 K = 1, N, NB
+ KB = MIN( N-K+1, NB )
+*
+* Update the upper triangle of A(1:k+kb-1,1:k+kb-1)
+*
+ CALL ZTRMM( 'L', UPLO, 'N', 'N',
+ $ K-1, KB, CONE, B, LDB, A( 1, K ), LDA )
+ CALL ZHEMM( 'R', UPLO, K-1, KB, HALF, A( K, K ),
+ $ LDA, B( 1, K ), LDB, CONE, A( 1, K ),
+ $ LDA )
+ CALL ZHER2K( UPLO, 'N', K-1, KB, CONE,
+ $ A( 1, K ), LDA, B( 1, K ), LDB, ONE, A,
+ $ LDA )
+ CALL ZHEMM( 'R', UPLO, K-1, KB, HALF, A( K, K ),
+ $ LDA, B( 1, K ), LDB, CONE, A( 1, K ),
+ $ LDA )
+ CALL ZTRMM( 'R', UPLO, 'C',
+ $ 'N', K-1, KB, CONE, B( K, K ), LDB,
+ $ A( 1, K ), LDA )
+ CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
+ $ B( K, K ), LDB, INFO )
+ 30 CONTINUE
+ ELSE
+*
+* Compute L**H*A*L
+*
+ DO 40 K = 1, N, NB
+ KB = MIN( N-K+1, NB )
+*
+* Update the lower triangle of A(1:k+kb-1,1:k+kb-1)
+*
+ CALL ZTRMM( 'R', UPLO, 'N', 'N',
+ $ KB, K-1, CONE, B, LDB, A( K, 1 ), LDA )
+ CALL ZHEMM( 'L', UPLO, KB, K-1, HALF, A( K, K ),
+ $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
+ $ LDA )
+ CALL ZHER2K( UPLO, 'C', K-1, KB,
+ $ CONE, A( K, 1 ), LDA, B( K, 1 ), LDB,
+ $ ONE, A, LDA )
+ CALL ZHEMM( 'L', UPLO, KB, K-1, HALF, A( K, K ),
+ $ LDA, B( K, 1 ), LDB, CONE, A( K, 1 ),
+ $ LDA )
+ CALL ZTRMM( 'L', UPLO, 'C',
+ $ 'N', KB, K-1, CONE, B( K, K ), LDB,
+ $ A( K, 1 ), LDA )
+ CALL ZHEGS2( ITYPE, UPLO, KB, A( K, K ), LDA,
+ $ B( K, K ), LDB, INFO )
+ 40 CONTINUE
+ END IF
+ END IF
+ END IF
+ RETURN
+*
+* End of ZHEGST
+*
+ END
diff --git a/fortran/zhegv.f b/fortran/zhegv.f
new file mode 100644
index 0000000..761b5d0
--- /dev/null
+++ b/fortran/zhegv.f
@@ -0,0 +1,321 @@
+*> \brief \b ZHEGV
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download ZHEGV + dependencies
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZHEGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
+* LWORK, RWORK, INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER JOBZ, UPLO
+* INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
+* ..
+* .. Array Arguments ..
+* DOUBLE PRECISION RWORK( * ), W( * )
+* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZHEGV computes all the eigenvalues, and optionally, the eigenvectors
+*> of a complex generalized Hermitian-definite eigenproblem, of the form
+*> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
+*> Here A and B are assumed to be Hermitian and B is also
+*> positive definite.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] ITYPE
+*> \verbatim
+*> ITYPE is INTEGER
+*> Specifies the problem type to be solved:
+*> = 1: A*x = (lambda)*B*x
+*> = 2: A*B*x = (lambda)*x
+*> = 3: B*A*x = (lambda)*x
+*> \endverbatim
+*>
+*> \param[in] JOBZ
+*> \verbatim
+*> JOBZ is CHARACTER*1
+*> = 'N': Compute eigenvalues only;
+*> = 'V': Compute eigenvalues and eigenvectors.
+*> \endverbatim
+*>
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> = 'U': Upper triangles of A and B are stored;
+*> = 'L': Lower triangles of A and B are stored.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrices A and B. N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*> A is COMPLEX*16 array, dimension (LDA, N)
+*> On entry, the Hermitian matrix A. If UPLO = 'U', the
+*> leading N-by-N upper triangular part of A contains the
+*> upper triangular part of the matrix A. If UPLO = 'L',
+*> the leading N-by-N lower triangular part of A contains
+*> the lower triangular part of the matrix A.
+*>
+*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
+*> matrix Z of eigenvectors. The eigenvectors are normalized
+*> as follows:
+*> if ITYPE = 1 or 2, Z**H*B*Z = I;
+*> if ITYPE = 3, Z**H*inv(B)*Z = I.
+*> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
+*> or the lower triangle (if UPLO='L') of A, including the
+*> diagonal, is destroyed.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is COMPLEX*16 array, dimension (LDB, N)
+*> On entry, the Hermitian positive definite matrix B.
+*> If UPLO = 'U', the leading N-by-N upper triangular part of B
+*> contains the upper triangular part of the matrix B.
+*> If UPLO = 'L', the leading N-by-N lower triangular part of B
+*> contains the lower triangular part of the matrix B.
+*>
+*> On exit, if INFO <= N, the part of B containing the matrix is
+*> overwritten by the triangular factor U or L from the Cholesky
+*> factorization B = U**H*U or B = L*L**H.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] W
+*> \verbatim
+*> W is DOUBLE PRECISION array, dimension (N)
+*> If INFO = 0, the eigenvalues in ascending order.
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
+*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+*> \endverbatim
+*>
+*> \param[in] LWORK
+*> \verbatim
+*> LWORK is INTEGER
+*> The length of the array WORK. LWORK >= max(1,2*N-1).
+*> For optimal efficiency, LWORK >= (NB+1)*N,
+*> where NB is the blocksize for ZHETRD 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] RWORK
+*> \verbatim
+*> RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> > 0: ZPOTRF or ZHEEV returned an error code:
+*> <= N: if INFO = i, ZHEEV failed to converge;
+*> i off-diagonal elements of an intermediate
+*> tridiagonal form did not converge to zero;
+*> > N: if INFO = N + i, for 1 <= i <= N, then the leading
+*> minor of order i of B is not positive definite.
+*> The factorization of B could not be completed and
+*> no eigenvalues or eigenvectors were computed.
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date December 2016
+*
+*> \ingroup complex16HEeigen
+*
+* =====================================================================
+ SUBROUTINE ZHEGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK,
+ $ LWORK, RWORK, INFO )
+*
+* -- LAPACK driver routine (version 3.7.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* December 2016
+*
+* .. Scalar Arguments ..
+ CHARACTER JOBZ, UPLO
+ INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
+* ..
+* .. Array Arguments ..
+ DOUBLE PRECISION RWORK( * ), W( * )
+ COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 ONE
+ PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL LQUERY, UPPER, WANTZ
+ CHARACTER TRANS
+ INTEGER LWKOPT, NB, NEIG
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ INTEGER ILAENV
+ EXTERNAL LSAME, ILAENV
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA, ZHEEV, ZHEGST, ZPOTRF, ZTRMM, ZTRSM
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC MAX
+* ..
+* .. Executable Statements ..
+*
+* Test the input parameters.
+*
+ WANTZ = LSAME( JOBZ, 'V' )
+ UPPER = LSAME( UPLO, 'U' )
+ LQUERY = ( LWORK.EQ.-1 )
+*
+ INFO = 0
+ IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
+ INFO = -1
+ ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
+ INFO = -2
+ ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -8
+ END IF
+*
+ IF( INFO.EQ.0 ) THEN
+ NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
+ LWKOPT = MAX( 1, ( NB + 1 )*N )
+ WORK( 1 ) = LWKOPT
+*
+ IF( LWORK.LT.MAX( 1, 2*N - 1 ) .AND. .NOT.LQUERY ) THEN
+ INFO = -11
+ END IF
+ END IF
+*
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZHEGV ', -INFO )
+ RETURN
+ ELSE IF( LQUERY ) THEN
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 )
+ $ RETURN
+*
+* Form a Cholesky factorization of B.
+*
+ CALL ZPOTRF( UPLO, N, B, LDB, INFO )
+ IF( INFO.NE.0 ) THEN
+ INFO = N + INFO
+ RETURN
+ END IF
+*
+* Transform problem to standard eigenvalue problem and solve.
+*
+ CALL ZHEGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
+ CALL ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO )
+*
+ IF( WANTZ ) THEN
+*
+* Backtransform eigenvectors to the original problem.
+*
+ NEIG = N
+ IF( INFO.GT.0 )
+ $ NEIG = INFO - 1
+ IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
+*
+* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
+* backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
+*
+ IF( UPPER ) THEN
+ TRANS = 'N'
+ ELSE
+ TRANS = 'C'
+ END IF
+*
+ CALL ZTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
+ $ B, LDB, A, LDA )
+*
+ ELSE IF( ITYPE.EQ.3 ) THEN
+*
+* For B*A*x=(lambda)*x;
+* backtransform eigenvectors: x = L*y or U**H *y
+*
+ IF( UPPER ) THEN
+ TRANS = 'C'
+ ELSE
+ TRANS = 'N'
+ END IF
+*
+ CALL ZTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, NEIG, ONE,
+ $ B, LDB, A, LDA )
+ END IF
+ END IF
+*
+ WORK( 1 ) = LWKOPT
+*
+ RETURN
+*
+* End of ZHEGV
+*
+ END
diff --git a/fortran/zhemm.f b/fortran/zhemm.f
new file mode 100644
index 0000000..d63778b
--- /dev/null
+++ b/fortran/zhemm.f
@@ -0,0 +1,371 @@
+*> \brief \b ZHEMM
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+*
+* .. Scalar Arguments ..
+* COMPLEX*16 ALPHA,BETA
+* INTEGER LDA,LDB,LDC,M,N
+* CHARACTER SIDE,UPLO
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZHEMM 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 an hermitian 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 hermitian 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 hermitian matrix A is to be
+*> referenced as follows:
+*>
+*> UPLO = 'U' or 'u' Only the upper triangular part of the
+*> hermitian matrix is to be referenced.
+*>
+*> UPLO = 'L' or 'l' Only the lower triangular part of the
+*> hermitian 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 COMPLEX*16
+*> On entry, ALPHA specifies the scalar alpha.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*> A is COMPLEX*16 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 hermitian 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 hermitian 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 hermitian
+*> 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 hermitian 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 hermitian 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 hermitian
+*> matrix and the strictly upper triangular part of A is not
+*> referenced.
+*> Note that the imaginary parts of the diagonal elements need
+*> not be set, they are assumed to be zero.
+*> \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 COMPLEX*16 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 COMPLEX*16
+*> 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 COMPLEX*16 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.
+*
+*> \date December 2016
+*
+*> \ingroup complex16_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 ZHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+*
+* -- Reference BLAS level3 routine (version 3.7.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* December 2016
+*
+* .. Scalar Arguments ..
+ COMPLEX*16 ALPHA,BETA
+ INTEGER LDA,LDB,LDC,M,N
+ CHARACTER SIDE,UPLO
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*)
+* ..
+*
+* =====================================================================
+*
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DBLE,DCONJG,MAX
+* ..
+* .. Local Scalars ..
+ COMPLEX*16 TEMP1,TEMP2
+ INTEGER I,INFO,J,K,NROWA
+ LOGICAL UPPER
+* ..
+* .. Parameters ..
+ COMPLEX*16 ONE
+ PARAMETER (ONE= (1.0D+0,0.0D+0))
+ COMPLEX*16 ZERO
+ PARAMETER (ZERO= (0.0D+0,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('ZHEMM ',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)*DCONJG(A(K,I))
+ 50 CONTINUE
+ IF (BETA.EQ.ZERO) THEN
+ C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2
+ ELSE
+ C(I,J) = BETA*C(I,J) + TEMP1*DBLE(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)*DCONJG(A(K,I))
+ 80 CONTINUE
+ IF (BETA.EQ.ZERO) THEN
+ C(I,J) = TEMP1*DBLE(A(I,I)) + ALPHA*TEMP2
+ ELSE
+ C(I,J) = BETA*C(I,J) + TEMP1*DBLE(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*DBLE(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*DCONJG(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*DCONJG(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 ZHEMM .
+*
+ END
diff --git a/fortran/ztrsv.f b/fortran/ztrsv.f
new file mode 100644
index 0000000..ba7aa35
--- /dev/null
+++ b/fortran/ztrsv.f
@@ -0,0 +1,375 @@
+*> \brief \b ZTRSV
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE ZTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+*
+* .. Scalar Arguments ..
+* INTEGER INCX,LDA,N
+* CHARACTER DIAG,TRANS,UPLO
+* ..
+* .. Array Arguments ..
+* COMPLEX*16 A(LDA,*),X(*)
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*>
+*> ZTRSV 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.
+*>
+*> 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] A
+*> \verbatim
+*> A is COMPLEX*16 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 COMPLEX*16 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.
+*
+*> \date December 2016
+*
+*> \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 ZTRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+*
+* -- Reference BLAS level2 routine (version 3.7.0) --
+* -- Reference BLAS is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* December 2016
+*
+* .. Scalar Arguments ..
+ INTEGER INCX,LDA,N
+ CHARACTER DIAG,TRANS,UPLO
+* ..
+* .. Array Arguments ..
+ COMPLEX*16 A(LDA,*),X(*)
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX*16 ZERO
+ PARAMETER (ZERO= (0.0D+0,0.0D+0))
+* ..
+* .. Local Scalars ..
+ COMPLEX*16 TEMP
+ INTEGER I,INFO,IX,J,JX,KX
+ LOGICAL NOCONJ,NOUNIT
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC DCONJG,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('ZTRSV ',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 A are
+* accessed sequentially with one pass through A.
+*
+ IF (LSAME(TRANS,'N')) THEN
+*
+* Form x := inv( A )*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ IF (INCX.EQ.1) THEN
+ DO 20 J = N,1,-1
+ IF (X(J).NE.ZERO) THEN
+ IF (NOUNIT) X(J) = X(J)/A(J,J)
+ TEMP = X(J)
+ DO 10 I = J - 1,1,-1
+ X(I) = X(I) - TEMP*A(I,J)
+ 10 CONTINUE
+ END IF
+ 20 CONTINUE
+ ELSE
+ JX = KX + (N-1)*INCX
+ DO 40 J = N,1,-1
+ IF (X(JX).NE.ZERO) THEN
+ IF (NOUNIT) X(JX) = X(JX)/A(J,J)
+ TEMP = X(JX)
+ IX = JX
+ DO 30 I = J - 1,1,-1
+ IX = IX - INCX
+ X(IX) = X(IX) - TEMP*A(I,J)
+ 30 CONTINUE
+ END IF
+ JX = JX - INCX
+ 40 CONTINUE
+ END IF
+ ELSE
+ IF (INCX.EQ.1) THEN
+ DO 60 J = 1,N
+ IF (X(J).NE.ZERO) THEN
+ IF (NOUNIT) X(J) = X(J)/A(J,J)
+ TEMP = X(J)
+ DO 50 I = J + 1,N
+ X(I) = X(I) - TEMP*A(I,J)
+ 50 CONTINUE
+ END IF
+ 60 CONTINUE
+ ELSE
+ JX = KX
+ DO 80 J = 1,N
+ IF (X(JX).NE.ZERO) THEN
+ IF (NOUNIT) X(JX) = X(JX)/A(J,J)
+ TEMP = X(JX)
+ IX = JX
+ DO 70 I = J + 1,N
+ IX = IX + INCX
+ X(IX) = X(IX) - TEMP*A(I,J)
+ 70 CONTINUE
+ END IF
+ JX = JX + INCX
+ 80 CONTINUE
+ END IF
+ END IF
+ ELSE
+*
+* Form x := inv( A**T )*x or x := inv( A**H )*x.
+*
+ IF (LSAME(UPLO,'U')) THEN
+ IF (INCX.EQ.1) THEN
+ DO 110 J = 1,N
+ TEMP = X(J)
+ IF (NOCONJ) THEN
+ DO 90 I = 1,J - 1
+ TEMP = TEMP - A(I,J)*X(I)
+ 90 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/A(J,J)
+ ELSE
+ DO 100 I = 1,J - 1
+ TEMP = TEMP - DCONJG(A(I,J))*X(I)
+ 100 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
+ END IF
+ X(J) = TEMP
+ 110 CONTINUE
+ ELSE
+ JX = KX
+ DO 140 J = 1,N
+ IX = KX
+ TEMP = X(JX)
+ IF (NOCONJ) THEN
+ DO 120 I = 1,J - 1
+ TEMP = TEMP - A(I,J)*X(IX)
+ IX = IX + INCX
+ 120 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/A(J,J)
+ ELSE
+ DO 130 I = 1,J - 1
+ TEMP = TEMP - DCONJG(A(I,J))*X(IX)
+ IX = IX + INCX
+ 130 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
+ END IF
+ X(JX) = TEMP
+ JX = JX + INCX
+ 140 CONTINUE
+ END IF
+ ELSE
+ IF (INCX.EQ.1) THEN
+ DO 170 J = N,1,-1
+ TEMP = X(J)
+ IF (NOCONJ) THEN
+ DO 150 I = N,J + 1,-1
+ TEMP = TEMP - A(I,J)*X(I)
+ 150 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/A(J,J)
+ ELSE
+ DO 160 I = N,J + 1,-1
+ TEMP = TEMP - DCONJG(A(I,J))*X(I)
+ 160 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
+ END IF
+ X(J) = TEMP
+ 170 CONTINUE
+ ELSE
+ KX = KX + (N-1)*INCX
+ JX = KX
+ DO 200 J = N,1,-1
+ IX = KX
+ TEMP = X(JX)
+ IF (NOCONJ) THEN
+ DO 180 I = N,J + 1,-1
+ TEMP = TEMP - A(I,J)*X(IX)
+ IX = IX - INCX
+ 180 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/A(J,J)
+ ELSE
+ DO 190 I = N,J + 1,-1
+ TEMP = TEMP - DCONJG(A(I,J))*X(IX)
+ IX = IX - INCX
+ 190 CONTINUE
+ IF (NOUNIT) TEMP = TEMP/DCONJG(A(J,J))
+ END IF
+ X(JX) = TEMP
+ JX = JX - INCX
+ 200 CONTINUE
+ END IF
+ END IF
+ END IF
+*
+ RETURN
+*
+* End of ZTRSV .
+*
+ END