 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CLAGHE( 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 REAL D( * ) 00013 COMPLEX A( LDA, * ), WORK( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * CLAGHE generates a complex hermitian matrix A, by pre- and post- 00020 * multiplying a real diagonal matrix D with a random unitary matrix: 00021 * A = U*D*U'. The semi-bandwidth may then be reduced to k by additional 00022 * 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) REAL array, dimension (N) 00035 * The diagonal elements of the diagonal matrix D. 00036 * 00037 * A (output) COMPLEX array, dimension (LDA,N) 00038 * The generated n by n hermitian 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 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 ZERO, ONE, HALF 00060 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ), 00061 $ ONE = ( 1.0E+0, 0.0E+0 ), 00062 $ HALF = ( 0.5E+0, 0.0E+0 ) ) 00063 * .. 00064 * .. Local Scalars .. 00065 INTEGER I, J 00066 REAL WN 00067 COMPLEX ALPHA, TAU, WA, WB 00068 * .. 00069 * .. External Subroutines .. 00070 EXTERNAL CAXPY, CGEMV, CGERC, CHEMV, CHER2, CLARNV, 00071 $ CSCAL, XERBLA 00072 * .. 00073 * .. External Functions .. 00074 REAL SCNRM2 00075 COMPLEX CDOTC 00076 EXTERNAL SCNRM2, CDOTC 00077 * .. 00078 * .. Intrinsic Functions .. 00079 INTRINSIC ABS, CONJG, MAX, REAL 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( 'CLAGHE', -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 hermitian matrix 00110 * 00111 DO 40 I = N - 1, 1, -1 00112 * 00113 * generate random reflection 00114 * 00115 CALL CLARNV( 3, ISEED, N-I+1, WORK ) 00116 WN = SCNRM2( 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 CSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 00123 WORK( 1 ) = ONE 00124 TAU = REAL( 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 * u 00131 * 00132 CALL CHEMV( 'Lower', N-I+1, TAU, A( I, I ), LDA, WORK, 1, ZERO, 00133 $ WORK( N+1 ), 1 ) 00134 * 00135 * compute v := y - 1/2 * tau * ( y, u ) * u 00136 * 00137 ALPHA = -HALF*TAU*CDOTC( N-I+1, WORK( N+1 ), 1, WORK, 1 ) 00138 CALL CAXPY( N-I+1, ALPHA, WORK, 1, WORK( N+1 ), 1 ) 00139 * 00140 * apply the transformation as a rank-2 update to A(i:n,i:n) 00141 * 00142 CALL CHER2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1, 00143 $ A( I, I ), LDA ) 00144 40 CONTINUE 00145 * 00146 * Reduce number of subdiagonals to K 00147 * 00148 DO 60 I = 1, N - 1 - K 00149 * 00150 * generate reflection to annihilate A(k+i+1:n,i) 00151 * 00152 WN = SCNRM2( N-K-I+1, A( K+I, I ), 1 ) 00153 WA = ( WN / ABS( A( K+I, I ) ) )*A( K+I, I ) 00154 IF( WN.EQ.ZERO ) THEN 00155 TAU = ZERO 00156 ELSE 00157 WB = A( K+I, I ) + WA 00158 CALL CSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 ) 00159 A( K+I, I ) = ONE 00160 TAU = REAL( WB / WA ) 00161 END IF 00162 * 00163 * apply reflection to A(k+i:n,i+1:k+i-1) from the left 00164 * 00165 CALL CGEMV( 'Conjugate transpose', N-K-I+1, K-1, ONE, 00166 $ A( K+I, I+1 ), LDA, A( K+I, I ), 1, ZERO, WORK, 1 ) 00167 CALL CGERC( N-K-I+1, K-1, -TAU, A( K+I, I ), 1, WORK, 1, 00168 $ A( K+I, I+1 ), LDA ) 00169 * 00170 * apply reflection to A(k+i:n,k+i:n) from the left and the right 00171 * 00172 * compute y := tau * A * u 00173 * 00174 CALL CHEMV( 'Lower', N-K-I+1, TAU, A( K+I, K+I ), LDA, 00175 $ A( K+I, I ), 1, ZERO, WORK, 1 ) 00176 * 00177 * compute v := y - 1/2 * tau * ( y, u ) * u 00178 * 00179 ALPHA = -HALF*TAU*CDOTC( N-K-I+1, WORK, 1, A( K+I, I ), 1 ) 00180 CALL CAXPY( N-K-I+1, ALPHA, A( K+I, I ), 1, WORK, 1 ) 00181 * 00182 * apply hermitian rank-2 update to A(k+i:n,k+i:n) 00183 * 00184 CALL CHER2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1, 00185 $ A( K+I, K+I ), LDA ) 00186 * 00187 A( K+I, I ) = -WA 00188 DO 50 J = K + I + 1, N 00189 A( J, I ) = ZERO 00190 50 CONTINUE 00191 60 CONTINUE 00192 * 00193 * Store full hermitian matrix 00194 * 00195 DO 80 J = 1, N 00196 DO 70 I = J + 1, N 00197 A( J, I ) = CONJG( A( I, J ) ) 00198 70 CONTINUE 00199 80 CONTINUE 00200 RETURN 00201 * 00202 * End of CLAGHE 00203 * 00204 END
1.7.3 
