Added (S,D,C,Z) (SY,HE) routines, drivers for new rook code

Close #82

Added routines for new factorization code for symmetric indefinite
( or Hermitian indefinite ) matrices with bounded Bunch-Kaufman
( rook ) pivoting algorithm.

New more efficient storage format for factors U ( or L ),
block-diagonal matrix D, and pivoting information stored in IPIV:

factor L is stored explicitly in lower triangle of A;
diagonal of D is stored on the diagonal of A;
subdiagonal elements of D are stored in array E;
IPIV format is the same as in *_ROOK routines, but differs
from SY Bunch-Kaufman routines (e.g. *SYTRF).
The factorization output of these new rook _RK routines is not
compatible
with the existing _ROOK routines and vice versa. This new factorization
format is designed in such a way, that there is a possibility in the
future
to write new Bunch-Kaufman routines that conform to this new
factorization format.
Then the future Bunch-Kaufman routines could share solver
*TRS_3,inversion *TRI_3
and condition estimator *CON_3.

To convert between the factorization formats in both ways the following
routines
are developed:

CONVERSION ROUTINES BETWEEN FACTORIZATION FORMATS

DOUBLE PRECISION (symmetric indefinite matrices):

new file:   SRC/dsyconvf.f
new file:   SRC/dsyconvf_rook.f
REAL (symmetric indefinite matrices):

new file:   SRC/csyconvf.f
new file:   SRC/csyconvf_rook.f
COMPLEX*16 (symmetric indefinite and Hermitian indefinite matrices):

new file:   SRC/zsyconvf.f
new file:   SRC/zsyconvf_rook.f
COMPLEX (symmetric indefinite and Hermitian indefinite matrices):

new file:   SRC/ssyconvf.f
new file:   SRC/ssyconvf_rook.f
*SYCONVF routine converts between old Bunch-Kaufman storage format (
denote (L1,D1,IPIV1) )
that is used by *SYTRF and new rook storage format ( denote (L2,D2,
IPIV2))
that is used by *SYTRF_RK

*SYCONVF_ROOK routine between old rook storage format ( denote
(L1,D1,IPIV2) )
that is used by *SYTRF_ROOK and new rook storage format ( denote
(L2,D2, IPIV2))
that is used by *SYTRF_RK

ROUTINES AND DRIVERS

DOUBLE PRECISION (symmetric indefinite matrices):

new file:   SRC/dsytf2_rk.f BLAS2 unblocked factorization
new file:   SRC/dlasyf_rk.f BLAS3 auxiliary blocked partial
factorization
new file:   SRC/dsytrf_rk.f BLAS3 blocked factorization
new file:   SRC/dsytrs_3.f BLAS3 solver
new file:   SRC/dsycon_3.f BLAS3 condition number estimator
new file:   SRC/dsytri_3.f BLAS3 inversion, sets the size of work array
and calls *sytri_3x
new file:   SRC/dsytri_3x.f BLAS3 auxiliary inversion, actually
computes blocked inversion
new file:   SRC/dsysv_rk.f BLAS3 solver driver
REAL (symmetric indefinite matrices):

new file:   SRC/ssytf2_rk.f BLAS2 unblocked factorization
new file:   SRC/slasyf_rk.f BLAS3 auxiliary blocked partial
factorization
new file:   SRC/ssytrf_rk.f BLAS3 blocked factorization
new file:   SRC/ssytrs_3.f BLAS3 solver
new file:   SRC/ssycon_3.f BLAS3 condition number estimator
new file:   SRC/ssytri_3.f BLAS3 inversion, sets the size of work array
and calls *sytri_3x
new file:   SRC/ssytri_3x.f BLAS3 auxiliary inversion, actually
computes blocked inversion
new file:   SRC/ssysv_rk.f BLAS3 solver driver
COMPLEX*16 (symmetric indefinite matrices):

