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.

1 files changed, 302 insertions, 0 deletions

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