LAPACK 3.3.1
Linear Algebra PACKage

dlaror.f

Go to the documentation of this file.
00001       SUBROUTINE DLAROR( 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 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       CHARACTER          INIT, SIDE
00009       INTEGER            INFO, LDA, M, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            ISEED( 4 )
00013       DOUBLE PRECISION   A( LDA, * ), X( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DLAROR pre- or post-multiplies an M by N matrix A by a random
00020 *  orthogonal matrix U, overwriting A.  A may optionally be initialized
00021 *  to the identity matrix before multiplying by U.  U is generated using
00022 *  the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  SIDE    (input) CHARACTER*1
00028 *          Specifies whether A is multiplied on the left or right by U.
00029 *          = 'L':         Multiply A on the left (premultiply) by U
00030 *          = 'R':         Multiply A on the right (postmultiply) by U'
00031 *          = 'C' or 'T':  Multiply A on the left by U and the right
00032 *                          by U' (Here, U' means U-transpose.)
00033 *
00034 *  INIT    (input) CHARACTER*1
00035 *          Specifies whether or not A should be initialized to the
00036 *          identity matrix.
00037 *          = 'I':  Initialize A to (a section of) the identity matrix
00038 *                   before applying U.
00039 *          = 'N':  No initialization.  Apply U to the input matrix A.
00040 *
00041 *          INIT = 'I' may be used to generate square or rectangular
00042 *          orthogonal matrices:
00043 *
00044 *          For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
00045 *          to each other, as will the columns.
00046 *
00047 *          If M < N, SIDE = 'R' produces a dense matrix whose rows are
00048 *          orthogonal and whose columns are not, while SIDE = 'L'
00049 *          produces a matrix whose rows are orthogonal, and whose first
00050 *          M columns are orthogonal, and whose remaining columns are
00051 *          zero.
00052 *
00053 *          If M > N, SIDE = 'L' produces a dense matrix whose columns
00054 *          are orthogonal and whose rows are not, while SIDE = 'R'
00055 *          produces a matrix whose columns are orthogonal, and whose
00056 *          first M rows are orthogonal, and whose remaining rows are
00057 *          zero.
00058 *
00059 *  M       (input) INTEGER
00060 *          The number of rows of A.
00061 *
00062 *  N       (input) INTEGER
00063 *          The number of columns of A.
00064 *
00065 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
00066 *          On entry, the array A.
00067 *          On exit, overwritten by U A ( if SIDE = 'L' ),
00068 *           or by A U ( if SIDE = 'R' ),
00069 *           or by U A U' ( if SIDE = 'C' or 'T').
00070 *
00071 *  LDA     (input) INTEGER
00072 *          The leading dimension of the array A.  LDA >= max(1,M).
00073 *
00074 *  ISEED   (input/output) INTEGER array, dimension (4)
00075 *          On entry ISEED specifies the seed of the random number
00076 *          generator. The array elements should be between 0 and 4095;
00077 *          if not they will be reduced mod 4096.  Also, ISEED(4) must
00078 *          be odd.  The random number generator uses a linear
00079 *          congruential sequence limited to small integers, and so
00080 *          should produce machine independent random numbers. The
00081 *          values of ISEED are changed on exit, and can be used in the
00082 *          next call to DLAROR to continue the same random number
00083 *          sequence.
00084 *
00085 *  X       (workspace) DOUBLE PRECISION array, dimension (3*MAX( M, N ))
00086 *          Workspace of length
00087 *              2*M + N if SIDE = 'L',
00088 *              2*N + M if SIDE = 'R',
00089 *              3*N     if SIDE = 'C' or 'T'.
00090 *
00091 *  INFO    (output) INTEGER
00092 *          An error flag.  It is set to:
00093 *          = 0:  normal return
00094 *          < 0:  if INFO = -k, the k-th argument had an illegal value
00095 *          = 1:  if the random numbers generated by DLARND are bad.
00096 *
00097 *  =====================================================================
00098 *
00099 *     .. Parameters ..
00100       DOUBLE PRECISION   ZERO, ONE, TOOSML
00101       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
00102      $                   TOOSML = 1.0D-20 )
00103 *     ..
00104 *     .. Local Scalars ..
00105       INTEGER            IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
00106       DOUBLE PRECISION   FACTOR, XNORM, XNORMS
00107 *     ..
00108 *     .. External Functions ..
00109       LOGICAL            LSAME
00110       DOUBLE PRECISION   DLARND, DNRM2
00111       EXTERNAL           LSAME, DLARND, DNRM2
00112 *     ..
00113 *     .. External Subroutines ..
00114       EXTERNAL           DGEMV, DGER, DLASET, DSCAL, XERBLA
00115 *     ..
00116 *     .. Intrinsic Functions ..
00117       INTRINSIC          ABS, SIGN
00118 *     ..
00119 *     .. Executable Statements ..
00120 *
00121       IF( N.EQ.0 .OR. M.EQ.0 )
00122      $   RETURN
00123 *
00124       ITYPE = 0
00125       IF( LSAME( SIDE, 'L' ) ) THEN
00126          ITYPE = 1
00127       ELSE IF( LSAME( SIDE, 'R' ) ) THEN
00128          ITYPE = 2
00129       ELSE IF( LSAME( SIDE, 'C' ) .OR. LSAME( SIDE, 'T' ) ) THEN
00130          ITYPE = 3
00131       END IF
00132 *
00133 *     Check for argument errors.
00134 *
00135       INFO = 0
00136       IF( ITYPE.EQ.0 ) THEN
00137          INFO = -1
00138       ELSE IF( M.LT.0 ) THEN
00139          INFO = -3
00140       ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
00141          INFO = -4
00142       ELSE IF( LDA.LT.M ) THEN
00143          INFO = -6
00144       END IF
00145       IF( INFO.NE.0 ) THEN
00146          CALL XERBLA( 'DLAROR', -INFO )
00147          RETURN
00148       END IF
00149 *
00150       IF( ITYPE.EQ.1 ) THEN
00151          NXFRM = M
00152       ELSE
00153          NXFRM = N
00154       END IF
00155 *
00156 *     Initialize A to the identity matrix if desired
00157 *
00158       IF( LSAME( INIT, 'I' ) )
00159      $   CALL DLASET( 'Full', M, N, ZERO, ONE, A, LDA )
00160 *
00161 *     If no rotation possible, multiply by random +/-1
00162 *
00163 *     Compute rotation by computing Householder transformations
00164 *     H(2), H(3), ..., H(nhouse)
00165 *
00166       DO 10 J = 1, NXFRM
00167          X( J ) = ZERO
00168    10 CONTINUE
00169 *
00170       DO 30 IXFRM = 2, NXFRM
00171          KBEG = NXFRM - IXFRM + 1
00172 *
00173 *        Generate independent normal( 0, 1 ) random numbers
00174 *
00175          DO 20 J = KBEG, NXFRM
00176             X( J ) = DLARND( 3, ISEED )
00177    20    CONTINUE
00178 *
00179 *        Generate a Householder transformation from the random vector X
00180 *
00181          XNORM = DNRM2( IXFRM, X( KBEG ), 1 )
00182          XNORMS = SIGN( XNORM, X( KBEG ) )
00183          X( KBEG+NXFRM ) = SIGN( ONE, -X( KBEG ) )
00184          FACTOR = XNORMS*( XNORMS+X( KBEG ) )
00185          IF( ABS( FACTOR ).LT.TOOSML ) THEN
00186             INFO = 1
00187             CALL XERBLA( 'DLAROR', INFO )
00188             RETURN
00189          ELSE
00190             FACTOR = ONE / FACTOR
00191          END IF
00192          X( KBEG ) = X( KBEG ) + XNORMS
00193 *
00194 *        Apply Householder transformation to A
00195 *
00196          IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
00197 *
00198 *           Apply H(k) from the left.
00199 *
00200             CALL DGEMV( 'T', IXFRM, N, ONE, A( KBEG, 1 ), LDA,
00201      $                  X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
00202             CALL DGER( IXFRM, N, -FACTOR, X( KBEG ), 1, X( 2*NXFRM+1 ),
00203      $                 1, A( KBEG, 1 ), LDA )
00204 *
00205          END IF
00206 *
00207          IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
00208 *
00209 *           Apply H(k) from the right.
00210 *
00211             CALL DGEMV( 'N', M, IXFRM, ONE, A( 1, KBEG ), LDA,
00212      $                  X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
00213             CALL DGER( M, IXFRM, -FACTOR, X( 2*NXFRM+1 ), 1, X( KBEG ),
00214      $                 1, A( 1, KBEG ), LDA )
00215 *
00216          END IF
00217    30 CONTINUE
00218 *
00219       X( 2*NXFRM ) = SIGN( ONE, DLARND( 3, ISEED ) )
00220 *
00221 *     Scale the matrix A by D.
00222 *
00223       IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
00224          DO 40 IROW = 1, M
00225             CALL DSCAL( N, X( NXFRM+IROW ), A( IROW, 1 ), LDA )
00226    40    CONTINUE
00227       END IF
00228 *
00229       IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
00230          DO 50 JCOL = 1, N
00231             CALL DSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
00232    50    CONTINUE
00233       END IF
00234       RETURN
00235 *
00236 *     End of DLAROR
00237 *
00238       END
 All Files Functions