new file:   SRC/zsytf2_rk.f BLAS2 unblocked factorization
new file:   SRC/zlasyf_rk.f BLAS3 auxiliary blocked partial
factorization
new file:   SRC/zsytrf_rk.f BLAS3 blocked factorization
new file:   SRC/zsytrs_3.f BLAS3 solver
new file:   SRC/zsycon_3.f BLAS3 condition number estimator
new file:   SRC/zsytri_3.f BLAS3 inversion, sets the size of work array
and calls *sytri_3x
new file:   SRC/zsytri_3x.f BLAS3 auxiliary inversion, actually
computes blocked inversion
new file:   SRC/zsysv_rk.f BLAS3 solver driver
COMPLEX*16 (Hermitian indefinite matrices):

new file:   SRC/zhetf2_rk.f BLAS2 unblocked factorization
new file:   SRC/zlahef_rk.f BLAS3 auxiliary blocked partial
factorization
new file:   SRC/zhetrf_rk.f BLAS3 blocked factorization
new file:   SRC/zhetrs_3.f BLAS3 solver
new file:   SRC/zhecon_3.f BLAS3 condition number estimator
new file:   SRC/zhetri_3.f BLAS3 inversion, sets the size of work array
and calls *sytri_3x
new file:   SRC/zhetri_3x.f BLAS3 auxiliary inversion, actually
computes blocked inversion
new file:   SRC/zhesv_rk.f BLAS3 solver driver
COMPLEX (symmetric indefinite matrices):

new file:   SRC/csytf2_rk.f BLAS2 unblocked factorization
new file:   SRC/clasyf_rk.f BLAS3 auxiliary blocked partial
factorization
new file:   SRC/csytrf_rk.f BLAS3 blocked factorization
new file:   SRC/csytrs_3.f BLAS3 solver
new file:   SRC/csycon_3.f BLAS3 condition number estimator
new file:   SRC/csytri_3.f BLAS3 inversion, sets the size of work array
and calls *sytri_3x
new file:   SRC/csytri_3x.f BLAS3 auxiliary inversion, actually
computes blocked inversion
new file:   SRC/csysv_rk.f BLAS3 solver driver
COMPLEX (Hermitian indefinite matrices):

new file:   SRC/chetf2_rk.f BLAS2 unblocked factorization
new file:   SRC/clahef_rk.f BLAS3 auxiliary blocked partial
factorization
new file:   SRC/chetrf_rk.f BLAS3 blocked factorization
new file:   SRC/chetrs_3.f BLAS3 solver
new file:   SRC/checon_3.f BLAS3 condition number estimator
new file:   SRC/chetri_3.f BLAS3 inversion, sets the size of work array
and calls *sytri_3x
new file:   SRC/chetri_3x.f BLAS3 auxiliary inversion, actually
computes blocked inversion
new file:   SRC/chesv_rk.f BLAS3 solver driver
MISC

modified:   SRC/CMakeLists.txt
modified:   SRC/Makefile
TEST CODE

modified:   TESTING/LIN/CMakeLists.txt
modified:   TESTING/LIN/Makefile

modified:   TESTING/LIN/aladhd.f
modified:   TESTING/LIN/alaerh.f
modified:   TESTING/LIN/alahd.f
DOUBLE PRECISION (symmetric indefinite matrices):

modified:   TESTING/LIN/dchkaa.f
modified:   TESTING/LIN/derrsy.f
modified:   TESTING/LIN/derrsyx.f
modified:   TESTING/LIN/derrvx.f
modified:   TESTING/LIN/derrvxx.f

modified:   TESTING/dtest.in

new file:   TESTING/LIN/dchksy_rk.f
new file:   TESTING/LIN/ddrvsy_rk.f
new file:   TESTING/LIN/dsyt01_3.f
REAL (symmetric indefinite matrices):

modified:   TESTING/LIN/schkaa.f
modified:   TESTING/LIN/serrsy.f
modified:   TESTING/LIN/serrsyx.f
modified:   TESTING/LIN/serrvx.f
modified:   TESTING/LIN/serrvxx.f

modified:   TESTING/stest.in

new file:   TESTING/LIN/schksy_rk.f
new file:   TESTING/LIN/sdrvsy_rk.f
new file:   TESTING/LIN/ssyt01_3.f
COMPLEX*16 (symmetric indefinite and Hermitian indefinite matrices):

modified:   TESTING/LIN/zchkaa.f
modified:   TESTING/LIN/zerrsy.f
modified:   TESTING/LIN/zerrsyx.f
modified:   TESTING/LIN/zerrhe.f
modified:   TESTING/LIN/zerrhex.f
modified:   TESTING/LIN/zerrvx.f
modified:   TESTING/LIN/zerrvxx.f

