LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clagge.f
Go to the documentation of this file.
00001 *> \brief \b CLAGGE
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 CLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       INTEGER            INFO, KL, KU, LDA, M, N
00015 *       ..
00016 *       .. Array Arguments ..
00017 *       INTEGER            ISEED( 4 )
00018 *       REAL               D( * )
00019 *       COMPLEX            A( LDA, * ), WORK( * )
00020 *       ..
00021 *  
00022 *
00023 *> \par Purpose:
00024 *  =============
00025 *>
00026 *> \verbatim
00027 *>
00028 *> CLAGGE generates a complex general m by n matrix A, by pre- and post-
00029 *> multiplying a real diagonal matrix D with random unitary matrices:
00030 *> A = U*D*V. The lower and upper bandwidths may then be reduced to
00031 *> kl and ku by additional unitary transformations.
00032 *> \endverbatim
00033 *
00034 *  Arguments:
00035 *  ==========
00036 *
00037 *> \param[in] M
00038 *> \verbatim
00039 *>          M is INTEGER
00040 *>          The number of rows of the matrix A.  M >= 0.
00041 *> \endverbatim
00042 *>
00043 *> \param[in] N
00044 *> \verbatim
00045 *>          N is INTEGER
00046 *>          The number of columns of the matrix A.  N >= 0.
00047 *> \endverbatim
00048 *>
00049 *> \param[in] KL
00050 *> \verbatim
00051 *>          KL is INTEGER
00052 *>          The number of nonzero subdiagonals within the band of A.
00053 *>          0 <= KL <= M-1.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] KU
00057 *> \verbatim
00058 *>          KU is INTEGER
00059 *>          The number of nonzero superdiagonals within the band of A.
00060 *>          0 <= KU <= N-1.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] D
00064 *> \verbatim
00065 *>          D is REAL array, dimension (min(M,N))
00066 *>          The diagonal elements of the diagonal matrix D.
00067 *> \endverbatim
00068 *>
00069 *> \param[out] A
00070 *> \verbatim
00071 *>          A is COMPLEX array, dimension (LDA,N)
00072 *>          The generated m by n matrix A.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] LDA
00076 *> \verbatim
00077 *>          LDA is INTEGER
00078 *>          The leading dimension of the array A.  LDA >= M.
00079 *> \endverbatim
00080 *>
00081 *> \param[in,out] ISEED
00082 *> \verbatim
00083 *>          ISEED is INTEGER array, dimension (4)
00084 *>          On entry, the seed of the random number generator; the array
00085 *>          elements must be between 0 and 4095, and ISEED(4) must be
00086 *>          odd.
00087 *>          On exit, the seed is updated.
00088 *> \endverbatim
00089 *>
00090 *> \param[out] WORK
00091 *> \verbatim
00092 *>          WORK is COMPLEX array, dimension (M+N)
00093 *> \endverbatim
00094 *>
00095 *> \param[out] INFO
00096 *> \verbatim
00097 *>          INFO is INTEGER
00098 *>          = 0: successful exit
00099 *>          < 0: if INFO = -i, the i-th argument had an illegal value
00100 *> \endverbatim
00101 *
00102 *  Authors:
00103 *  ========
00104 *
00105 *> \author Univ. of Tennessee 
00106 *> \author Univ. of California Berkeley 
00107 *> \author Univ. of Colorado Denver 
00108 *> \author NAG Ltd. 
00109 *
00110 *> \date November 2011
00111 *
00112 *> \ingroup complex_matgen
00113 *
00114 *  =====================================================================
00115       SUBROUTINE CLAGGE( M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO )
00116 *
00117 *  -- LAPACK auxiliary routine (version 3.4.0) --
00118 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00119 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00120 *     November 2011
00121 *
00122 *     .. Scalar Arguments ..
00123       INTEGER            INFO, KL, KU, LDA, M, N
00124 *     ..
00125 *     .. Array Arguments ..
00126       INTEGER            ISEED( 4 )
00127       REAL               D( * )
00128       COMPLEX            A( LDA, * ), WORK( * )
00129 *     ..
00130 *
00131 *  =====================================================================
00132 *
00133 *     .. Parameters ..
00134       COMPLEX            ZERO, ONE
00135       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ),
00136      $                   ONE = ( 1.0E+0, 0.0E+0 ) )
00137 *     ..
00138 *     .. Local Scalars ..
00139       INTEGER            I, J
00140       REAL               WN
00141       COMPLEX            TAU, WA, WB
00142 *     ..
00143 *     .. External Subroutines ..
00144       EXTERNAL           CGEMV, CGERC, CLACGV, CLARNV, CSCAL, XERBLA
00145 *     ..
00146 *     .. Intrinsic Functions ..
00147       INTRINSIC          ABS, MAX, MIN, REAL
00148 *     ..
00149 *     .. External Functions ..
00150       REAL               SCNRM2
00151       EXTERNAL           SCNRM2
00152 *     ..
00153 *     .. Executable Statements ..
00154 *
00155 *     Test the input arguments
00156 *
00157       INFO = 0
00158       IF( M.LT.0 ) THEN
00159          INFO = -1
00160       ELSE IF( N.LT.0 ) THEN
00161          INFO = -2
00162       ELSE IF( KL.LT.0 .OR. KL.GT.M-1 ) THEN
00163          INFO = -3
00164       ELSE IF( KU.LT.0 .OR. KU.GT.N-1 ) THEN
00165          INFO = -4
00166       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00167          INFO = -7
00168       END IF
00169       IF( INFO.LT.0 ) THEN
00170          CALL XERBLA( 'CLAGGE', -INFO )
00171          RETURN
00172       END IF
00173 *
00174 *     initialize A to diagonal matrix
00175 *
00176       DO 20 J = 1, N
00177          DO 10 I = 1, M
00178             A( I, J ) = ZERO
00179    10    CONTINUE
00180    20 CONTINUE
00181       DO 30 I = 1, MIN( M, N )
00182          A( I, I ) = D( I )
00183    30 CONTINUE
00184 *
00185 *     pre- and post-multiply A by random unitary matrices
00186 *
00187       DO 40 I = MIN( M, N ), 1, -1
00188          IF( I.LT.M ) THEN
00189 *
00190 *           generate random reflection
00191 *
00192             CALL CLARNV( 3, ISEED, M-I+1, WORK )
00193             WN = SCNRM2( M-I+1, WORK, 1 )
00194             WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 )
00195             IF( WN.EQ.ZERO ) THEN
00196                TAU = ZERO
00197             ELSE
00198                WB = WORK( 1 ) + WA
00199                CALL CSCAL( M-I, ONE / WB, WORK( 2 ), 1 )
00200                WORK( 1 ) = ONE
00201                TAU = REAL( WB / WA )
00202             END IF
00203 *
00204 *           multiply A(i:m,i:n) by random reflection from the left
00205 *
00206             CALL CGEMV( 'Conjugate transpose', M-I+1, N-I+1, ONE,
00207      $                  A( I, I ), LDA, WORK, 1, ZERO, WORK( M+1 ), 1 )
00208             CALL CGERC( M-I+1, N-I+1, -TAU, WORK, 1, WORK( M+1 ), 1,
00209      $                  A( I, I ), LDA )
00210          END IF
00211          IF( I.LT.N ) THEN
00212 *
00213 *           generate random reflection
00214 *
00215             CALL CLARNV( 3, ISEED, N-I+1, WORK )
00216             WN = SCNRM2( N-I+1, WORK, 1 )
00217             WA = ( WN / ABS( WORK( 1 ) ) )*WORK( 1 )
00218             IF( WN.EQ.ZERO ) THEN
00219                TAU = ZERO
00220             ELSE
00221                WB = WORK( 1 ) + WA
00222                CALL CSCAL( N-I, ONE / WB, WORK( 2 ), 1 )
00223                WORK( 1 ) = ONE
00224                TAU = REAL( WB / WA )
00225             END IF
00226 *
00227 *           multiply A(i:m,i:n) by random reflection from the right
00228 *
00229             CALL CGEMV( 'No transpose', M-I+1, N-I+1, ONE, A( I, I ),
00230      $                  LDA, WORK, 1, ZERO, WORK( N+1 ), 1 )
00231             CALL CGERC( M-I+1, N-I+1, -TAU, WORK( N+1 ), 1, WORK, 1,
00232      $                  A( I, I ), LDA )
00233          END IF
00234    40 CONTINUE
00235 *
00236 *     Reduce number of subdiagonals to KL and number of superdiagonals
00237 *     to KU
00238 *
00239       DO 70 I = 1, MAX( M-1-KL, N-1-KU )
00240          IF( KL.LE.KU ) THEN
00241 *
00242 *           annihilate subdiagonal elements first (necessary if KL = 0)
00243 *
00244             IF( I.LE.MIN( M-1-KL, N ) ) THEN
00245 *
00246 *              generate reflection to annihilate A(kl+i+1:m,i)
00247 *
00248                WN = SCNRM2( M-KL-I+1, A( KL+I, I ), 1 )
00249                WA = ( WN / ABS( A( KL+I, I ) ) )*A( KL+I, I )
00250                IF( WN.EQ.ZERO ) THEN
00251                   TAU = ZERO
00252                ELSE
00253                   WB = A( KL+I, I ) + WA
00254                   CALL CSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 )
00255                   A( KL+I, I ) = ONE
00256                   TAU = REAL( WB / WA )
00257                END IF
00258 *
00259 *              apply reflection to A(kl+i:m,i+1:n) from the left
00260 *
00261                CALL CGEMV( 'Conjugate transpose', M-KL-I+1, N-I, ONE,
00262      $                     A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO,
00263      $                     WORK, 1 )
00264                CALL CGERC( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK,
00265      $                     1, A( KL+I, I+1 ), LDA )
00266                A( KL+I, I ) = -WA
00267             END IF
00268 *
00269             IF( I.LE.MIN( N-1-KU, M ) ) THEN
00270 *
00271 *              generate reflection to annihilate A(i,ku+i+1:n)
00272 *
00273                WN = SCNRM2( N-KU-I+1, A( I, KU+I ), LDA )
00274                WA = ( WN / ABS( A( I, KU+I ) ) )*A( I, KU+I )
00275                IF( WN.EQ.ZERO ) THEN
00276                   TAU = ZERO
00277                ELSE
00278                   WB = A( I, KU+I ) + WA
00279                   CALL CSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA )
00280                   A( I, KU+I ) = ONE
00281                   TAU = REAL( WB / WA )
00282                END IF
00283 *
00284 *              apply reflection to A(i+1:m,ku+i:n) from the right
00285 *
00286                CALL CLACGV( N-KU-I+1, A( I, KU+I ), LDA )
00287                CALL CGEMV( 'No transpose', M-I, N-KU-I+1, ONE,
00288      $                     A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO,
00289      $                     WORK, 1 )
00290                CALL CGERC( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ),
00291      $                     LDA, A( I+1, KU+I ), LDA )
00292                A( I, KU+I ) = -WA
00293             END IF
00294          ELSE
00295 *
00296 *           annihilate superdiagonal elements first (necessary if
00297 *           KU = 0)
00298 *
00299             IF( I.LE.MIN( N-1-KU, M ) ) THEN
00300 *
00301 *              generate reflection to annihilate A(i,ku+i+1:n)
00302 *
00303                WN = SCNRM2( N-KU-I+1, A( I, KU+I ), LDA )
00304                WA = ( WN / ABS( A( I, KU+I ) ) )*A( I, KU+I )
00305                IF( WN.EQ.ZERO ) THEN
00306                   TAU = ZERO
00307                ELSE
00308                   WB = A( I, KU+I ) + WA
00309                   CALL CSCAL( N-KU-I, ONE / WB, A( I, KU+I+1 ), LDA )
00310                   A( I, KU+I ) = ONE
00311                   TAU = REAL( WB / WA )
00312                END IF
00313 *
00314 *              apply reflection to A(i+1:m,ku+i:n) from the right
00315 *
00316                CALL CLACGV( N-KU-I+1, A( I, KU+I ), LDA )
00317                CALL CGEMV( 'No transpose', M-I, N-KU-I+1, ONE,
00318      $                     A( I+1, KU+I ), LDA, A( I, KU+I ), LDA, ZERO,
00319      $                     WORK, 1 )
00320                CALL CGERC( M-I, N-KU-I+1, -TAU, WORK, 1, A( I, KU+I ),
00321      $                     LDA, A( I+1, KU+I ), LDA )
00322                A( I, KU+I ) = -WA
00323             END IF
00324 *
00325             IF( I.LE.MIN( M-1-KL, N ) ) THEN
00326 *
00327 *              generate reflection to annihilate A(kl+i+1:m,i)
00328 *
00329                WN = SCNRM2( M-KL-I+1, A( KL+I, I ), 1 )
00330                WA = ( WN / ABS( A( KL+I, I ) ) )*A( KL+I, I )
00331                IF( WN.EQ.ZERO ) THEN
00332                   TAU = ZERO
00333                ELSE
00334                   WB = A( KL+I, I ) + WA
00335                   CALL CSCAL( M-KL-I, ONE / WB, A( KL+I+1, I ), 1 )
00336                   A( KL+I, I ) = ONE
00337                   TAU = REAL( WB / WA )
00338                END IF
00339 *
00340 *              apply reflection to A(kl+i:m,i+1:n) from the left
00341 *
00342                CALL CGEMV( 'Conjugate transpose', M-KL-I+1, N-I, ONE,
00343      $                     A( KL+I, I+1 ), LDA, A( KL+I, I ), 1, ZERO,
00344      $                     WORK, 1 )
00345                CALL CGERC( M-KL-I+1, N-I, -TAU, A( KL+I, I ), 1, WORK,
00346      $                     1, A( KL+I, I+1 ), LDA )
00347                A( KL+I, I ) = -WA
00348             END IF
00349          END IF
00350 *
00351          DO 50 J = KL + I + 1, M
00352             A( J, I ) = ZERO
00353    50    CONTINUE
00354 *
00355          DO 60 J = KU + I + 1, N
00356             A( I, J ) = ZERO
00357    60    CONTINUE
00358    70 CONTINUE
00359       RETURN
00360 *
00361 *     End of CLAGGE
00362 *
00363       END
 All Files Functions