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