modified:   TESTING/ztest.in

new file:   TESTING/LIN/zchksy_rk.f
new file:   TESTING/LIN/zdrvsy_rk.f
new file:   TESTING/LIN/zsyt01_3.f
new file:   TESTING/LIN/zchkhe_rk.f
new file:   TESTING/LIN/zdrvhe_rk.f
new file:   TESTING/LIN/zhet01_3.f
COMPLEX (symmetric indefinite and Hermitian indefinite matrices):

modified:   TESTING/LIN/cchkaa.f
modified:   TESTING/LIN/cerrsy.f
modified:   TESTING/LIN/cerrsyx.f
modified:   TESTING/LIN/cerrhe.f
modified:   TESTING/LIN/cerrhex.f
modified:   TESTING/LIN/cerrvx.f
modified:   TESTING/LIN/cerrvxx.f

modified:   TESTING/ctest.in

new file:   TESTING/LIN/cchksy_rk.f
new file:   TESTING/LIN/cdrvsy_rk.f
new file:   TESTING/LIN/csyt01_3.f
new file:   TESTING/LIN/cchkhe_rk.f
new file:   TESTING/LIN/cdrvhe_rk.f
new file:   TESTING/LIN/chet01_3.f

1 files changed, 544 insertions, 0 deletions

diff --git a/SRC/ssyconvf_rook.f b/SRC/ssyconvf_rook.f
new file mode 100644
index 00000000..69f04f6d
--- /dev/null
+++ b/SRC/ssyconvf_rook.f
@@ -0,0 +1,544 @@
+*> \brief \b SSYCONVF_ROOK
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+*            http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download SSYCONVF_ROOK + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssyconvf_rook.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssyconvf_rook.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssyconvf_rook.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SSYCONVF_ROOK( UPLO, WAY, N, A, LDA, IPIV, E, INFO )
+*
+*       .. Scalar Arguments ..
+*       CHARACTER          UPLO, WAY
+*       INTEGER            INFO, LDA, N
+*       ..
+*       .. Array Arguments ..
+*       INTEGER            IPIV( * )
+*       REAL               A( LDA, * ), E( * )
+*       ..
+*
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*> If parameter WAY = 'C':
+*> SSYCONVF_ROOK converts the factorization output format used in
+*> SSYTRF_ROOK provided on entry in parameter A into the factorization
+*> output format used in SSYTRF_RK (or SSYTRF_BK) that is stored
+*> on exit in parameters A and E. IPIV format for SSYTRF_ROOK and
+*> SSYTRF_RK (or SSYTRF_BK) is the same and is not converted.
+*>
+*> If parameter WAY = 'R':
+*> SSYCONVF_ROOK performs the conversion in reverse direction, i.e.
+*> converts the factorization output format used in SSYTRF_RK
+*> (or SSYTRF_BK) provided on entry in parametes A and E into
+*> the factorization output format used in SSYTRF_ROOK that is stored
+*> on exit in parameter A. IPIV format for SSYTRF_ROOK and
+*> SSYTRF_RK (or SSYTRF_BK) is the same and is not converted.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*>          UPLO is CHARACTER*1
+*>          Specifies whether the details of the factorization are
+*>          stored as an upper or lower triangular matrix A.
+*>          = 'U':  Upper triangular
+*>          = 'L':  Lower triangular
+*> \endverbatim
+*>
+*> \param[in] WAY
+*> \verbatim
+*>          WAY is CHARACTER*1
+*>          = 'C': Convert
+*>          = 'R': Revert
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the matrix A.  N >= 0.
+*> \endverbatim
+*>
+*> \param[in,out] A
+*> \verbatim
+*>          A is REAL array, dimension (LDA,N)
+*>
+*>          1) If WAY ='C':
+*>
+*>          On entry, contains factorization details in format used in
+*>          SSYTRF_ROOK:
+*>            a) all elements of the symmetric block diagonal
+*>               matrix D on the diagonal of A and on superdiagonal
+*>               (or subdiagonal) of A, and
+*>            b) If UPLO = 'U': multipliers used to obtain factor U
+*>               in the superdiagonal part of A.
+*>               If UPLO = 'L': multipliers used to obtain factor L
+*>               in the superdiagonal part of A.
+*>
+*>          On exit, contains factorization details in format used in
+*>          SSYTRF_RK or SSYTRF_BK:
+*>            a) ONLY diagonal elements of the symmetric block diagonal
+*>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
+*>               (superdiagonal (or subdiagonal) elements of D
+*>                are stored on exit in array E), and
+*>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
+*>               If UPLO = 'L': factor L in the subdiagonal part of A.
+*>
+*>          2) If WAY = 'R':
+*>
+*>          On entry, contains factorization details in format used in
+*>          SSYTRF_RK or SSYTRF_BK:
+*>            a) ONLY diagonal elements of the symmetric block diagonal
+*>               matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
+*>               (superdiagonal (or subdiagonal) elements of D
+*>                are stored on exit in array E), and
+*>            b) If UPLO = 'U': factor U in the superdiagonal part of A.
+*>               If UPLO = 'L': factor L in the subdiagonal part of A.
+*>
+*>          On exit, contains factorization details in format used in
+*>          SSYTRF_ROOK:
+*>            a) all elements of the symmetric block diagonal
+*>               matrix D on the diagonal of A and on superdiagonal
+*>               (or subdiagonal) of A, and
+*>            b) If UPLO = 'U': multipliers used to obtain factor U
+*>               in the superdiagonal part of A.
+*>               If UPLO = 'L': multipliers used to obtain factor L
+*>               in the superdiagonal part of A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*>          LDA is INTEGER
+*>          The leading dimension of the array A.  LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in,out] E
+*> \verbatim
+*>          E is REAL array, dimension (N)
+*>
+*>          1) If WAY ='C':
+*>
+*>          On entry, just a workspace.
+*>
+*>          On exit, contains the superdiagonal (or subdiagonal)
+*>          elements of the symmetric block diagonal matrix D
+*>          with 1-by-1 or 2-by-2 diagonal blocks, where
+*>          If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0;
+*>          If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0.
+*>
+*>          2) If WAY = 'R':
+*>
+*>          On entry, contains the superdiagonal (or subdiagonal)
+*>          elements of the symmetric block diagonal matrix D
+*>          with 1-by-1 or 2-by-2 diagonal blocks, where
+*>          If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not referenced;
+*>          If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.
+*>
+*>          On exit, is not changed
+*> \endverbatim
+*.
+*> \param[in] IPIV
+*> \verbatim
+*>          IPIV is INTEGER array, dimension (N)
+*>          On entry, details of the interchanges and the block
+*>          structure of D as determined:
+*>          1) by SSYTRF_ROOK, if WAY ='C';
+*>          2) by SSYTRF_RK (or SSYTRF_BK), if WAY ='R'.
+*>          The IPIV format is the same for all these routines.
+*>
+*>          On exit, is not changed.
+*> \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 November 2016
+*
+*> \ingroup singleSYcomputational
+*
+*> \par Contributors:
+*  ==================
+*>
+*> \verbatim
+*>
+*>  November 2016,  Igor Kozachenko,
+*>                  Computer Science Division,
+*>                  University of California, Berkeley
+*>
+*> \endverbatim
+*  =====================================================================
+      SUBROUTINE SSYCONVF_ROOK( UPLO, WAY, N, A, LDA, E, IPIV, 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..--
+*     November 2016
+*
+*     .. Scalar Arguments ..
+      CHARACTER          UPLO, WAY
+      INTEGER            INFO, LDA, N
+*     ..
+*     .. Array Arguments ..
+      INTEGER            IPIV( * )
+      REAL               A( LDA, * ), E( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ZERO
+      PARAMETER          ( ZERO = 0.0E+0 )
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*
+*     .. External Subroutines ..
+      EXTERNAL           SSWAP, XERBLA
+*     .. Local Scalars ..
+      LOGICAL            UPPER, CONVERT
+      INTEGER            I, IP, IP2
+*     ..
+*     .. 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( 'SSYCONVF_ROOK', -INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( UPPER ) THEN
+*
+*        Begin A is UPPER
+*
+         IF ( CONVERT ) THEN
+*
+*           Convert A (A is upper)
+*
+*
+*           Convert VALUE
+*
+*           Assign superdiagonal entries of D to array E and zero out
+*           corresponding entries in input storage A
+*
+            I = N
+            E( 1 ) = ZERO
+            DO WHILE ( I.GT.1 )
+               IF( IPIV( I ).LT.0 ) THEN
+                  E( I ) = A( I-1, I )
+                  E( I-1 ) = ZERO
+                  A( I-1, I ) = ZERO
+                  I = I - 1
+               ELSE
+                  E( I ) = ZERO
+               END IF
+               I = I - 1
+            END DO
+*
+*           Convert PERMUTATIONS
+*
+*           Apply permutaions to submatrices of upper part of A
+*           in factorization order where i decreases from N to 1
+*
+            I = N
+            DO WHILE ( I.GE.1 )
+               IF( IPIV( I ).GT.0 ) THEN
+*
+*                 1-by-1 pivot interchange
+*
+*                 Swap rows i and IPIV(i) in A(1:i,N-i:N)
+*
+                  IP = IPIV( I )
+                  IF( I.LT.N ) THEN
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( N-I, A( I, I+1 ), LDA,
+     $                              A( IP, I+1 ), LDA )
+                     END IF
+                  END IF
+*
+               ELSE
+*
+*                 2-by-2 pivot interchange
+*
+*                 Swap rows i and IPIV(i) and i-1 and IPIV(i-1)
+*                 in A(1:i,N-i:N)
+*
+                  IP = -IPIV( I )
+                  IP2 = -IPIV( I-1 )
+                  IF( I.LT.N ) THEN
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( N-I, A( I, I+1 ), LDA,
+     $                              A( IP, I+1 ), LDA )
+                     END IF
+                     IF( IP2.NE.(I-1) ) THEN
+                        CALL SSWAP( N-I, A( I-1, I+1 ), LDA,
+     $                              A( IP2, I+1 ), LDA )
+                     END IF
+                  END IF
+                  I = I - 1
+*
+               END IF
+               I = I - 1
+            END DO
+*
+         ELSE
+*
+*           Revert A (A is upper)
+*
+*
+*           Revert PERMUTATIONS
+*
+*           Apply permutaions to submatrices of upper part of A
+*           in reverse factorization order where i increases from 1 to N
+*
+            I = 1
+            DO WHILE ( I.LE.N )
+               IF( IPIV( I ).GT.0 ) THEN
+*
+*                 1-by-1 pivot interchange
+*
+*                 Swap rows i and IPIV(i) in A(1:i,N-i:N)
+*
+                  IP = IPIV( I )
+                  IF( I.LT.N ) THEN
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( N-I, A( IP, I+1 ), LDA,
+     $                              A( I, I+1 ), LDA )
+                     END IF
+                  END IF
+*
+               ELSE
+*
+*                 2-by-2 pivot interchange
+*
+*                 Swap rows i-1 and IPIV(i-1) and i and IPIV(i)
+*                 in A(1:i,N-i:N)
+*
+                  I = I + 1
+                  IP = -IPIV( I )
+                  IP2 = -IPIV( I-1 )
+                  IF( I.LT.N ) THEN
+                     IF( IP2.NE.(I-1) ) THEN
+                        CALL SSWAP( N-I, A( IP2, I+1 ), LDA,
+     $                              A( I-1, I+1 ), LDA )
+                     END IF
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( N-I, A( IP, I+1 ), LDA,
+     $                              A( I, I+1 ), LDA )
+                     END IF
+                  END IF
+*
+               END IF
+               I = I + 1
+            END DO
+*
+*           Revert VALUE
+*           Assign superdiagonal entries of D from array E to
+*           superdiagonal entries of A.
+*
+            I = N
+            DO WHILE ( I.GT.1 )
+               IF( IPIV( I ).LT.0 ) THEN
+                  A( I-1, I ) = E( I )
+                  I = I - 1
+               END IF
+               I = I - 1
+            END DO
+*
+*        End A is UPPER
+*
+         END IF
+*
+      ELSE
+*
+*        Begin A is LOWER
+*
+         IF ( CONVERT ) THEN
+*
+*           Convert A (A is lower)
+*
+*
+*           Convert VALUE
+*           Assign subdiagonal entries of D to array E and zero out
+*           corresponding entries in input storage A
+*
+            I = 1
+            E( N ) = ZERO
+            DO WHILE ( I.LE.N )
+               IF( I.LT.N .AND. IPIV(I).LT.0 ) THEN
+                  E( I ) = A( I+1, I )
+                  E( I+1 ) = ZERO
+                  A( I+1, I ) = ZERO
+                  I = I + 1
+               ELSE
+                  E( I ) = ZERO
+               END IF
+               I = I + 1
+            END DO
+*
+*           Convert PERMUTATIONS
+*
+*           Apply permutaions to submatrices of lower part of A
+*           in factorization order where i increases from 1 to N
+*
+            I = 1
+            DO WHILE ( I.LE.N )
+               IF( IPIV( I ).GT.0 ) THEN
+*
+*                 1-by-1 pivot interchange
+*
+*                 Swap rows i and IPIV(i) in A(i:N,1:i-1)
+*
+                  IP = IPIV( I )
+                  IF ( I.GT.1 ) THEN
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( I-1, A( I, 1 ), LDA,
+     $                              A( IP, 1 ), LDA )
+                     END IF
+                  END IF
+*
+               ELSE
+*
+*                 2-by-2 pivot interchange
+*
+*                 Swap rows i and IPIV(i) and i+1 and IPIV(i+1)
+*                 in A(i:N,1:i-1)
+*
+                  IP = -IPIV( I )
+                  IP2 = -IPIV( I+1 )
+                  IF ( I.GT.1 ) THEN
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( I-1, A( I, 1 ), LDA,
+     $                              A( IP, 1 ), LDA )
+                     END IF
+                     IF( IP2.NE.(I+1) ) THEN
+                        CALL SSWAP( I-1, A( I+1, 1 ), LDA,
+     $                              A( IP2, 1 ), LDA )
+                     END IF
+                  END IF
+                  I = I + 1
+*
+               END IF
+               I = I + 1
+            END DO
+*
+         ELSE
+*
+*           Revert A (A is lower)
+*
+*
+*           Revert PERMUTATIONS
+*
+*           Apply permutaions to submatrices of lower part of A
+*           in reverse factorization order where i decreases from N to 1
+*
+            I = N
+            DO WHILE ( I.GE.1 )
+               IF( IPIV( I ).GT.0 ) THEN
+*
+*                 1-by-1 pivot interchange
+*
+*                 Swap rows i and IPIV(i) in A(i:N,1:i-1)
+*
+                  IP = IPIV( I )
+                  IF ( I.GT.1 ) THEN
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( I-1, A( IP, 1 ), LDA,
+     $                              A( I, 1 ), LDA )
+                     END IF
+                  END IF
+*
+               ELSE
+*
+*                 2-by-2 pivot interchange
+*
+*                 Swap rows i+1 and IPIV(i+1) and i and IPIV(i)
+*                 in A(i:N,1:i-1)
+*
+                  I = I - 1
+                  IP = -IPIV( I )
+                  IP2 = -IPIV( I+1 )
+                  IF ( I.GT.1 ) THEN
+                     IF( IP2.NE.(I+1) ) THEN
+                        CALL SSWAP( I-1, A( IP2, 1 ), LDA,
+     $                              A( I+1, 1 ), LDA )
+                     END IF
+                     IF( IP.NE.I ) THEN
+                        CALL SSWAP( I-1, A( IP, 1 ), LDA,
+     $                              A( I, 1 ), LDA )
+                     END IF
+                  END IF
+*
+               END IF
+               I = I - 1
+            END DO
+*
+*           Revert VALUE
+*           Assign subdiagonal entries of D from array E to
+*           subgiagonal entries of A.
+*
+            I = 1
+            DO WHILE ( I.LE.N-1 )
+               IF( IPIV( I ).LT.0 ) THEN
+                  A( I + 1, I ) = E( I )
+                  I = I + 1
+               END IF
+               I = I + 1
+            END DO
+*
+         END IF
+*
+*        End A is LOWER
+*
+      END IF
+
+      RETURN
+*
+*     End of SSYCONVF_ROOK
+*
+      END