diff options
Diffstat (limited to 'SRC/ssyconv.f')
| -rw-r--r-- | SRC/ssyconv.f | 302 |
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 |