![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAGSY 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE SLAGSY( N, K, D, A, LDA, ISEED, WORK, INFO ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER INFO, K, LDA, N 00015 * .. 00016 * .. Array Arguments .. 00017 * INTEGER ISEED( 4 ) 00018 * REAL A( LDA, * ), D( * ), WORK( * ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> SLAGSY generates a real symmetric matrix A, by pre- and post- 00028 *> multiplying a real diagonal matrix D with a random orthogonal matrix: 00029 *> A = U*D*U'. The semi-bandwidth may then be reduced to k by additional 00030 *> orthogonal transformations. 00031 *> \endverbatim 00032 * 00033 * Arguments: 00034 * ========== 00035 * 00036 *> \param[in] N 00037 *> \verbatim 00038 *> N is INTEGER 00039 *> The order of the matrix A. N >= 0. 00040 *> \endverbatim 00041 *> 00042 *> \param[in] K 00043 *> \verbatim 00044 *> K is INTEGER 00045 *> The number of nonzero subdiagonals within the band of A. 00046 *> 0 <= K <= N-1. 00047 *> \endverbatim 00048 *> 00049 *> \param[in] D 00050 *> \verbatim 00051 *> D is REAL array, dimension (N) 00052 *> The diagonal elements of the diagonal matrix D. 00053 *> \endverbatim 00054 *> 00055 *> \param[out] A 00056 *> \verbatim 00057 *> A is REAL array, dimension (LDA,N) 00058 *> The generated n by n symmetric matrix A (the full matrix is 00059 *> stored). 00060 *> \endverbatim 00061 *> 00062 *> \param[in] LDA 00063 *> \verbatim 00064 *> LDA is INTEGER 00065 *> The leading dimension of the array A. LDA >= N. 00066 *> \endverbatim 00067 *> 00068 *> \param[in,out] ISEED 00069 *> \verbatim 00070 *> ISEED is INTEGER array, dimension (4) 00071 *> On entry, the seed of the random number generator; the array 00072 *> elements must be between 0 and 4095, and ISEED(4) must be 00073 *> odd. 00074 *> On exit, the seed is updated. 00075 *> \endverbatim 00076 *> 00077 *> \param[out] WORK 00078 *> \verbatim 00079 *> WORK is REAL array, dimension (2*N) 00080 *> \endverbatim 00081 *> 00082 *> \param[out] INFO 00083 *> \verbatim 00084 *> INFO is INTEGER 00085 *> = 0: successful exit 00086 *> < 0: if INFO = -i, the i-th argument had an illegal value 00087 *> \endverbatim 00088 * 00089 * Authors: 00090 * ======== 00091 * 00092 *> \author Univ. of Tennessee 00093 *> \author Univ. of California Berkeley 00094 *> \author Univ. of Colorado Denver 00095 *> \author NAG Ltd. 00096 * 00097 *> \date November 2011 00098 * 00099 *> \ingroup real_matgen 00100 * 00101 * ===================================================================== 00102 SUBROUTINE SLAGSY( N, K, D, A, LDA, ISEED, WORK, INFO ) 00103 * 00104 * -- LAPACK auxiliary routine (version 3.4.0) -- 00105 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00106 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00107 * November 2011 00108 * 00109 * .. Scalar Arguments .. 00110 INTEGER INFO, K, LDA, N 00111 * .. 00112 * .. Array Arguments .. 00113 INTEGER ISEED( 4 ) 00114 REAL A( LDA, * ), D( * ), WORK( * ) 00115 * .. 00116 * 00117 * ===================================================================== 00118 * 00119 * .. Parameters .. 00120 REAL ZERO, ONE, HALF 00121 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, HALF = 0.5E+0 ) 00122 * .. 00123 * .. Local Scalars .. 00124 INTEGER I, J 00125 REAL ALPHA, TAU, WA, WB, WN 00126 * .. 00127 * .. External Subroutines .. 00128 EXTERNAL SAXPY, SGEMV, SGER, SLARNV, SSCAL, SSYMV, 00129 $ SSYR2, XERBLA 00130 * .. 00131 * .. External Functions .. 00132 REAL SDOT, SNRM2 00133 EXTERNAL SDOT, SNRM2 00134 * .. 00135 * .. Intrinsic Functions .. 00136 INTRINSIC MAX, SIGN 00137 * .. 00138 * .. Executable Statements .. 00139 * 00140 * Test the input arguments 00141 * 00142 INFO = 0 00143 IF( N.LT.0 ) THEN 00144 INFO = -1 00145 ELSE IF( K.LT.0 .OR. K.GT.N-1 ) THEN 00146 INFO = -2 00147 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00148 INFO = -5 00149 END IF 00150 IF( INFO.LT.0 ) THEN 00151 CALL XERBLA( 'SLAGSY', -INFO ) 00152 RETURN 00153 END IF 00154 * 00155 * initialize lower triangle of A to diagonal matrix 00156 * 00157 DO 20 J = 1, N 00158 DO 10 I = J + 1, N 00159 A( I, J ) = ZERO 00160 10 CONTINUE 00161 20 CONTINUE 00162 DO 30 I = 1, N 00163 A( I, I ) = D( I ) 00164 30 CONTINUE 00165 * 00166 * Generate lower triangle of symmetric matrix 00167 * 00168 DO 40 I = N - 1, 1, -1 00169 * 00170 * generate random reflection 00171 * 00172 CALL SLARNV( 3, ISEED, N-I+1, WORK ) 00173 WN = SNRM2( N-I+1, WORK, 1 ) 00174 WA = SIGN( WN, WORK( 1 ) ) 00175 IF( WN.EQ.ZERO ) THEN 00176 TAU = ZERO 00177 ELSE 00178 WB = WORK( 1 ) + WA 00179 CALL SSCAL( N-I, ONE / WB, WORK( 2 ), 1 ) 00180 WORK( 1 ) = ONE 00181 TAU = WB / WA 00182 END IF 00183 * 00184 * apply random reflection to A(i:n,i:n) from the left 00185 * and the right 00186 * 00187 * compute y := tau * A * u 00188 * 00189 CALL SSYMV( 'Lower', N-I+1, TAU, A( I, I ), LDA, WORK, 1, ZERO, 00190 $ WORK( N+1 ), 1 ) 00191 * 00192 * compute v := y - 1/2 * tau * ( y, u ) * u 00193 * 00194 ALPHA = -HALF*TAU*SDOT( N-I+1, WORK( N+1 ), 1, WORK, 1 ) 00195 CALL SAXPY( N-I+1, ALPHA, WORK, 1, WORK( N+1 ), 1 ) 00196 * 00197 * apply the transformation as a rank-2 update to A(i:n,i:n) 00198 * 00199 CALL SSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1, 00200 $ A( I, I ), LDA ) 00201 40 CONTINUE 00202 * 00203 * Reduce number of subdiagonals to K 00204 * 00205 DO 60 I = 1, N - 1 - K 00206 * 00207 * generate reflection to annihilate A(k+i+1:n,i) 00208 * 00209 WN = SNRM2( N-K-I+1, A( K+I, I ), 1 ) 00210 WA = SIGN( WN, A( K+I, I ) ) 00211 IF( WN.EQ.ZERO ) THEN 00212 TAU = ZERO 00213 ELSE 00214 WB = A( K+I, I ) + WA 00215 CALL SSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 ) 00216 A( K+I, I ) = ONE 00217 TAU = WB / WA 00218 END IF 00219 * 00220 * apply reflection to A(k+i:n,i+1:k+i-1) from the left 00221 * 00222 CALL SGEMV( 'Transpose', N-K-I+1, K-1, ONE, A( K+I, I+1 ), LDA, 00223 $ A( K+I, I ), 1, ZERO, WORK, 1 ) 00224 CALL SGER( N-K-I+1, K-1, -TAU, A( K+I, I ), 1, WORK, 1, 00225 $ A( K+I, I+1 ), LDA ) 00226 * 00227 * apply reflection to A(k+i:n,k+i:n) from the left and the right 00228 * 00229 * compute y := tau * A * u 00230 * 00231 CALL SSYMV( 'Lower', N-K-I+1, TAU, A( K+I, K+I ), LDA, 00232 $ A( K+I, I ), 1, ZERO, WORK, 1 ) 00233 * 00234 * compute v := y - 1/2 * tau * ( y, u ) * u 00235 * 00236 ALPHA = -HALF*TAU*SDOT( N-K-I+1, WORK, 1, A( K+I, I ), 1 ) 00237 CALL SAXPY( N-K-I+1, ALPHA, A( K+I, I ), 1, WORK, 1 ) 00238 * 00239 * apply symmetric rank-2 update to A(k+i:n,k+i:n) 00240 * 00241 CALL SSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1, 00242 $ A( K+I, K+I ), LDA ) 00243 * 00244 A( K+I, I ) = -WA 00245 DO 50 J = K + I + 1, N 00246 A( J, I ) = ZERO 00247 50 CONTINUE 00248 60 CONTINUE 00249 * 00250 * Store full symmetric matrix 00251 * 00252 DO 80 J = 1, N 00253 DO 70 I = J + 1, N 00254 A( J, I ) = A( I, J ) 00255 70 CONTINUE 00256 80 CONTINUE 00257 RETURN 00258 * 00259 * End of SLAGSY 00260 * 00261 END