Add SYTRS2 routine - A BLAS 3 version of SYTRS

Add SYCONV routine: convert back and forth the factorization returned by SYTRF to be able to call SYTRS2.
Modify SYSV that now is calling SYTRS2 instead of SYTRS (and also SYCONV to convert and revert the factorization returned by SYTRF).
Modify testing to have TRS but also TRS2 tested in the LIN testing for SY.

13 files changed, 2348 insertions, 17 deletions

diff --git a/SRC/Makefile b/SRC/Makefile
index ce43806d..67315be6 100644
--- a/SRC/Makefile
+++ b/SRC/Makefile
@@ -131,7 +131,7 @@ SLASRC = \
    ssptrf.o ssptri.o ssptrs.o sstegr.o sstein.o sstev.o  sstevd.o sstevr.o \
    sstevx.o ssycon.o ssyev.o  ssyevd.o ssyevr.o ssyevx.o ssygs2.o \
    ssygst.o ssygv.o  ssygvd.o ssygvx.o ssyrfs.o ssysv.o  ssysvx.o \
-   ssytd2.o ssytf2.o ssytrd.o ssytrf.o ssytri.o ssytrs.o stbcon.o \
+   ssytd2.o ssytf2.o ssytrd.o ssytrf.o ssytri.o ssytrs.o ssytrs2.o ssyconv.o stbcon.o \
    stbrfs.o stbtrs.o stgevc.o stgex2.o stgexc.o stgsen.o \
    stgsja.o stgsna.o stgsy2.o stgsyl.o stpcon.o stprfs.o stptri.o \
    stptrs.o \
@@ -195,7 +195,7 @@ CLASRC = \
    cspsvx.o csptrf.o csptri.o csptrs.o csrscl.o cstedc.o \
    cstegr.o cstein.o csteqr.o csycon.o csymv.o \
    csyr.o   csyrfs.o csysv.o  csysvx.o csytf2.o csytrf.o csytri.o \
-   csytrs.o ctbcon.o ctbrfs.o ctbtrs.o ctgevc.o ctgex2.o \
+   csytrs.o csytrs2.o csyconv.o ctbcon.o ctbrfs.o ctbtrs.o ctgevc.o ctgex2.o \
    ctgexc.o ctgsen.o ctgsja.o ctgsna.o ctgsy2.o ctgsyl.o ctpcon.o \
    ctprfs.o ctptri.o \
    ctptrs.o ctrcon.o ctrevc.o ctrexc.o ctrrfs.o ctrsen.o ctrsna.o \
@@ -264,7 +264,7 @@ DLASRC = \
    dstevx.o dsycon.o dsyev.o  dsyevd.o dsyevr.o \
    dsyevx.o dsygs2.o dsygst.o dsygv.o  dsygvd.o dsygvx.o dsyrfs.o \
    dsysv.o  dsysvx.o \
-   dsytd2.o dsytf2.o dsytrd.o dsytrf.o dsytri.o dsytrs.o dtbcon.o \
+   dsytd2.o dsytf2.o dsytrd.o dsytrf.o dsytri.o dsytrs.o dsytrs2.o dsyconv.o dtbcon.o \
    dtbrfs.o dtbtrs.o dtgevc.o dtgex2.o dtgexc.o dtgsen.o \
    dtgsja.o dtgsna.o dtgsy2.o dtgsyl.o dtpcon.o dtprfs.o dtptri.o \
    dtptrs.o \
@@ -332,7 +332,7 @@ ZLASRC = \
    zspsvx.o zsptrf.o zsptri.o zsptrs.o zdrscl.o zstedc.o \
    zstegr.o zstein.o zsteqr.o zsycon.o zsymv.o \
    zsyr.o   zsyrfs.o zsysv.o  zsysvx.o zsytf2.o zsytrf.o zsytri.o \
-   zsytrs.o ztbcon.o ztbrfs.o ztbtrs.o ztgevc.o ztgex2.o \
+   zsytrs.o zsytrs2.o zsyconv.o ztbcon.o ztbrfs.o ztbtrs.o ztgevc.o ztgex2.o \
    ztgexc.o ztgsen.o ztgsja.o ztgsna.o ztgsy2.o ztgsyl.o ztpcon.o \
    ztprfs.o ztptri.o \
    ztptrs.o ztrcon.o ztrevc.o ztrexc.o ztrrfs.o ztrsen.o ztrsna.o \
