001:       SUBROUTINE ZLAGHE( N, K, D, A, LDA, ISEED, WORK, INFO )
002: *
003: *  -- LAPACK auxiliary test routine (version 3.1) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       INTEGER            INFO, K, LDA, N
009: *     ..
010: *     .. Array Arguments ..
011:       INTEGER            ISEED( 4 )
012:       DOUBLE PRECISION   D( * )
013:       COMPLEX*16         A( LDA, * ), WORK( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  ZLAGHE generates a complex hermitian matrix A, by pre- and post-
020: *  multiplying a real diagonal matrix D with a random unitary matrix:
021: *  A = U*D*U''. The semi-bandwidth may then be reduced to k by additional
022: *  unitary transformations.
023: *
024: *  Arguments
025: *  =========
026: *
027: *  N       (input) INTEGER
028: *          The order of the matrix A.  N >= 0.
029: *
030: *  K       (input) INTEGER
031: *          The number of nonzero subdiagonals within the band of A.
032: *          0 <= K <= N-1.
033: *
034: *  D       (input) DOUBLE PRECISION array, dimension (N)
035: *          The diagonal elements of the diagonal matrix D.
036: *
037: *  A       (output) COMPLEX*16 array, dimension (LDA,N)
038: *          The generated n by n hermitian matrix A (the full matrix is
039: *          stored).
040: *
041: *  LDA     (input) INTEGER
042: *          The leading dimension of the array A.  LDA >= N.
043: *
044: *  ISEED   (input/output) INTEGER array, dimension (4)
045: *          On entry, the seed of the random number generator; the array
046: *          elements must be between 0 and 4095, and ISEED(4) must be
047: *          odd.
048: *          On exit, the seed is updated.
049: *
050: *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
051: *
052: *  INFO    (output) INTEGER
053: *          = 0: successful exit
054: *          < 0: if INFO = -i, the i-th argument had an illegal value
055: *
056: *  =====================================================================
057: *
058: *     .. Parameters ..
059:       COMPLEX*16         ZERO, ONE, HALF
060:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
061:      $                   ONE = ( 1.0D+0, 0.0D+0 ),
062:      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
063: *     ..
064: *     .. Local Scalars ..
065:       INTEGER            I, J
066:       DOUBLE PRECISION   WN
067:       COMPLEX*16         ALPHA, TAU, WA, WB, TMP
068: *     ..
069: *     .. External Subroutines ..
070:       EXTERNAL           XERBLA, ZAXPY, ZGEMV, ZGERC, ZHEMV, ZHER2,
071:      $                   ZLARNV, ZSCAL
072: *     ..
073: *     .. External Functions ..
074:       DOUBLE PRECISION   DZNRM2
075: *      COMPLEX*16         ZDOTC
076: *      EXTERNAL           DZNRM2, ZDOTC
077:       EXTERNAL           DZNRM2
078: *     ..
079: *     .. Intrinsic Functions ..
080:       INTRINSIC          ABS, DBLE, DCONJG, MAX
081: *     ..
082: *     .. Executable Statements ..
083: *
084: *     Test the input arguments
085: *
086:       INFO = 0
087:       IF( N.LT.0 ) THEN
088:          INFO = -1
089:       ELSE IF( K.LT.0 .OR. K.GT.N-1 ) THEN
090:          INFO = -2
091:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
092:          INFO = -5
093:       END IF
094:       IF( INFO.LT.0 ) THEN
095:          CALL XERBLA( 'ZLAGHE', -INFO )
096:          RETURN
097:       END IF
098: *
099: *     initialize lower triangle of A to diagonal matrix
100: *
101:       DO 20 J = 1, N
102:          DO 10 I = J + 1, N
103:             A( I, J ) = ZERO
104:    10    CONTINUE
105:    20 CONTINUE
106:       DO 30 I = 1, N
107:          A( I, I ) = D( I )
108:    30 CONTINUE
109: *
110: *     Generate lower triangle of hermitian matrix
111: *
112:       DO 40 I = N - 1, 1, -1
113: *
114: *        generate random reflection
115: *
116:          CALL ZLARNV( 3, ISEED, N-I+1, WORK )
117:          WN = DZNRM2( N-I+1, WORK, 1 )
118:          WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 )
119:          IF( WN.EQ.ZERO ) THEN
120:             TAU = ZERO
121:          ELSE
122:             WB = WORK( 1 ) + WA
123:             CALL ZSCAL( N-I, ONE / WB, WORK( 2 ), 1 )
124:             WORK( 1 ) = ONE
125:             TAU = DBLE( WB / WA )
126:          END IF
127: *
128: *        apply random reflection to A(i:n,i:n) from the left
129: *        and the right
130: *
131: *        compute  y := tau * A * u
132: *
133:          CALL ZHEMV( 'Lower', N-I+1, TAU, A( I, I ), LDA, WORK, 1, ZERO,
134:      $               WORK( N+1 ), 1 )
135: *
136: *        compute  v := y - 1/2 * tau * ( y, u ) * u
137: *
138:          CALL PLASMA_ZDOTCSUB( N-I+1, WORK( N+1 ), 1, WORK, 1, TMP )
139:          ALPHA = -HALF*TAU*TMP
140:          CALL ZAXPY( N-I+1, ALPHA, WORK, 1, WORK( N+1 ), 1 )
141: *
142: *        apply the transformation as a rank-2 update to A(i:n,i:n)
143: *
144:          CALL ZHER2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1,
145:      $               A( I, I ), LDA )
146:    40 CONTINUE
147: *
148: *     Reduce number of subdiagonals to K
149: *
150:       DO 60 I = 1, N - 1 - K
151: *
152: *        generate reflection to annihilate A(k+i+1:n,i)
153: *
154:          WN = DZNRM2( N-K-I+1, A( K+I, I ), 1 )
155:          WA = ( WN / ABS( A( K+I, I ) ) )*A( K+I, I )
156:          IF( WN.EQ.ZERO ) THEN
157:             TAU = ZERO
158:          ELSE
159:             WB = A( K+I, I ) + WA
160:             CALL ZSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 )
161:             A( K+I, I ) = ONE
162:             TAU = DBLE( WB / WA )
163:          END IF
164: *
165: *        apply reflection to A(k+i:n,i+1:k+i-1) from the left
166: *
167:          CALL ZGEMV( 'Conjugate transpose', N-K-I+1, K-1, ONE,
168:      $               A( K+I, I+1 ), LDA, A( K+I, I ), 1, ZERO, WORK, 1 )
169:          CALL ZGERC( N-K-I+1, K-1, -TAU, A( K+I, I ), 1, WORK, 1,
170:      $               A( K+I, I+1 ), LDA )
171: *
172: *        apply reflection to A(k+i:n,k+i:n) from the left and the right
173: *
174: *        compute  y := tau * A * u
175: *
176:          CALL ZHEMV( 'Lower', N-K-I+1, TAU, A( K+I, K+I ), LDA,
177:      $               A( K+I, I ), 1, ZERO, WORK, 1 )
178: *
179: *        compute  v := y - 1/2 * tau * ( y, u ) * u
180: *
181:          CALL PLASMA_ZDOTCSUB( N-K-I+1, WORK, 1, A( K+I, I ), 1, TMP )
182:          ALPHA = -HALF*TAU*TMP
183:          CALL ZAXPY( N-K-I+1, ALPHA, A( K+I, I ), 1, WORK, 1 )
184: *
185: *        apply hermitian rank-2 update to A(k+i:n,k+i:n)
186: *
187:          CALL ZHER2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1,
188:      $               A( K+I, K+I ), LDA )
189: *
190:          A( K+I, I ) = -WA
191:          DO 50 J = K + I + 1, N
192:             A( J, I ) = ZERO
193:    50    CONTINUE
194:    60 CONTINUE
195: *
196: *     Store full hermitian matrix
197: *
198:       DO 80 J = 1, N
199:          DO 70 I = J + 1, N
200:             A( J, I ) = DCONJG( A( I, J ) )
201:    70    CONTINUE
202:    80 CONTINUE
203:       RETURN
204: *
205: *     End of ZLAGHE
206: *
207:       END
208: