SUBROUTINE CGESVD_F95( A, S, U, VT, WW, JOB, INFO ) ! ! -- LAPACK95 interface driver routine (version 3.0) -- ! UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK ! September, 2000 ! ! .. USE STATEMENTS .. USE LA_PRECISION, ONLY: WP => SP USE LA_AUXMOD, ONLY: ERINFO, LSAME USE F77_LAPACK, ONLY: GESVD_F77 => LA_GESVD ! .. IMPLICIT STATEMENT .. IMPLICIT NONE ! .. SCALAR ARGUMENTS .. CHARACTER(LEN=1), OPTIONAL, INTENT(IN) :: JOB INTEGER, INTENT(OUT), OPTIONAL :: INFO ! .. ARRAY ARGUMENTS .. COMPLEX(WP), INTENT(INOUT) :: A(:,:) REAL(WP), INTENT(OUT) :: S(:) REAL(WP), INTENT(OUT), OPTIONAL :: WW(:) COMPLEX(WP), INTENT(OUT), OPTIONAL, TARGET :: U(:,:), VT(:,:) !---------------------------------------------------------------------- ! ! Purpose ! ======= ! ! LA_GESVD and LA_GESDD compute the singular values and, ! optionally, the left and/or right singular vectors from the singular ! value decomposition (SVD) of a real or complex m by n matrix A. The ! SVD of A is written ! A = U * SIGMA * V^H ! where SIGMA is an m by n matrix which is zero except for its ! min(m, n) diagonal elements, U is an m by m orthogonal (unitary) ! matrix, and V is an n by n orthogonal (unitary) matrix. The diagonal ! elements of SIGMA , i.e., the values ! ! sigma(i)= SIGMA(i,i), i = 1, 2,..., min(m, n) ! are the singular values of A; they are real and non-negative, and are ! returned in descending order. The first min(m, n) columns of U and V ! are the left and right singular vectors of A, respectively. ! LA_GESDD solves the same problem as LA_GESVD but uses a divide and ! conquer method if singular vectors are desired. For large matrices it ! is usually much faster than LA_GESVD when singular vectors are ! desired, but uses more workspace. ! ! Note: The routine returns V^H , not V . ! ! ======== ! ! SUBROUTINE LA_GESVD / LA_GESDD( A, S, U=u, VT=vt, & ! WW=ww, JOB=job, INFO=info ) ! (), INTENT(INOUT) :: A(:,:) ! REAL(), INTENT(OUT) :: S(:) ! (), INTENT(OUT), OPTIONAL :: U(:,:), VT(:,:) ! REAL(), INTENT(OUT), OPTIONAL :: WW(:) ! CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: JOB ! INTEGER, INTENT(OUT), OPTIONAL :: INFO ! where ! ::= REAL | COMPLEX ! ::= KIND(1.0) | KIND(1.0D0) ! ! Arguments ! ========= ! ! A (input/output) REAL or COMPLEX array, shape (:, :) with ! size(A, 1) = m and size(A, 2) = n. ! On entry, the matrix A. ! On exit, if JOB = 'U' and U is not present, then A is ! overwritten with the first min(m, n) columns of U (the left ! singular vectors, stored columnwise). ! If JOB = 'V' and VT is not present, then A is overwritten with ! the first min(m, n) rows of V^H (the right singular vectors, ! stored rowwise). ! In all cases the original contents of A are destroyed. ! S (output) REAL array, shape (:) with size(S) = min(m, n). ! The singular values of A, sorted so that S(i) >= S(i+1). ! U Optional (output) REAL or COMPLEX array, shape (:, :) with ! size(U, 1) = m and size(U, 2) = m or min(m, n). ! If size(U, 2) = m, U contains the m by m matrix U . ! If size(U; 2) = min(m, n), U contains the first min(m, n) ! columns of U (the left singular vectors, stored columnwise). ! VT Optional (output) REAL or COMPLEX array, shape (:, :) with ! size(VT, 1) = n or min(m, n) and size(VT, 2) = n. ! If size(VT, 1) = n , VT contains the n by n matrix V^H . ! If size(VT, 1) = min(m, n), VT contains the first min(m, n) ! rows of V^H (the right singular vectors, stored rowwise). ! WW Optional (output) REAL array, shape (:) with size(WW) = ! min(m, n) - 1 ! If INFO > 0, WW contains the unconverged superdiagonal elements ! of an upper bidiagonal matrix B whose diagonal is in SIGMA (not ! necessarily sorted). B has the same singular values as A. ! Note: WW is a dummy argument for LA_GESDD. ! JOB Optional (input) CHARACTER(LEN=1). ! = 'N': neither columns of U nor rows of V^H are returned in ! array A. ! = 'U': if U is not present, the first min(m, n) columns of U ! (the left singular vectors) are returned in array A; ! = 'V': if VT is not present, the first min(m, n) rows of V^H ! (the right singular vectors) are returned in array A; ! Default value: 'N'. ! INFO Optional (output) INTEGER. ! = 0: successful exit. ! < 0: if INFO = -i, the i-th argument had an illegal value. ! > 0: The algorithm did not converge. ! If INFO is not present and an error occurs, then the program is ! terminated with an error message. !---------------------------------------------------------------------- ! .. LOCAL PARAMETERS .. CHARACTER(LEN=8), PARAMETER :: SRNAME = 'LA_GESVD' ! .. LOCAL SCALARS .. CHARACTER(LEN=1) :: LJOB CHARACTER(LEN=1) :: LJOBU, LJOBVT INTEGER, SAVE :: LWORK = 0 INTEGER :: N, M, LINFO, LD, ISTAT, ISTAT1, S1U, S2U, S1VT, S2VT, & NN, MN, SWW ! .. LOCAL ARRAYS .. COMPLEX(WP), TARGET :: LLU(1,1), LLVT(1,1) COMPLEX(WP), POINTER :: WORK(:) REAL(WP), POINTER :: RWORK(:) COMPLEX(WP) :: WORKMIN(1) ! .. INTRINSIC FUNCTIONS .. INTRINSIC MIN, MAX, PRESENT, SIZE ! .. EXECUTABLE STATEMENTS .. LINFO = 0; ISTAT = 0; M = SIZE(A,1); N = SIZE(A,2) LD = MAX(1,M); MN = MIN(M,N) IF( PRESENT(JOB) )THEN; LJOB = JOB; ELSE; LJOB = 'N'; ENDIF IF( PRESENT(U) )THEN; S1U = SIZE(U,1); S2U = SIZE(U,2) ELSE; S1U = 1; S2U = 1; END IF IF( PRESENT(VT) )THEN; S1VT = SIZE(VT,1); S2VT = SIZE(VT,2) ELSE; S1VT = 1; S2VT = 1; END IF IF( PRESENT(WW) )THEN; SWW = SIZE(WW); ELSE; SWW = MN-1; ENDIF ! .. TEST THE ARGUMENTS IF( M < 0 .OR. N < 0 )THEN; LINFO = -1 ELSE IF( SIZE( S ) /= MN )THEN; LINFO = -2 ELSE IF( PRESENT(U) .AND. ( S1U /= M .OR. & ( S2U /= M .AND. S2U /= MN ) ) )THEN; LINFO = -3 ELSE IF( PRESENT(VT) .AND. ( ( S1VT /= N .AND. S1VT /= MN ) & .OR. S2VT /= N ) )THEN; LINFO = -4 ELSE IF( SWW /= MN-1 .AND. MN > 0 ) THEN; LINFO = -5 ELSE IF( PRESENT(JOB) .AND. ( .NOT. ( LSAME(LJOB,'U') .OR. & LSAME(LJOB,'V') .OR. LSAME(LJOB,'N') ) .OR. & LSAME(LJOB,'U') .AND. PRESENT(U) .OR. & LSAME(LJOB,'V') .AND. PRESENT(VT)) )THEN; LINFO = -6 ELSE IF( PRESENT(U) )THEN IF( S2U == M )THEN; LJOBU = 'A'; ELSE; LJOBU = 'S'; ENDIF ELSE; IF( LSAME(LJOB,'U') ) THEN; LJOBU = 'O' ELSE; LJOBU = 'N'; ENDIF ENDIF IF( PRESENT(VT) )THEN IF( S1VT == N )THEN; LJOBVT = 'A'; ELSE; LJOBVT = 'S'; ENDIF ELSE; IF( LSAME(LJOB,'V') )THEN; LJOBVT = 'O' ELSE; LJOBVT = 'N'; ENDIF ENDIF ALLOCATE( RWORK( MAX(1, 5*MIN(M,N))), STAT=ISTAT ) IF( ISTAT == 0 )THEN LWORK = -1 IF( PRESENT(U) ) THEN IF ( PRESENT(VT) )THEN CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), & & VT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO ) ELSE CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), & & LLVT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO ) ENDIF ELSE IF ( PRESENT(VT) )THEN CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), & & VT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO ) ELSE CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), & & LLVT, MAX(1,S1VT), WORKMIN, LWORK, RWORK, LINFO ) ENDIF ENDIF LWORK = WORKMIN(1) ! THE NEXT LINE SHOULD BE REMOVED LWORK = LWORK + 1 ALLOCATE(WORK(LWORK), STAT = ISTAT) IF( ISTAT /= 0 )THEN DEALLOCATE(WORK,STAT=ISTAT1) NN = 2*MIN(M,N)+MAX(M,N) LWORK = MAX( 1, NN); ALLOCATE(WORK(LWORK), STAT=ISTAT) IF( ISTAT == 0) CALL ERINFO( -200, SRNAME, LINFO ) END IF END IF IF( ISTAT == 0 ) THEN IF( PRESENT(U) ) THEN IF ( PRESENT(VT) )THEN CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), & VT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO ) ELSE CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, U, MAX(1,S1U), & & LLVT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO ) ENDIF ELSE IF ( PRESENT(VT) )THEN CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), & & VT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO ) ELSE CALL GESVD_F77( LJOBU, LJOBVT, M, N, A, LD, S, LLU, MAX(1,S1U), & & LLVT, MAX(1,S1VT), WORK, LWORK, RWORK, LINFO ) ENDIF ENDIF LWORK = INT(WORK(1)+1) IF( LINFO > 0 .AND. PRESENT(WW) ) WW(1:MN-1) = RWORK(1:MN-1) ELSE; LINFO = -100; ENDIF DEALLOCATE(WORK, RWORK, STAT=ISTAT1) ENDIF CALL ERINFO(LINFO,SRNAME,INFO,ISTAT) END SUBROUTINE CGESVD_F95