LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlaror.f
Go to the documentation of this file.
00001 *> \brief \b ZLAROR
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 ZLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
00012 * 
00013 *       .. Scalar Arguments ..
00014 *       CHARACTER          INIT, SIDE
00015 *       INTEGER            INFO, LDA, M, N
00016 *       ..
00017 *       .. Array Arguments ..
00018 *       INTEGER            ISEED( 4 )
00019 *       COMPLEX*16         A( LDA, * ), X( * )
00020 *       ..
00021 *  
00022 *
00023 *> \par Purpose:
00024 *  =============
00025 *>
00026 *> \verbatim
00027 *>
00028 *>    ZLAROR pre- or post-multiplies an M by N matrix A by a random
00029 *>    unitary matrix U, overwriting A. A may optionally be
00030 *>    initialized to the identity matrix before multiplying by U.
00031 *>    U is generated using the method of G.W. Stewart
00032 *>    ( SIAM J. Numer. Anal. 17, 1980, pp. 403-409 ).
00033 *>    (BLAS-2 version)
00034 *> \endverbatim
00035 *
00036 *  Arguments:
00037 *  ==========
00038 *
00039 *> \param[in] SIDE
00040 *> \verbatim
00041 *>          SIDE is CHARACTER*1
00042 *>           SIDE specifies whether A is multiplied on the left or right
00043 *>           by U.
00044 *>       SIDE = 'L'   Multiply A on the left (premultiply) by U
00045 *>       SIDE = 'R'   Multiply A on the right (postmultiply) by UC>       SIDE = 'C'   Multiply A on the left by U and the right by UC>       SIDE = 'T'   Multiply A on the left by U and the right by U'
00046 *>           Not modified.
00047 *> \endverbatim
00048 *>
00049 *> \param[in] INIT
00050 *> \verbatim
00051 *>          INIT is CHARACTER*1
00052 *>           INIT specifies whether or not A should be initialized to
00053 *>           the identity matrix.
00054 *>              INIT = 'I'   Initialize A to (a section of) the
00055 *>                           identity matrix before applying U.
00056 *>              INIT = 'N'   No initialization.  Apply U to the
00057 *>                           input matrix A.
00058 *>
00059 *>           INIT = 'I' may be used to generate square (i.e., unitary)
00060 *>           or rectangular orthogonal matrices (orthogonality being
00061 *>           in the sense of ZDOTC):
00062 *>
00063 *>           For square matrices, M=N, and SIDE many be either 'L' or
00064 *>           'R'; the rows will be orthogonal to each other, as will the
00065 *>           columns.
00066 *>           For rectangular matrices where M < N, SIDE = 'R' will
00067 *>           produce a dense matrix whose rows will be orthogonal and
00068 *>           whose columns will not, while SIDE = 'L' will produce a
00069 *>           matrix whose rows will be orthogonal, and whose first M
00070 *>           columns will be orthogonal, the remaining columns being
00071 *>           zero.
00072 *>           For matrices where M > N, just use the previous
00073 *>           explanation, interchanging 'L' and 'R' and "rows" and
00074 *>           "columns".
00075 *>
00076 *>           Not modified.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] M
00080 *> \verbatim
00081 *>          M is INTEGER
00082 *>           Number of rows of A. Not modified.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] N
00086 *> \verbatim
00087 *>          N is INTEGER
00088 *>           Number of columns of A. Not modified.
00089 *> \endverbatim
00090 *>
00091 *> \param[in,out] A
00092 *> \verbatim
00093 *>           A is COMPLEX*16 array, dimension ( LDA, N )
00094 *>           Input and output array. Overwritten by U A ( if SIDE = 'L' )
00095 *>           or by A U ( if SIDE = 'R' )
00096 *>           or by U A U* ( if SIDE = 'C')
00097 *>           or by U A U' ( if SIDE = 'T') on exit.
00098 *> \endverbatim
00099 *>
00100 *> \param[in] LDA
00101 *> \verbatim
00102 *>          LDA is INTEGER
00103 *>           Leading dimension of A. Must be at least MAX ( 1, M ).
00104 *>           Not modified.
00105 *> \endverbatim
00106 *>
00107 *> \param[in,out] ISEED
00108 *> \verbatim
00109 *>          ISEED is INTEGER array, dimension ( 4 )
00110 *>           On entry ISEED specifies the seed of the random number
00111 *>           generator. The array elements should be between 0 and 4095;
00112 *>           if not they will be reduced mod 4096.  Also, ISEED(4) must
00113 *>           be odd.  The random number generator uses a linear
00114 *>           congruential sequence limited to small integers, and so
00115 *>           should produce machine independent random numbers. The
00116 *>           values of ISEED are changed on exit, and can be used in the
00117 *>           next call to ZLAROR to continue the same random number
00118 *>           sequence.
00119 *>           Modified.
00120 *> \endverbatim
00121 *>
00122 *> \param[out] X
00123 *> \verbatim
00124 *>          X is COMPLEX*16 array, dimension ( 3*MAX( M, N ) )
00125 *>           Workspace. Of length:
00126 *>               2*M + N if SIDE = 'L',
00127 *>               2*N + M if SIDE = 'R',
00128 *>               3*N     if SIDE = 'C' or 'T'.
00129 *>           Modified.
00130 *> \endverbatim
00131 *>
00132 *> \param[out] INFO
00133 *> \verbatim
00134 *>          INFO is INTEGER
00135 *>           An error flag.  It is set to:
00136 *>            0  if no error.
00137 *>            1  if ZLARND returned a bad random number (installation
00138 *>               problem)
00139 *>           -1  if SIDE is not L, R, C, or T.
00140 *>           -3  if M is negative.
00141 *>           -4  if N is negative or if SIDE is C or T and N is not equal
00142 *>               to M.
00143 *>           -6  if LDA is less than M.
00144 *> \endverbatim
00145 *
00146 *  Authors:
00147 *  ========
00148 *
00149 *> \author Univ. of Tennessee 
00150 *> \author Univ. of California Berkeley 
00151 *> \author Univ. of Colorado Denver 
00152 *> \author NAG Ltd. 
00153 *
00154 *> \date November 2011
00155 *
00156 *> \ingroup complex16_matgen
00157 *
00158 *  =====================================================================
00159       SUBROUTINE ZLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
00160 *
00161 *  -- LAPACK auxiliary routine (version 3.4.0) --
00162 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00163 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00164 *     November 2011
00165 *
00166 *     .. Scalar Arguments ..
00167       CHARACTER          INIT, SIDE
00168       INTEGER            INFO, LDA, M, N
00169 *     ..
00170 *     .. Array Arguments ..
00171       INTEGER            ISEED( 4 )
00172       COMPLEX*16         A( LDA, * ), X( * )
00173 *     ..
00174 *
00175 *  =====================================================================
00176 *
00177 *     .. Parameters ..
00178       DOUBLE PRECISION   ZERO, ONE, TOOSML
00179       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
00180      $                   TOOSML = 1.0D-20 )
00181       COMPLEX*16         CZERO, CONE
00182       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00183      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00184 *     ..
00185 *     .. Local Scalars ..
00186       INTEGER            IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
00187       DOUBLE PRECISION   FACTOR, XABS, XNORM
00188       COMPLEX*16         CSIGN, XNORMS
00189 *     ..
00190 *     .. External Functions ..
00191       LOGICAL            LSAME
00192       DOUBLE PRECISION   DZNRM2
00193       COMPLEX*16         ZLARND
00194       EXTERNAL           LSAME, DZNRM2, ZLARND
00195 *     ..
00196 *     .. External Subroutines ..
00197       EXTERNAL           XERBLA, ZGEMV, ZGERC, ZLACGV, ZLASET, ZSCAL
00198 *     ..
00199 *     .. Intrinsic Functions ..
00200       INTRINSIC          ABS, DCMPLX, DCONJG
00201 *     ..
00202 *     .. Executable Statements ..
00203 *
00204       INFO = 0
00205       IF( N.EQ.0 .OR. M.EQ.0 )
00206      $   RETURN
00207 *
00208       ITYPE = 0
00209       IF( LSAME( SIDE, 'L' ) ) THEN
00210          ITYPE = 1
00211       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
00212          ITYPE = 2
00213       ELSE IF( LSAME( SIDE, 'C' ) ) THEN
00214          ITYPE = 3
00215       ELSE IF( LSAME( SIDE, 'T' ) ) THEN
00216          ITYPE = 4
00217       END IF
00218 *
00219 *     Check for argument errors.
00220 *
00221       IF( ITYPE.EQ.0 ) THEN
00222          INFO = -1
00223       ELSE IF( M.LT.0 ) THEN
00224          INFO = -3
00225       ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
00226          INFO = -4
00227       ELSE IF( LDA.LT.M ) THEN
00228          INFO = -6
00229       END IF
00230       IF( INFO.NE.0 ) THEN
00231          CALL XERBLA( 'ZLAROR', -INFO )
00232          RETURN
00233       END IF
00234 *
00235       IF( ITYPE.EQ.1 ) THEN
00236          NXFRM = M
00237       ELSE
00238          NXFRM = N
00239       END IF
00240 *
00241 *     Initialize A to the identity matrix if desired
00242 *
00243       IF( LSAME( INIT, 'I' ) )
00244      $   CALL ZLASET( 'Full', M, N, CZERO, CONE, A, LDA )
00245 *
00246 *     If no rotation possible, still multiply by
00247 *     a random complex number from the circle |x| = 1
00248 *
00249 *      2)      Compute Rotation by computing Householder
00250 *              Transformations H(2), H(3), ..., H(n).  Note that the
00251 *              order in which they are computed is irrelevant.
00252 *
00253       DO 10 J = 1, NXFRM
00254          X( J ) = CZERO
00255    10 CONTINUE
00256 *
00257       DO 30 IXFRM = 2, NXFRM
00258          KBEG = NXFRM - IXFRM + 1
00259 *
00260 *        Generate independent normal( 0, 1 ) random numbers
00261 *
00262          DO 20 J = KBEG, NXFRM
00263             X( J ) = ZLARND( 3, ISEED )
00264    20    CONTINUE
00265 *
00266 *        Generate a Householder transformation from the random vector X
00267 *
00268          XNORM = DZNRM2( IXFRM, X( KBEG ), 1 )
00269          XABS = ABS( X( KBEG ) )
00270          IF( XABS.NE.CZERO ) THEN
00271             CSIGN = X( KBEG ) / XABS
00272          ELSE
00273             CSIGN = CONE
00274          END IF
00275          XNORMS = CSIGN*XNORM
00276          X( NXFRM+KBEG ) = -CSIGN
00277          FACTOR = XNORM*( XNORM+XABS )
00278          IF( ABS( FACTOR ).LT.TOOSML ) THEN
00279             INFO = 1
00280             CALL XERBLA( 'ZLAROR', -INFO )
00281             RETURN
00282          ELSE
00283             FACTOR = ONE / FACTOR
00284          END IF
00285          X( KBEG ) = X( KBEG ) + XNORMS
00286 *
00287 *        Apply Householder transformation to A
00288 *
00289          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN
00290 *
00291 *           Apply H(k) on the left of A
00292 *
00293             CALL ZGEMV( 'C', IXFRM, N, CONE, A( KBEG, 1 ), LDA,
00294      $                  X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 )
00295             CALL ZGERC( IXFRM, N, -DCMPLX( FACTOR ), X( KBEG ), 1,
00296      $                  X( 2*NXFRM+1 ), 1, A( KBEG, 1 ), LDA )
00297 *
00298          END IF
00299 *
00300          IF( ITYPE.GE.2 .AND. ITYPE.LE.4 ) THEN
00301 *
00302 *           Apply H(k)* (or H(k)') on the right of A
00303 *
00304             IF( ITYPE.EQ.4 ) THEN
00305                CALL ZLACGV( IXFRM, X( KBEG ), 1 )
00306             END IF
00307 *
00308             CALL ZGEMV( 'N', M, IXFRM, CONE, A( 1, KBEG ), LDA,
00309      $                  X( KBEG ), 1, CZERO, X( 2*NXFRM+1 ), 1 )
00310             CALL ZGERC( M, IXFRM, -DCMPLX( FACTOR ), X( 2*NXFRM+1 ), 1,
00311      $                  X( KBEG ), 1, A( 1, KBEG ), LDA )
00312 *
00313          END IF
00314    30 CONTINUE
00315 *
00316       X( 1 ) = ZLARND( 3, ISEED )
00317       XABS = ABS( X( 1 ) )
00318       IF( XABS.NE.ZERO ) THEN
00319          CSIGN = X( 1 ) / XABS
00320       ELSE
00321          CSIGN = CONE
00322       END IF
00323       X( 2*NXFRM ) = CSIGN
00324 *
00325 *     Scale the matrix A by D.
00326 *
00327       IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 .OR. ITYPE.EQ.4 ) THEN
00328          DO 40 IROW = 1, M
00329             CALL ZSCAL( N, DCONJG( X( NXFRM+IROW ) ), A( IROW, 1 ),
00330      $                  LDA )
00331    40    CONTINUE
00332       END IF
00333 *
00334       IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
00335          DO 50 JCOL = 1, N
00336             CALL ZSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
00337    50    CONTINUE
00338       END IF
00339 *
00340       IF( ITYPE.EQ.4 ) THEN
00341          DO 60 JCOL = 1, N
00342             CALL ZSCAL( M, DCONJG( X( NXFRM+JCOL ) ), A( 1, JCOL ), 1 )
00343    60    CONTINUE
00344       END IF
00345       RETURN
00346 *
00347 *     End of ZLAROR
00348 *
00349       END
 All Files Functions