LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlagsy.f
Go to the documentation of this file.
00001 *> \brief \b DLAGSY
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 DLAGSY( 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 *       DOUBLE PRECISION   A( LDA, * ), D( * ), WORK( * )
00019 *       ..
00020 *  
00021 *
00022 *> \par Purpose:
00023 *  =============
00024 *>
00025 *> \verbatim
00026 *>
00027 *> DLAGSY 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 DOUBLE PRECISION array, dimension (N)
00052 *>          The diagonal elements of the diagonal matrix D.
00053 *> \endverbatim
00054 *>
00055 *> \param[out] A
00056 *> \verbatim
00057 *>          A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 double_matgen
00100 *
00101 *  =====================================================================
00102       SUBROUTINE DLAGSY( 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       DOUBLE PRECISION   A( LDA, * ), D( * ), WORK( * )
00115 *     ..
00116 *
00117 *  =====================================================================
00118 *
00119 *     .. Parameters ..
00120       DOUBLE PRECISION   ZERO, ONE, HALF
00121       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 )
00122 *     ..
00123 *     .. Local Scalars ..
00124       INTEGER            I, J
00125       DOUBLE PRECISION   ALPHA, TAU, WA, WB, WN
00126 *     ..
00127 *     .. External Subroutines ..
00128       EXTERNAL           DAXPY, DGEMV, DGER, DLARNV, DSCAL, DSYMV,
00129      $                   DSYR2, XERBLA
00130 *     ..
00131 *     .. External Functions ..
00132       DOUBLE PRECISION   DDOT, DNRM2
00133       EXTERNAL           DDOT, DNRM2
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( 'DLAGSY', -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 DLARNV( 3, ISEED, N-I+1, WORK )
00173          WN = DNRM2( 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 DSCAL( 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 DSYMV( '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*DDOT( N-I+1, WORK( N+1 ), 1, WORK, 1 )
00195          CALL DAXPY( 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 DSYR2( '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 = DNRM2( 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 DSCAL( 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 DGEMV( '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 DGER( 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 DSYMV( '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*DDOT( N-K-I+1, WORK, 1, A( K+I, I ), 1 )
00237          CALL DAXPY( 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 DSYR2( '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 DLAGSY
00260 *
00261       END
 All Files Functions