144 lines
		
	
	
		
			4.2 KiB
		
	
	
	
		
			Fortran
		
	
	
	
			
		
		
	
	
			144 lines
		
	
	
		
			4.2 KiB
		
	
	
	
		
			Fortran
		
	
	
	
      SUBROUTINE CLAUU2F( UPLO, N, A, LDA, INFO )
 | 
						|
*
 | 
						|
*  -- LAPACK auxiliary routine (version 3.1) --
 | 
						|
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 | 
						|
*     November 2006
 | 
						|
*
 | 
						|
*     .. Scalar Arguments ..
 | 
						|
      CHARACTER          UPLO
 | 
						|
      INTEGER            INFO, LDA, N
 | 
						|
*     ..
 | 
						|
*     .. Array Arguments ..
 | 
						|
      COMPLEX            A( LDA, * )
 | 
						|
*     ..
 | 
						|
*
 | 
						|
*  Purpose
 | 
						|
*  =======
 | 
						|
*
 | 
						|
*  CLAUU2 computes the product U * U' or L' * L, where the triangular
 | 
						|
*  factor U or L is stored in the upper or lower triangular part of
 | 
						|
*  the array A.
 | 
						|
*
 | 
						|
*  If UPLO = 'U' or 'u' then the upper triangle of the result is stored,
 | 
						|
*  overwriting the factor U in A.
 | 
						|
*  If UPLO = 'L' or 'l' then the lower triangle of the result is stored,
 | 
						|
*  overwriting the factor L in A.
 | 
						|
*
 | 
						|
*  This is the unblocked form of the algorithm, calling Level 2 BLAS.
 | 
						|
*
 | 
						|
*  Arguments
 | 
						|
*  =========
 | 
						|
*
 | 
						|
*  UPLO    (input) CHARACTER*1
 | 
						|
*          Specifies whether the triangular factor stored in the array A
 | 
						|
*          is upper or lower triangular:
 | 
						|
*          = 'U':  Upper triangular
 | 
						|
*          = 'L':  Lower triangular
 | 
						|
*
 | 
						|
*  N       (input) INTEGER
 | 
						|
*          The order of the triangular factor U or L.  N >= 0.
 | 
						|
*
 | 
						|
*  A       (input/output) COMPLEX array, dimension (LDA,N)
 | 
						|
*          On entry, the triangular factor U or L.
 | 
						|
*          On exit, if UPLO = 'U', the upper triangle of A is
 | 
						|
*          overwritten with the upper triangle of the product U * U';
 | 
						|
*          if UPLO = 'L', the lower triangle of A is overwritten with
 | 
						|
*          the lower triangle of the product L' * L.
 | 
						|
*
 | 
						|
*  LDA     (input) INTEGER
 | 
						|
*          The leading dimension of the array A.  LDA >= max(1,N).
 | 
						|
*
 | 
						|
*  INFO    (output) INTEGER
 | 
						|
*          = 0: successful exit
 | 
						|
*          < 0: if INFO = -k, the k-th argument had an illegal value
 | 
						|
*
 | 
						|
*  =====================================================================
 | 
						|
*
 | 
						|
*     .. Parameters ..
 | 
						|
      COMPLEX            ONE
 | 
						|
      PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
 | 
						|
*     ..
 | 
						|
*     .. Local Scalars ..
 | 
						|
      LOGICAL            UPPER
 | 
						|
      INTEGER            I
 | 
						|
      REAL               AII
 | 
						|
*     ..
 | 
						|
*     .. External Functions ..
 | 
						|
      LOGICAL            LSAME
 | 
						|
      COMPLEX            CDOTC
 | 
						|
      EXTERNAL           LSAME, CDOTC
 | 
						|
*     ..
 | 
						|
*     .. External Subroutines ..
 | 
						|
      EXTERNAL           CGEMV, CLACGV, CSSCAL, XERBLA
 | 
						|
*     ..
 | 
						|
*     .. Intrinsic Functions ..
 | 
						|
      INTRINSIC          CMPLX, MAX, REAL
 | 
						|
*     ..
 | 
						|
*     .. Executable Statements ..
 | 
						|
*
 | 
						|
*     Test the input parameters.
 | 
						|
*
 | 
						|
      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( LDA.LT.MAX( 1, N ) ) THEN
 | 
						|
         INFO = -4
 | 
						|
      END IF
 | 
						|
      IF( INFO.NE.0 ) THEN
 | 
						|
         CALL XERBLA( 'CLAUU2', -INFO )
 | 
						|
         RETURN
 | 
						|
      END IF
 | 
						|
*
 | 
						|
*     Quick return if possible
 | 
						|
*
 | 
						|
      IF( N.EQ.0 )
 | 
						|
     $   RETURN
 | 
						|
*
 | 
						|
      IF( UPPER ) THEN
 | 
						|
*
 | 
						|
*        Compute the product U * U'.
 | 
						|
*
 | 
						|
         DO 10 I = 1, N
 | 
						|
            AII = A( I, I )
 | 
						|
            IF( I.LT.N ) THEN
 | 
						|
               A( I, I ) = AII*AII + REAL( CDOTC( N-I, A( I, I+1 ), LDA,
 | 
						|
     $                     A( I, I+1 ), LDA ) )
 | 
						|
               CALL CLACGV( N-I, A( I, I+1 ), LDA )
 | 
						|
               CALL CGEMV( 'No transpose', I-1, N-I, ONE, A( 1, I+1 ),
 | 
						|
     $                     LDA, A( I, I+1 ), LDA, CMPLX( AII ),
 | 
						|
     $                     A( 1, I ), 1 )
 | 
						|
               CALL CLACGV( N-I, A( I, I+1 ), LDA )
 | 
						|
            ELSE
 | 
						|
               CALL CSSCAL( I, AII, A( 1, I ), 1 )
 | 
						|
            END IF
 | 
						|
   10    CONTINUE
 | 
						|
*
 | 
						|
      ELSE
 | 
						|
*
 | 
						|
*        Compute the product L' * L.
 | 
						|
*
 | 
						|
         DO 20 I = 1, N
 | 
						|
            AII = A( I, I )
 | 
						|
            IF( I.LT.N ) THEN
 | 
						|
               A( I, I ) = AII*AII + REAL( CDOTC( N-I, A( I+1, I ), 1,
 | 
						|
     $                     A( I+1, I ), 1 ) )
 | 
						|
               CALL CLACGV( I-1, A( I, 1 ), LDA )
 | 
						|
               CALL CGEMV( 'Conjugate transpose', N-I, I-1, ONE,
 | 
						|
     $                     A( I+1, 1 ), LDA, A( I+1, I ), 1,
 | 
						|
     $                     CMPLX( AII ), A( I, 1 ), LDA )
 | 
						|
               CALL CLACGV( I-1, A( I, 1 ), LDA )
 | 
						|
            ELSE
 | 
						|
               CALL CSSCAL( I, AII, A( I, 1 ), LDA )
 | 
						|
            END IF
 | 
						|
   20    CONTINUE
 | 
						|
      END IF
 | 
						|
*
 | 
						|
      RETURN
 | 
						|
*
 | 
						|
*     End of CLAUU2
 | 
						|
*
 | 
						|
      END
 |