diff --git a/SRC/csyconv.f b/SRC/csyconv.f
new file mode 100644
index 00000000..74306022
--- /dev/null
+++ b/SRC/csyconv.f
@@ -0,0 +1,302 @@
+      SUBROUTINE CSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO, WAY
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CSYCONV convert A given by TRF into L and D and vice-versa.
+*  Get Non-diag elements of D (returned in workspace) and 
+*  apply or reverse permutation done in TRF.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+* 
+*  WAY     (input) CHARACTER*1
+*          = 'C': Convert 
+*          = 'R': Revert
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input) COMPLEX array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by CSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by CSYTRF.
+*
+* WORK     (workspace) COMPLEX array, dimension (N)
+*
+* LWORK    (input) INTEGER
+*          The length of WORK.  LWORK >=1. 
+*          LWORK = N
+*
+*          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.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            ZERO
+      PARAMETER          ( ZERO = (0.0E+0,0.0E+0) )
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Local Scalars ..
+      LOGICAL            UPPER, CONVERT
+      INTEGER            I, IP
+      COMPLEX            TEMP
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      CONVERT = LSAME( WAY, 'C' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN
+         INFO = -2
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYCONV', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*      A is UPPER
+*
+*      Convert A (A is upper)
+*
+*        Convert VALUE
+*
+         IF ( CONVERT ) THEN
+            I=N
+            WORK(1)=ZERO
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I-1,I)
+                  A(I-1,I)=ZERO
+                  I=I-1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Convert PERMUTATIONS
+*  
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0) THEN
+               IP=IPIV(I)
+               IF( I .LT. N) THEN
+                  DO 12 J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+ 12            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+               IF( I .LT. N) THEN
+             DO 13 J= I+1,N
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I-1,J)
+                 A(I-1,J)=TEMP
+ 13            CONTINUE
+                ENDIF
+                I=I-1
+           ENDIF
+           I=I-1
+        END DO
+
+         ELSE
+*
+*      Revert A (A is upper)
+*
+*
+*        Revert PERMUTATIONS
+*  
+            I=1
+            DO WHILE ( I .LE. N )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF( I .LT. N) THEN
+                  DO J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+                  END DO
+                  ENDIF
+               ELSE
+                 IP=-IPIV(I)
+                 I=I+1
+                 IF( I .LT. N) THEN
+                    DO J= I+1,N
+                       TEMP=A(IP,J)
+                       A(IP,J)=A(I-1,J)
+                       A(I-1,J)=TEMP
+                    END DO
+                 ENDIF
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Revert VALUE
+*
+            I=N
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  A(I-1,I)=WORK(I)
+                  I=I-1
+               ENDIF
+               I=I-1
+            END DO
+         END IF
+      ELSE
+*
+*      A is LOWER
+*
+         IF ( CONVERT ) THEN
+*
+*      Convert A (A is lower)
+*
+*
+*        Convert VALUE
+*
+            I=1
+            WORK(N)=ZERO
+            DO WHILE ( I .LE. N )
+               IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I+1,I)
+                  A(I+1,I)=ZERO
+                  I=I+1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Convert PERMUTATIONS
+*
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+               IP=IPIV(I)
+               IF (I .GT. 1) THEN
+               DO 22 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I,J)
+                 A(I,J)=TEMP
+ 22            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+              IF ( I .GT. 1 ) THEN
+              DO 23 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I+1,J)
+                 A(I+1,J)=TEMP
+ 23           CONTINUE
+              ENDIF
+              I=I+1
+           ENDIF
+           I=I+1
+        END DO
+         ELSE
+*
+*      Revert A (A is lower)
+*
+*
+*        Revert PERMUTATIONS
+*
+            I=N
+            DO WHILE ( I .GE. 1 )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I,J)
+                        A(I,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ELSE
+                  IP=-IPIV(I)
+                  I=I-1
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I+1,J)
+                        A(I+1,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Revert VALUE
+*
+            I=1
+            DO WHILE ( I .LE. N-1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  A(I+1,I)=WORK(I)
+                  I=I+1
+               ENDIF
+               I=I+1
+            END DO
+         END IF
+      END IF
+
+      RETURN
+*
+*     End of CSYCONV
+*
+      END
diff --git a/SRC/csysv.f b/SRC/csysv.f
index 855405f2..a9832e53 100644
--- a/SRC/csysv.f
+++ b/SRC/csysv.f
@@ -113,7 +113,7 @@
       EXTERNAL           ILAENV, LSAME
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           CSYTRF, CSYTRS, XERBLA
+      EXTERNAL           CSYCONV, CSYTRF, CSYTRS2, XERBLA
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MAX
@@ -160,9 +160,17 @@
       CALL CSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN
 *
+*        Convert A
+*
+         CALL CSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
 *        Solve the system A*X = B, overwriting B with X.
 *
-         CALL CSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         CALL CSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO )
+*
+*        Revert A
+*
+         CALL CSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
 *
       END IF
 *
diff --git a/SRC/csytrs2.f b/SRC/csytrs2.f
new file mode 100644
index 00000000..bd670c1b
--- /dev/null
+++ b/SRC/csytrs2.f
@@ -0,0 +1,272 @@
+      SUBROUTINE CSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
+     $                    WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  CSYTRS2 solves a system of linear equations A*X = B with a COMPLEX
+*  symmetric matrix A using the factorization A = U*D*U**T or
+*  A = L*D*L**T computed by CSYTRF and converted by CSYCONV.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  NRHS    (input) INTEGER
+*          The number of right hand sides, i.e., the number of columns
+*          of the matrix B.  NRHS >= 0.
+*
+*  A       (input) COMPLEX array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by CSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by CSYTRF.
+*
+*  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
+*          On entry, the right hand side matrix B.
+*          On exit, the solution matrix X.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,N).
+*
+*  WORK    (workspace) COMPLEX array, dimension (N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            ONE
+      PARAMETER          ( ONE = (1.0E+0,0.0E+0) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, J, K, KP
+      COMPLEX            AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CSCAL, CSWAP, CTRSM, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+      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
+      ELSE IF( NRHS.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'CSYTRS2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. NRHS.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Solve A*X = B, where A = U*D*U'.
+*
+*       P' * B  
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( KP.EQ.-IPIV( K-1 ) )
+     $         CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         END IF
+        END DO
+*
+*  Compute (U \P' * B) -> B    [ (U \P' * B) ]
+*
+        CALL CTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (U \P' * B) ]
+*       
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL CSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSEIF ( I .GT. 1) THEN
+               IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
+                  AKM1K = WORK(I)
+                  AKM1 = A( I-1, I-1 ) / AKM1K
+                  AK = A( I, I ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 15 J = 1, NRHS
+                     BKM1 = B( I-1, J ) / AKM1K
+                     BK = B( I, J ) / AKM1K
+                     B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 15              CONTINUE
+               I = I - 1
+               ENDIF
+            ENDIF
+            I = I - 1
+         END DO
+*
+*      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
+*
+         CALL CTRSM('L','U','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
+*
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L'.
+*
+*       P' * B  
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K and -IPIV(K+1).
+            KP = -IPIV( K+1 )
+            IF( KP.EQ.-IPIV( K ) )
+     $         CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+*  Compute (L \P' * B) -> B    [ (L \P' * B) ]
+*
+        CALL CTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (L \P' * B) ]
+*       
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL CSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSE
+                  AKM1K = WORK(I)
+                  AKM1 = A( I, I ) / AKM1K
+                  AK = A( I+1, I+1 ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 25 J = 1, NRHS
+                     BKM1 = B( I, J ) / AKM1K
+                     BK = B( I+1, J ) / AKM1K
+                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 25              CONTINUE
+                  I = I + 1
+            ENDIF
+            I = I + 1
+         END DO
+*
+*  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
+* 
+        CALL CTRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
+*
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
+     $         CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         ENDIF
+        END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of CSYTRS2
+*
+      END
diff --git a/SRC/dsyconv.f b/SRC/dsyconv.f
new file mode 100644
index 00000000..0666098a
--- /dev/null
+++ b/SRC/dsyconv.f
@@ -0,0 +1,302 @@
+      SUBROUTINE DSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO, WAY
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      DOUBLE PRECISION   A( LDA, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DSYCONV convert A given by TRF into L and D and vice-versa.
+*  Get Non-diag elements of D (returned in workspace) and 
+*  apply or reverse permutation done in TRF.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+* 
+*  WAY     (input) CHARACTER*1
+*          = 'C': Convert 
+*          = 'R': Revert
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by DSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by DSYTRF.
+*
+* WORK     (workspace) DOUBLE PRECISION array, dimension (N)
+*
+* LWORK    (input) INTEGER
+*          The length of WORK.  LWORK >=1. 
+*          LWORK = N
+*
+*          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.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ZERO
+      PARAMETER          ( ZERO = 0.0D+0 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Local Scalars ..
+      LOGICAL            UPPER, CONVERT
+      INTEGER            I, IP
+      DOUBLE PRECISION   TEMP
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      CONVERT = LSAME( WAY, 'C' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN
+         INFO = -2
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DSYCONV', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*      A is UPPER
+*
+*      Convert A (A is upper)
+*
+*        Convert VALUE
+*
+         IF ( CONVERT ) THEN
+            I=N
+            WORK(1)=ZERO
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I-1,I)
+                  A(I-1,I)=ZERO
+                  I=I-1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Convert PERMUTATIONS
+*  
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0) THEN
+               IP=IPIV(I)
+               IF( I .LT. N) THEN
+                  DO 12 J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+ 12            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+               IF( I .LT. N) THEN
+             DO 13 J= I+1,N
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I-1,J)
+                 A(I-1,J)=TEMP
+ 13            CONTINUE
+                ENDIF
+                I=I-1
+           ENDIF
+           I=I-1
+        END DO
+
+         ELSE
+*
+*      Revert A (A is upper)
+*
+*
+*        Revert PERMUTATIONS
+*  
+            I=1
+            DO WHILE ( I .LE. N )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF( I .LT. N) THEN
+                  DO J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+                  END DO
+                  ENDIF
+               ELSE
+                 IP=-IPIV(I)
+                 I=I+1
+                 IF( I .LT. N) THEN
+                    DO J= I+1,N
+                       TEMP=A(IP,J)
+                       A(IP,J)=A(I-1,J)
+                       A(I-1,J)=TEMP
+                    END DO
+                 ENDIF
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Revert VALUE
+*
+            I=N
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  A(I-1,I)=WORK(I)
+                  I=I-1
+               ENDIF
+               I=I-1
+            END DO
+         END IF
+      ELSE
+*
+*      A is LOWER
+*
+         IF ( CONVERT ) THEN
+*
+*      Convert A (A is lower)
+*
+*
+*        Convert VALUE
+*
+            I=1
+            WORK(N)=ZERO
+            DO WHILE ( I .LE. N )
+               IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I+1,I)
+                  A(I+1,I)=ZERO
+                  I=I+1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Convert PERMUTATIONS
+*
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+               IP=IPIV(I)
+               IF (I .GT. 1) THEN
+               DO 22 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I,J)
+                 A(I,J)=TEMP
+ 22            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+              IF (I .GT. 1) THEN
+              DO 23 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I+1,J)
+                 A(I+1,J)=TEMP
+ 23           CONTINUE
+              ENDIF
+              I=I+1
+           ENDIF
+           I=I+1
+        END DO
+         ELSE
+*
+*      Revert A (A is lower)
+*
+*
+*        Revert PERMUTATIONS
+*
+            I=N
+            DO WHILE ( I .GE. 1 )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I,J)
+                        A(I,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ELSE
+                  IP=-IPIV(I)
+                  I=I-1
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I+1,J)
+                        A(I+1,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Revert VALUE
+*
+            I=1
+            DO WHILE ( I .LE. N-1 )
+               IF( IPIV(I) .LT. ZERO ) THEN
+                  A(I+1,I)=WORK(I)
+                  I=I+1
+               ENDIF
+               I=I+1
+            END DO
+         END IF
+      END IF
+
+      RETURN
+*
+*     End of DSYCONV
+*
+      END
diff --git a/SRC/dsysv.f b/SRC/dsysv.f
index 3abe723e..f8dbac5b 100644
--- a/SRC/dsysv.f
+++ b/SRC/dsysv.f
@@ -4,7 +4,7 @@
 *  -- LAPACK driver routine (version 3.2) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2006
+*     May 2010
 *
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
@@ -105,7 +105,7 @@
 *
 *     .. Local Scalars ..
       LOGICAL            LQUERY
-      INTEGER            LWKOPT, NB
+      INTEGER            LWKOPT, NB, IINFO
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -113,7 +113,7 @@
       EXTERNAL           LSAME, ILAENV
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           DSYTRF, DSYTRS, XERBLA
+      EXTERNAL           DSYCONV, DSYTRF, DSYTRS2, XERBLA
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MAX
@@ -123,6 +123,7 @@
 *     Test the input parameters.
 *
       INFO = 0
+      IINFO = 0
       LQUERY = ( LWORK.EQ.-1 )
       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -160,9 +161,17 @@
       CALL DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN
 *
+*        Convert A
+*
+         CALL DSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
 *        Solve the system A*X = B, overwriting B with X.
 *
-         CALL DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         CALL DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO )
+*
+*        Revert A
+*
+         CALL DSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
 *
       END IF
 *
diff --git a/SRC/dsytrs2.f b/SRC/dsytrs2.f
new file mode 100644
index 00000000..658aed24
--- /dev/null
+++ b/SRC/dsytrs2.f
@@ -0,0 +1,272 @@
+      SUBROUTINE DSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
+     $                    WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  DSYTRS2 solves a system of linear equations A*X = B with a real
+*  symmetric matrix A using the factorization A = U*D*U**T or
+*  A = L*D*L**T computed by DSYTRF and converted by DSYCONV.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  NRHS    (input) INTEGER
+*          The number of right hand sides, i.e., the number of columns
+*          of the matrix B.  NRHS >= 0.
+*
+*  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by DSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by DSYTRF.
+*
+*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
+*          On entry, the right hand side matrix B.
+*          On exit, the solution matrix X.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,N).
+*
+*  WORK    (workspace) REAL array, dimension (N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, J, K, KP
+      DOUBLE PRECISION   AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DSCAL, DSWAP, DTRSM, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+      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
+      ELSE IF( NRHS.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'DSYTRS2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. NRHS.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Solve A*X = B, where A = U*D*U'.
+*
+*       P' * B  
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( KP.EQ.-IPIV( K-1 ) )
+     $         CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         END IF
+        END DO
+*
+*  Compute (U \P' * B) -> B    [ (U \P' * B) ]
+*
+        CALL DTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (U \P' * B) ]
+*       
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSEIF ( I .GT. 1) THEN
+               IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
+                  AKM1K = WORK(I)
+                  AKM1 = A( I-1, I-1 ) / AKM1K
+                  AK = A( I, I ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 15 J = 1, NRHS
+                     BKM1 = B( I-1, J ) / AKM1K
+                     BK = B( I, J ) / AKM1K
+                     B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 15              CONTINUE
+               I = I - 1
+               ENDIF
+            ENDIF
+            I = I - 1
+         END DO
+*
+*      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
+*
+         CALL DTRSM('L','U','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
+*
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
+     $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L'.
+*
+*       P' * B  
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K and -IPIV(K+1).
+            KP = -IPIV( K+1 )
+            IF( KP.EQ.-IPIV( K ) )
+     $         CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+*  Compute (L \P' * B) -> B    [ (L \P' * B) ]
+*
+        CALL DTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (L \P' * B) ]
+*       
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL DSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSE
+                  AKM1K = WORK(I)
+                  AKM1 = A( I, I ) / AKM1K
+                  AK = A( I+1, I+1 ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 25 J = 1, NRHS
+                     BKM1 = B( I, J ) / AKM1K
+                     BK = B( I+1, J ) / AKM1K
+                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 25              CONTINUE
+                  I = I + 1
+            ENDIF
+            I = I + 1
+         END DO
+*
+*  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
+* 
+        CALL DTRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
+*
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
+     $         CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         ENDIF
+        END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of DSYTRS2
+*
+      END
diff --git a/SRC/ssyconv.f b/SRC/ssyconv.f
new file mode 100644
index 00000000..5a966225
--- /dev/null
+++ b/SRC/ssyconv.f
@@ -0,0 +1,302 @@
+      SUBROUTINE SSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO, WAY
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      REAL               A( LDA, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SSYCONV convert A given by TRF into L and D and vice-versa.
+*  Get Non-diag elements of D (returned in workspace) and 
+*  apply or reverse permutation done in TRF.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+* 
+*  WAY     (input) CHARACTER*1
+*          = 'C': Convert 
+*          = 'R': Revert
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input) REAL array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by SSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by SSYTRF.
+*
+* WORK     (workspace) REAL array, dimension (N)
+*
+* LWORK    (input) INTEGER
+*          The length of WORK.  LWORK >=1. 
+*          LWORK = N
+*
+*          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.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO
+      PARAMETER          ( ZERO = 0.0E+0 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Local Scalars ..
+      LOGICAL            UPPER, CONVERT
+      INTEGER            I, IP
+      REAL               TEMP
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      CONVERT = LSAME( WAY, 'C' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN
+         INFO = -2
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SSYCONV', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*      A is UPPER
+*
+*      Convert A (A is upper)
+*
+*        Convert VALUE
+*
+         IF ( CONVERT ) THEN
+            I=N
+            WORK(1)=ZERO
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I-1,I)
+                  A(I-1,I)=ZERO
+                  I=I-1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Convert PERMUTATIONS
+*  
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0) THEN
+               IP=IPIV(I)
+               IF( I .LT. N) THEN
+                  DO 12 J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+ 12            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+               IF( I .LT. N) THEN
+             DO 13 J= I+1,N
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I-1,J)
+                 A(I-1,J)=TEMP
+ 13            CONTINUE
+                ENDIF
+                I=I-1
+           ENDIF
+           I=I-1
+        END DO
+
+         ELSE
+*
+*      Revert A (A is upper)
+*
+*
+*        Revert PERMUTATIONS
+*  
+            I=1
+            DO WHILE ( I .LE. N )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF( I .LT. N) THEN
+                  DO J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+                  END DO
+                  ENDIF
+               ELSE
+                 IP=-IPIV(I)
+                 I=I+1
+                 IF( I .LT. N) THEN
+                    DO J= I+1,N
+                       TEMP=A(IP,J)
+                       A(IP,J)=A(I-1,J)
+                       A(I-1,J)=TEMP
+                    END DO
+                 ENDIF
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Revert VALUE
+*
+            I=N
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  A(I-1,I)=WORK(I)
+                  I=I-1
+               ENDIF
+               I=I-1
+            END DO
+         END IF
+      ELSE
+*
+*      A is LOWER
+*
+         IF ( CONVERT ) THEN
+*
+*      Convert A (A is lower)
+*
+*
+*        Convert VALUE
+*
+            I=1
+            WORK(N)=ZERO
+            DO WHILE ( I .LE. N )
+               IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I+1,I)
+                  A(I+1,I)=ZERO
+                  I=I+1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Convert PERMUTATIONS
+*
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+               IP=IPIV(I)
+               IF (I .GT. 1) THEN
+               DO 22 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I,J)
+                 A(I,J)=TEMP
+ 22            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+              IF (I .GT. 1) THEN
+              DO 23 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I+1,J)
+                 A(I+1,J)=TEMP
+ 23           CONTINUE
+              ENDIF
+              I=I+1
+           ENDIF
+           I=I+1
+        END DO
+         ELSE
+*
+*      Revert A (A is lower)
+*
+*
+*        Revert PERMUTATIONS
+*
+            I=N
+            DO WHILE ( I .GE. 1 )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I,J)
+                        A(I,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ELSE
+                  IP=-IPIV(I)
+                  I=I-1
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I+1,J)
+                        A(I+1,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Revert VALUE
+*
+            I=1
+            DO WHILE ( I .LE. N-1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  A(I+1,I)=WORK(I)
+                  I=I+1
+               ENDIF
+               I=I+1
+            END DO
+         END IF
+      END IF
+
+      RETURN
+*
+*     End of SSYCONV
+*
+      END
diff --git a/SRC/ssysv.f b/SRC/ssysv.f
index bcefac89..a7262b65 100644
--- a/SRC/ssysv.f
+++ b/SRC/ssysv.f
@@ -4,7 +4,7 @@
 *  -- LAPACK driver routine (version 3.2) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
-*     November 2006
+*     May 2010
 *
 *     .. Scalar Arguments ..
       CHARACTER          UPLO
@@ -105,7 +105,7 @@
 *
 *     .. Local Scalars ..
       LOGICAL            LQUERY
-      INTEGER            LWKOPT, NB
+      INTEGER            LWKOPT, NB, IINFO
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -113,7 +113,7 @@
       EXTERNAL           ILAENV, LSAME
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           SSYTRF, SSYTRS, XERBLA
+      EXTERNAL           SSYCONV, SSYTRF, SSYTRS2, XERBLA
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MAX
@@ -123,6 +123,7 @@
 *     Test the input parameters.
 *
       INFO = 0
+      IINFO = 0
       LQUERY = ( LWORK.EQ.-1 )
       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -160,9 +161,17 @@
       CALL SSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN
 *
+*        Convert A
+*
+         CALL SSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
 *        Solve the system A*X = B, overwriting B with X.
 *
-         CALL SSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         CALL SSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO )
+*
+*        Revert A
+*
+         CALL SSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
 *
       END IF
 *
diff --git a/SRC/ssytrs2.f b/SRC/ssytrs2.f
new file mode 100644
index 00000000..cb98bbd9
--- /dev/null
+++ b/SRC/ssytrs2.f
@@ -0,0 +1,272 @@
+      SUBROUTINE SSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
+     $                    WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      REAL               A( LDA, * ), B( LDB, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  SSYTRS2 solves a system of linear equations A*X = B with a real
+*  symmetric matrix A using the factorization A = U*D*U**T or
+*  A = L*D*L**T computed by SSYTRF and converted by SSYCONV.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  NRHS    (input) INTEGER
+*          The number of right hand sides, i.e., the number of columns
+*          of the matrix B.  NRHS >= 0.
+*
+*  A       (input) REAL array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by SSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by SSYTRF.
+*
+*  B       (input/output) REAL array, dimension (LDB,NRHS)
+*          On entry, the right hand side matrix B.
+*          On exit, the solution matrix X.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,N).
+*
+*  WORK    (workspace) REAL array, dimension (N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE
+      PARAMETER          ( ONE = 1.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, J, K, KP
+      REAL               AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SSCAL, SSWAP, STRSM, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+      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
+      ELSE IF( NRHS.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'SSYTRS2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. NRHS.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Solve A*X = B, where A = U*D*U'.
+*
+*       P' * B  
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( KP.EQ.-IPIV( K-1 ) )
+     $         CALL SSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         END IF
+        END DO
+*
+*  Compute (U \P' * B) -> B    [ (U \P' * B) ]
+*
+        CALL STRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (U \P' * B) ]
+*       
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL SSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSEIF ( I .GT. 1) THEN
+               IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
+                  AKM1K = WORK(I)
+                  AKM1 = A( I-1, I-1 ) / AKM1K
+                  AK = A( I, I ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 15 J = 1, NRHS
+                     BKM1 = B( I-1, J ) / AKM1K
+                     BK = B( I, J ) / AKM1K
+                     B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 15              CONTINUE
+               I = I - 1
+               ENDIF
+            ENDIF
+            I = I - 1
+         END DO
+*
+*      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
+*
+         CALL STRSM('L','U','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
+*
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
+     $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L'.
+*
+*       P' * B  
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K and -IPIV(K+1).
+            KP = -IPIV( K+1 )
+            IF( KP.EQ.-IPIV( K ) )
+     $         CALL SSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+*  Compute (L \P' * B) -> B    [ (L \P' * B) ]
+*
+        CALL STRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (L \P' * B) ]
+*       
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL SSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSE
+                  AKM1K = WORK(I)
+                  AKM1 = A( I, I ) / AKM1K
+                  AK = A( I+1, I+1 ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 25 J = 1, NRHS
+                     BKM1 = B( I, J ) / AKM1K
+                     BK = B( I+1, J ) / AKM1K
+                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 25              CONTINUE
+                  I = I + 1
+            ENDIF
+            I = I + 1
+         END DO
+*
+*  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
+* 
+        CALL STRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
+*
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
+     $         CALL SSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         ENDIF
+        END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of SSYTRS2
+*
+      END
diff --git a/SRC/zsyconv.f b/SRC/zsyconv.f
new file mode 100644
index 00000000..205c650e
--- /dev/null
+++ b/SRC/zsyconv.f
@@ -0,0 +1,302 @@
+      SUBROUTINE ZSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO, WAY
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      DOUBLE COMPLEX     A( LDA, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZSYCONV convert A given by TRF into L and D and vice-versa.
+*  Get Non-diag elements of D (returned in workspace) and 
+*  apply or reverse permutation done in TRF.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+* 
+*  WAY     (input) CHARACTER*1
+*          = 'C': Convert 
+*          = 'R': Revert
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  A       (input) DOUBLE COMPLEX array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by ZSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by ZSYTRF.
+*
+* WORK     (workspace) DOUBLE COMPLEX array, dimension (N)
+*
+* LWORK    (input) INTEGER
+*          The length of WORK.  LWORK >=1. 
+*          LWORK = N
+*
+*          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.
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE COMPLEX     ZERO
+      PARAMETER          ( ZERO = (0.0D+0,0.0D+0) )
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Local Scalars ..
+      LOGICAL            UPPER, CONVERT
+      INTEGER            I, IP
+      DOUBLE COMPLEX     TEMP
+*     ..
+*     .. Executable Statements ..
+*
+      INFO = 0
+      UPPER = LSAME( UPLO, 'U' )
+      CONVERT = LSAME( WAY, 'C' )
+      IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+         INFO = -1
+      ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN
+         INFO = -2
+      ELSE IF( N.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZSYCONV', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*      A is UPPER
+*
+*      Convert A (A is upper)
+*
+*        Convert VALUE
+*
+         IF ( CONVERT ) THEN
+            I=N
+            WORK(1)=ZERO
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I-1,I)
+                  A(I-1,I)=ZERO
+                  I=I-1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Convert PERMUTATIONS
+*  
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0) THEN
+               IP=IPIV(I)
+               IF( I .LT. N) THEN
+                  DO 12 J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+ 12            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+               IF( I .LT. N) THEN
+             DO 13 J= I+1,N
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I-1,J)
+                 A(I-1,J)=TEMP
+ 13            CONTINUE
+                ENDIF
+                I=I-1
+           ENDIF
+           I=I-1
+        END DO
+
+         ELSE
+*
+*      Revert A (A is upper)
+*
+*
+*        Revert PERMUTATIONS
+*  
+            I=1
+            DO WHILE ( I .LE. N )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF( I .LT. N) THEN
+                  DO J= I+1,N
+                    TEMP=A(IP,J)
+                    A(IP,J)=A(I,J)
+                    A(I,J)=TEMP
+                  END DO
+                  ENDIF
+               ELSE
+                 IP=-IPIV(I)
+                 I=I+1
+                 IF( I .LT. N) THEN
+                    DO J= I+1,N
+                       TEMP=A(IP,J)
+                       A(IP,J)=A(I-1,J)
+                       A(I-1,J)=TEMP
+                    END DO
+                 ENDIF
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Revert VALUE
+*
+            I=N
+            DO WHILE ( I .GT. 1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  A(I-1,I)=WORK(I)
+                  I=I-1
+               ENDIF
+               I=I-1
+            END DO
+         END IF
+      ELSE
+*
+*      A is LOWER
+*
+         IF ( CONVERT ) THEN
+*
+*      Convert A (A is lower)
+*
+*
+*        Convert VALUE
+*
+            I=1
+            WORK(N)=ZERO
+            DO WHILE ( I .LE. N )
+               IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN
+                  WORK(I)=A(I+1,I)
+                  A(I+1,I)=ZERO
+                  I=I+1
+               ELSE
+                  WORK(I)=ZERO
+               ENDIF
+               I=I+1
+            END DO
+*
+*        Convert PERMUTATIONS
+*
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+               IP=IPIV(I)
+               IF (I .GT. 1) THEN
+               DO 22 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I,J)
+                 A(I,J)=TEMP
+ 22            CONTINUE
+               ENDIF
+            ELSE
+              IP=-IPIV(I)
+              IF (I .GT. 1) THEN
+              DO 23 J= 1,I-1
+                 TEMP=A(IP,J)
+                 A(IP,J)=A(I+1,J)
+                 A(I+1,J)=TEMP
+ 23           CONTINUE
+              ENDIF
+              I=I+1
+           ENDIF
+           I=I+1
+        END DO
+         ELSE
+*
+*      Revert A (A is lower)
+*
+*
+*        Revert PERMUTATIONS
+*
+            I=N
+            DO WHILE ( I .GE. 1 )
+               IF( IPIV(I) .GT. 0 ) THEN
+                  IP=IPIV(I)
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I,J)
+                        A(I,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ELSE
+                  IP=-IPIV(I)
+                  I=I-1
+                  IF (I .GT. 1) THEN
+                     DO J= 1,I-1
+                        TEMP=A(I+1,J)
+                        A(I+1,J)=A(IP,J)
+                        A(IP,J)=TEMP
+                     END DO
+                  ENDIF
+               ENDIF
+               I=I-1
+            END DO
+*
+*        Revert VALUE
+*
+            I=1
+            DO WHILE ( I .LE. N-1 )
+               IF( IPIV(I) .LT. 0 ) THEN
+                  A(I+1,I)=WORK(I)
+                  I=I+1
+               ENDIF
+               I=I+1
+            END DO
+         END IF
+      END IF
+
+      RETURN
+*
+*     End of ZSYCONV
+*
+      END
diff --git a/SRC/zsysv.f b/SRC/zsysv.f
index 2520e165..bf7baab7 100644
--- a/SRC/zsysv.f
+++ b/SRC/zsysv.f
@@ -105,7 +105,7 @@
 *
 *     .. Local Scalars ..
       LOGICAL            LQUERY
-      INTEGER            LWKOPT, NB
+      INTEGER            LWKOPT, NB, IINFO
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
@@ -113,7 +113,7 @@
       EXTERNAL           LSAME, ILAENV
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL           XERBLA, ZSYTRF, ZSYTRS
+      EXTERNAL           ZSYCONV, XERBLA, ZSYTRF, ZSYTRS2
 *     ..
 *     .. Intrinsic Functions ..
       INTRINSIC          MAX
@@ -123,6 +123,7 @@
 *     Test the input parameters.
 *
       INFO = 0
+      IINFO = 0
       LQUERY = ( LWORK.EQ.-1 )
       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
          INFO = -1
@@ -160,9 +161,17 @@
       CALL ZSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
       IF( INFO.EQ.0 ) THEN
 *
+*        Convert A
+*
+         CALL ZSYCONV( UPLO, 'C', N, A, LDA, IPIV, WORK, IINFO )
+*
 *        Solve the system A*X = B, overwriting B with X.
 *
-         CALL ZSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+         CALL ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, INFO )
+*
+*        Revert A
+*
+         CALL ZSYCONV( UPLO, 'R', N, A, LDA, IPIV, WORK, IINFO )
 *
       END IF
 *
diff --git a/SRC/zsytrs2.f b/SRC/zsytrs2.f
new file mode 100644
index 00000000..88d54408
--- /dev/null
+++ b/SRC/zsytrs2.f
@@ -0,0 +1,272 @@
+      SUBROUTINE ZSYTRS2( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, 
+     $                    WORK, INFO )
+*
+*  -- LAPACK PROTOTYPE routine (version 3.2) --
+*
+*  -- Written by Julie Langou of the Univ. of TN    --
+*     May 2010
+*
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO
+      INTEGER            INFO, LDA, LDB, N, NRHS
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      DOUBLE COMPLEX   A( LDA, * ), B( LDB, * ), WORK( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZSYTRS2 solves a system of linear equations A*X = B with a real
+*  symmetric matrix A using the factorization A = U*D*U**T or
+*  A = L*D*L**T computed by ZSYTRF and converted by ZSYCONV.
+*
+*  Arguments
+*  =========
+*
+*  UPLO    (input) CHARACTER*1
+*          Specifies whether the details of the factorization are stored
+*          as an upper or lower triangular matrix.
+*          = 'U':  Upper triangular, form is A = U*D*U**T;
+*          = 'L':  Lower triangular, form is A = L*D*L**T.
+*
+*  N       (input) INTEGER
+*          The order of the matrix A.  N >= 0.
+*
+*  NRHS    (input) INTEGER
+*          The number of right hand sides, i.e., the number of columns
+*          of the matrix B.  NRHS >= 0.
+*
+*  A       (input) DOUBLE COMPLEX array, dimension (LDA,N)
+*          The block diagonal matrix D and the multipliers used to
+*          obtain the factor U or L as computed by ZSYTRF.
+*
+*  LDA     (input) INTEGER
+*          The leading dimension of the array A.  LDA >= max(1,N).
+*
+*  IPIV    (input) INTEGER array, dimension (N)
+*          Details of the interchanges and the block structure of D
+*          as determined by ZSYTRF.
+*
+*  B       (input/output) DOUBLE COMPLEX array, dimension (LDB,NRHS)
+*          On entry, the right hand side matrix B.
+*          On exit, the solution matrix X.
+*
+*  LDB     (input) INTEGER
+*          The leading dimension of the array B.  LDB >= max(1,N).
+*
+*  WORK    (workspace) REAL array, dimension (N)
+*
+*  INFO    (output) INTEGER
+*          = 0:  successful exit
+*          < 0:  if INFO = -i, the i-th argument had an illegal value
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE COMPLEX     ONE
+      PARAMETER          ( ONE = (1.0D+0,0.0D+0) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            UPPER
+      INTEGER            I, J, K, KP
+      DOUBLE COMPLEX     AK, AKM1, AKM1K, BK, BKM1, DENOM
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZSCAL, ZSWAP, ZTRSM, XERBLA
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          MAX
+*     ..
+*     .. Executable Statements ..
+*
+      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
+      ELSE IF( NRHS.LT.0 ) THEN
+         INFO = -3
+      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+         INFO = -5
+      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+         INFO = -8
+      END IF
+      IF( INFO.NE.0 ) THEN
+         CALL XERBLA( 'ZSYTRS2', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 .OR. NRHS.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Solve A*X = B, where A = U*D*U'.
+*
+*       P' * B  
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( KP.EQ.-IPIV( K-1 ) )
+     $         CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         END IF
+        END DO
+*
+*  Compute (U \P' * B) -> B    [ (U \P' * B) ]
+*
+        CALL ZTRSM('L','U','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (U \P' * B) ]
+*       
+         I=N
+         DO WHILE ( I .GE. 1 )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSEIF ( I .GT. 1) THEN
+               IF ( IPIV(I-1) .EQ. IPIV(I) ) THEN
+                  AKM1K = WORK(I)
+                  AKM1 = A( I-1, I-1 ) / AKM1K
+                  AK = A( I, I ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 15 J = 1, NRHS
+                     BKM1 = B( I-1, J ) / AKM1K
+                     BK = B( I, J ) / AKM1K
+                     B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 15              CONTINUE
+               I = I - 1
+               ENDIF
+            ENDIF
+            I = I - 1
+         END DO
+*
+*      Compute (U' \ B) -> B   [ U' \ (D \ (U \P' * B) ) ]
+*
+         CALL ZTRSM('L','U','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (U' \ (D \ (U \P' * B) )) ]
+*
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K .LT. N .AND. KP.EQ.-IPIV( K+1 ) )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+      ELSE
+*
+*        Solve A*X = B, where A = L*D*L'.
+*
+*       P' * B  
+        K=1
+        DO WHILE ( K .LE. N )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K and -IPIV(K+1).
+            KP = -IPIV( K+1 )
+            IF( KP.EQ.-IPIV( K ) )
+     $         CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K+2
+         ENDIF
+        END DO
+*
+*  Compute (L \P' * B) -> B    [ (L \P' * B) ]
+*
+        CALL ZTRSM('L','L','N','U',N,NRHS,ONE,A,N,B,N)
+*
+*  Compute D \ B -> B   [ D \ (L \P' * B) ]
+*       
+         I=1
+         DO WHILE ( I .LE. N )
+            IF( IPIV(I) .GT. 0 ) THEN
+              CALL ZSCAL( NRHS, ONE / A( I, I ), B( I, 1 ), N )
+            ELSE
+                  AKM1K = WORK(I)
+                  AKM1 = A( I, I ) / AKM1K
+                  AK = A( I+1, I+1 ) / AKM1K
+                  DENOM = AKM1*AK - ONE
+                  DO 25 J = 1, NRHS
+                     BKM1 = B( I, J ) / AKM1K
+                     BK = B( I+1, J ) / AKM1K
+                     B( I, J ) = ( AK*BKM1-BK ) / DENOM
+                     B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ 25              CONTINUE
+                  I = I + 1
+            ENDIF
+            I = I + 1
+         END DO
+*
+*  Compute (L' \ B) -> B   [ L' \ (D \ (L \P' * B) ) ]
+* 
+        CALL ZTRSM('L','L','T','U',N,NRHS,ONE,A,N,B,N)
+*
+*       P * B  [ P * (L' \ (D \ (L \P' * B) )) ]
+*
+        K=N
+        DO WHILE ( K .GE. 1 )
+         IF( IPIV( K ).GT.0 ) THEN
+*           1 x 1 diagonal block
+*           Interchange rows K and IPIV(K).
+            KP = IPIV( K )
+            IF( KP.NE.K )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-1
+         ELSE
+*           2 x 2 diagonal block
+*           Interchange rows K-1 and -IPIV(K).
+            KP = -IPIV( K )
+            IF( K.GT.1 .AND. KP.EQ.-IPIV( K-1 ) )
+     $         CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+            K=K-2
+         ENDIF
+        END DO
+*
+      END IF
+*
+      RETURN
+*
+*     End of ZSYTRS2
+*
+      END