|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZLAGSY( N, K, D, A, LDA, ISEED, WORK, INFO ) 00002 * 00003 * -- LAPACK auxiliary test routine (version 3.1) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2006 00006 * 00007 * .. Scalar Arguments .. 00008 INTEGER INFO, K, LDA, N 00009 * .. 00010 * .. Array Arguments .. 00011 INTEGER ISEED( 4 ) 00012 DOUBLE PRECISION D( * ) 00013 COMPLEX*16 A( LDA, * ), WORK( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * ZLAGSY generates a complex symmetric matrix A, by pre- and post- 00020 * multiplying a real diagonal matrix D with a random unitary matrix: 00021 * A = U*D*U**T. The semi-bandwidth may then be reduced to k by 00022 * additional unitary transformations. 00023 * 00024 * Arguments 00025 * ========= 00026 * 00027 * N (input) INTEGER 00028 * The order of the matrix A. N >= 0. 00029 * 00030 * K (input) INTEGER 00031 * The number of nonzero subdiagonals within the band of A. 00032 * 0 <= K <= N-1. 00033 * 00034 * D (input) DOUBLE PRECISION array, dimension (N) 00035 * The diagonal elements of the diagonal matrix D. 00036 * 00037 * A (output) COMPLEX*16 array, dimension (LDA,N) 00038 * The generated n by n symmetric matrix A (the full matrix is 00039 * stored). 00040 * 00041 * LDA (input) INTEGER 00042 * The leading dimension of the array A. LDA >= N. 00043 * 00044 * ISEED (input/output) INTEGER array, dimension (4) 00045 * On entry, the seed of the random number generator; the array 00046 * elements must be between 0 and 4095, and ISEED(4) must be 00047 * odd. 00048 * On exit, the seed is updated. 00049 * 00050 * WORK (workspace) COMPLEX*16 array, dimension (2*N) 00051 * 00052 * INFO (output) INTEGER 00053 * = 0: successful exit 00054 * < 0: if INFO = -i, the i-th argument had an illegal value 00055 * 00056 * ===================================================================== 00057 * 00058 * .. Parameters .. 00059 COMPLEX*16 ZERO, ONE, HALF 00060 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 00061 $ ONE = ( 1.0D+0, 0.0D+0 ), 00062 $ HALF = ( 0.5D+0, 0.0D+0 ) ) 00063 * .. 00064 * .. Local Scalars .. 00065 INTEGER I, II, J, JJ 00066 DOUBLE PRECISION WN 00067 COMPLEX*16 ALPHA, TAU, WA, WB 00068 * .. 00069 * .. External Subroutines .. 00070 EXTERNAL XERBLA, ZAXPY, ZGEMV, ZGERC, ZLACGV, ZLARNV, 00071 $ ZSCAL, ZSYMV 00072 * .. 00073 * .. External Functions .. 00074 DOUBLE PRECISION DZNRM2 00075 COMPLEX*16 ZDOTC 00076 EXTERNAL DZNRM2, ZDOTC 00077 * .. 00078 * .. Intrinsic Functions .. 00079 INTRINSIC ABS, DBLE, MAX 00080 * .. 00081 * .. Executable Statements .. 00082 * 00083 * Test the input arguments 00084 * 00085 INFO = 0 00086 IF( N.LT.0 ) THEN 00087 INFO = -1 00088 ELSE IF( K.LT.0 .OR. K.GT.N-1 ) THEN 00089 INFO = -2 00090 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00091 INFO = -5 00092 END IF 00093 IF( INFO.LT.0 ) THEN 00094 CALL XERBLA( 'ZLAGSY', -INFO ) 00095 RETURN 00096 END IF 00097 * 00098 * initialize lower triangle of A to diagonal matrix 00099 * 00100 DO 20 J = 1, N 00101 DO 10 I = J + 1, N 00102 A( I, J ) = ZERO 00103 10 CONTINUE 00104 20 CONTINUE 00105 DO 30 I = 1, N 00106 A( I, I ) = D( I ) 00107 30 CONTINUE 00108 * 00109 * Generate lower triangle of symmetric matrix 00110 * 00111 DO 60 I = N - 1, 1, -1 00112 * 00113 * generate random reflection 00114 * 00115 CALL ZLARNV( 3, ISEED, N-I+1, WORK ) 00116 WN = DZNRM2( N-I+1, WORK, 1 ) 00117 WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 ) 00118 IF( WN.EQ.ZERO ) THEN 00119 TAU = ZERO 00120 ELSE 00121 WB = WORK( 1 ) + WA 00122 CALL ZSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 00123 WORK( 1 ) = ONE 00124 TAU = DBLE( WB / WA ) 00125 END IF 00126 * 00127 * apply random reflection to A(i:n,i:n) from the left 00128 * and the right 00129 * 00130 * compute y := tau * A * conjg(u) 00131 * 00132 CALL ZLACGV( N-I+1, WORK, 1 ) 00133 CALL ZSYMV( 'Lower', N-I+1, TAU, A( I, I ), LDA, WORK, 1, ZERO, 00134 $ WORK( N+1 ), 1 ) 00135 CALL ZLACGV( N-I+1, WORK, 1 ) 00136 * 00137 * compute v := y - 1/2 * tau * ( u, y ) * u 00138 * 00139 ALPHA = -HALF*TAU*ZDOTC( N-I+1, WORK, 1, WORK( N+1 ), 1 ) 00140 CALL ZAXPY( N-I+1, ALPHA, WORK, 1, WORK( N+1 ), 1 ) 00141 * 00142 * apply the transformation as a rank-2 update to A(i:n,i:n) 00143 * 00144 * CALL ZSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1, 00145 * $ A( I, I ), LDA ) 00146 * 00147 DO 50 JJ = I, N 00148 DO 40 II = JJ, N 00149 A( II, JJ ) = A( II, JJ ) - 00150 $ WORK( II-I+1 )*WORK( N+JJ-I+1 ) - 00151 $ WORK( N+II-I+1 )*WORK( JJ-I+1 ) 00152 40 CONTINUE 00153 50 CONTINUE 00154 60 CONTINUE 00155 * 00156 * Reduce number of subdiagonals to K 00157 * 00158 DO 100 I = 1, N - 1 - K 00159 * 00160 * generate reflection to annihilate A(k+i+1:n,i) 00161 * 00162 WN = DZNRM2( N-K-I+1, A( K+I, I ), 1 ) 00163 WA = ( WN / ABS( A( K+I, I ) ) )*A( K+I, I ) 00164 IF( WN.EQ.ZERO ) THEN 00165 TAU = ZERO 00166 ELSE 00167 WB = A( K+I, I ) + WA 00168 CALL ZSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 ) 00169 A( K+I, I ) = ONE 00170 TAU = DBLE( WB / WA ) 00171 END IF 00172 * 00173 * apply reflection to A(k+i:n,i+1:k+i-1) from the left 00174 * 00175 CALL ZGEMV( 'Conjugate transpose', N-K-I+1, K-1, ONE, 00176 $ A( K+I, I+1 ), LDA, A( K+I, I ), 1, ZERO, WORK, 1 ) 00177 CALL ZGERC( N-K-I+1, K-1, -TAU, A( K+I, I ), 1, WORK, 1, 00178 $ A( K+I, I+1 ), LDA ) 00179 * 00180 * apply reflection to A(k+i:n,k+i:n) from the left and the right 00181 * 00182 * compute y := tau * A * conjg(u) 00183 * 00184 CALL ZLACGV( N-K-I+1, A( K+I, I ), 1 ) 00185 CALL ZSYMV( 'Lower', N-K-I+1, TAU, A( K+I, K+I ), LDA, 00186 $ A( K+I, I ), 1, ZERO, WORK, 1 ) 00187 CALL ZLACGV( N-K-I+1, A( K+I, I ), 1 ) 00188 * 00189 * compute v := y - 1/2 * tau * ( u, y ) * u 00190 * 00191 ALPHA = -HALF*TAU*ZDOTC( N-K-I+1, A( K+I, I ), 1, WORK, 1 ) 00192 CALL ZAXPY( N-K-I+1, ALPHA, A( K+I, I ), 1, WORK, 1 ) 00193 * 00194 * apply symmetric rank-2 update to A(k+i:n,k+i:n) 00195 * 00196 * CALL ZSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1, 00197 * $ A( K+I, K+I ), LDA ) 00198 * 00199 DO 80 JJ = K + I, N 00200 DO 70 II = JJ, N 00201 A( II, JJ ) = A( II, JJ ) - A( II, I )*WORK( JJ-K-I+1 ) - 00202 $ WORK( II-K-I+1 )*A( JJ, I ) 00203 70 CONTINUE 00204 80 CONTINUE 00205 * 00206 A( K+I, I ) = -WA 00207 DO 90 J = K + I + 1, N 00208 A( J, I ) = ZERO 00209 90 CONTINUE 00210 100 CONTINUE 00211 * 00212 * Store full symmetric matrix 00213 * 00214 DO 120 J = 1, N 00215 DO 110 I = J + 1, N 00216 A( J, I ) = A( I, J ) 00217 110 CONTINUE 00218 120 CONTINUE 00219 RETURN 00220 * 00221 * End of ZLAGSY 00222 * 00223 END
1.7.3 
