LAPACK: TESTING/MATGEN/zlaror.f Source File

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