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