LAPACK 3.3.0

dlagsy.f

Go to the documentation of this file.
00001       SUBROUTINE DLAGSY( 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   A( LDA, * ), D( * ), WORK( * )
00013 *     ..
00014 *
00015 *  Purpose
00016 *  =======
00017 *
00018 *  DLAGSY generates a real symmetric matrix A, by pre- and post-
00019 *  multiplying a real diagonal matrix D with a random orthogonal matrix:
00020 *  A = U*D*U'. The semi-bandwidth may then be reduced to k by additional
00021 *  orthogonal transformations.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  N       (input) INTEGER
00027 *          The order of the matrix A.  N >= 0.
00028 *
00029 *  K       (input) INTEGER
00030 *          The number of nonzero subdiagonals within the band of A.
00031 *          0 <= K <= N-1.
00032 *
00033 *  D       (input) DOUBLE PRECISION array, dimension (N)
00034 *          The diagonal elements of the diagonal matrix D.
00035 *
00036 *  A       (output) DOUBLE PRECISION array, dimension (LDA,N)
00037 *          The generated n by n symmetric matrix A (the full matrix is
00038 *          stored).
00039 *
00040 *  LDA     (input) INTEGER
00041 *          The leading dimension of the array A.  LDA >= N.
00042 *
00043 *  ISEED   (input/output) INTEGER array, dimension (4)
00044 *          On entry, the seed of the random number generator; the array
00045 *          elements must be between 0 and 4095, and ISEED(4) must be
00046 *          odd.
00047 *          On exit, the seed is updated.
00048 *
00049 *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
00050 *
00051 *  INFO    (output) INTEGER
00052 *          = 0: successful exit
00053 *          < 0: if INFO = -i, the i-th argument had an illegal value
00054 *
00055 *  =====================================================================
00056 *
00057 *     .. Parameters ..
00058       DOUBLE PRECISION   ZERO, ONE, HALF
00059       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 )
00060 *     ..
00061 *     .. Local Scalars ..
00062       INTEGER            I, J
00063       DOUBLE PRECISION   ALPHA, TAU, WA, WB, WN
00064 *     ..
00065 *     .. External Subroutines ..
00066       EXTERNAL           DAXPY, DGEMV, DGER, DLARNV, DSCAL, DSYMV,
00067      $                   DSYR2, XERBLA
00068 *     ..
00069 *     .. External Functions ..
00070       DOUBLE PRECISION   DDOT, DNRM2
00071       EXTERNAL           DDOT, DNRM2
00072 *     ..
00073 *     .. Intrinsic Functions ..
00074       INTRINSIC          MAX, SIGN
00075 *     ..
00076 *     .. Executable Statements ..
00077 *
00078 *     Test the input arguments
00079 *
00080       INFO = 0
00081       IF( N.LT.0 ) THEN
00082          INFO = -1
00083       ELSE IF( K.LT.0 .OR. K.GT.N-1 ) THEN
00084          INFO = -2
00085       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00086          INFO = -5
00087       END IF
00088       IF( INFO.LT.0 ) THEN
00089          CALL XERBLA( 'DLAGSY', -INFO )
00090          RETURN
00091       END IF
00092 *
00093 *     initialize lower triangle of A to diagonal matrix
00094 *
00095       DO 20 J = 1, N
00096          DO 10 I = J + 1, N
00097             A( I, J ) = ZERO
00098    10    CONTINUE
00099    20 CONTINUE
00100       DO 30 I = 1, N
00101          A( I, I ) = D( I )
00102    30 CONTINUE
00103 *
00104 *     Generate lower triangle of symmetric matrix
00105 *
00106       DO 40 I = N - 1, 1, -1
00107 *
00108 *        generate random reflection
00109 *
00110          CALL DLARNV( 3, ISEED, N-I+1, WORK )
00111          WN = DNRM2( N-I+1, WORK, 1 )
00112          WA = SIGN( WN, WORK( 1 ) )
00113          IF( WN.EQ.ZERO ) THEN
00114             TAU = ZERO
00115          ELSE
00116             WB = WORK( 1 ) + WA
00117             CALL DSCAL( N-I, ONE / WB, WORK( 2 ), 1 )
00118             WORK( 1 ) = ONE
00119             TAU = WB / WA
00120          END IF
00121 *
00122 *        apply random reflection to A(i:n,i:n) from the left
00123 *        and the right
00124 *
00125 *        compute  y := tau * A * u
00126 *
00127          CALL DSYMV( 'Lower', N-I+1, TAU, A( I, I ), LDA, WORK, 1, ZERO,
00128      $               WORK( N+1 ), 1 )
00129 *
00130 *        compute  v := y - 1/2 * tau * ( y, u ) * u
00131 *
00132          ALPHA = -HALF*TAU*DDOT( N-I+1, WORK( N+1 ), 1, WORK, 1 )
00133          CALL DAXPY( N-I+1, ALPHA, WORK, 1, WORK( N+1 ), 1 )
00134 *
00135 *        apply the transformation as a rank-2 update to A(i:n,i:n)
00136 *
00137          CALL DSYR2( 'Lower', N-I+1, -ONE, WORK, 1, WORK( N+1 ), 1,
00138      $               A( I, I ), LDA )
00139    40 CONTINUE
00140 *
00141 *     Reduce number of subdiagonals to K
00142 *
00143       DO 60 I = 1, N - 1 - K
00144 *
00145 *        generate reflection to annihilate A(k+i+1:n,i)
00146 *
00147          WN = DNRM2( N-K-I+1, A( K+I, I ), 1 )
00148          WA = SIGN( WN, A( K+I, I ) )
00149          IF( WN.EQ.ZERO ) THEN
00150             TAU = ZERO
00151          ELSE
00152             WB = A( K+I, I ) + WA
00153             CALL DSCAL( N-K-I, ONE / WB, A( K+I+1, I ), 1 )
00154             A( K+I, I ) = ONE
00155             TAU = WB / WA
00156          END IF
00157 *
00158 *        apply reflection to A(k+i:n,i+1:k+i-1) from the left
00159 *
00160          CALL DGEMV( 'Transpose', N-K-I+1, K-1, ONE, A( K+I, I+1 ), LDA,
00161      $               A( K+I, I ), 1, ZERO, WORK, 1 )
00162          CALL DGER( N-K-I+1, K-1, -TAU, A( K+I, I ), 1, WORK, 1,
00163      $              A( K+I, I+1 ), LDA )
00164 *
00165 *        apply reflection to A(k+i:n,k+i:n) from the left and the right
00166 *
00167 *        compute  y := tau * A * u
00168 *
00169          CALL DSYMV( 'Lower', N-K-I+1, TAU, A( K+I, K+I ), LDA,
00170      $               A( K+I, I ), 1, ZERO, WORK, 1 )
00171 *
00172 *        compute  v := y - 1/2 * tau * ( y, u ) * u
00173 *
00174          ALPHA = -HALF*TAU*DDOT( N-K-I+1, WORK, 1, A( K+I, I ), 1 )
00175          CALL DAXPY( N-K-I+1, ALPHA, A( K+I, I ), 1, WORK, 1 )
00176 *
00177 *        apply symmetric rank-2 update to A(k+i:n,k+i:n)
00178 *
00179          CALL DSYR2( 'Lower', N-K-I+1, -ONE, A( K+I, I ), 1, WORK, 1,
00180      $               A( K+I, K+I ), LDA )
00181 *
00182          A( K+I, I ) = -WA
00183          DO 50 J = K + I + 1, N
00184             A( J, I ) = ZERO
00185    50    CONTINUE
00186    60 CONTINUE
00187 *
00188 *     Store full symmetric matrix
00189 *
00190       DO 80 J = 1, N
00191          DO 70 I = J + 1, N
00192             A( J, I ) = A( I, J )
00193    70    CONTINUE
00194    80 CONTINUE
00195       RETURN
00196 *
00197 *     End of DLAGSY
00198 *
00199       END
 All Files Functions