You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
302 lines
8.8 KiB
302 lines
8.8 KiB
2 years ago
|
*> \brief \b DLAROR
|
||
|
|
*
|
||
|
|
* =========== DOCUMENTATION ===========
|
||
|
|
*
|
||
|
|
* Online html documentation available at
|
||
|
|
* http://www.netlib.org/lapack/explore-html/
|
||
|
|
*
|
||
|
|
* Definition:
|
||
|
|
* ===========
|
||
|
|
*
|
||
|
|
* SUBROUTINE DLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
|
||
|
|
*
|
||
|
|
* .. Scalar Arguments ..
|
||
|
|
* CHARACTER INIT, SIDE
|
||
|
|
* INTEGER INFO, LDA, M, N
|
||
|
|
* ..
|
||
|
|
* .. Array Arguments ..
|
||
|
|
* INTEGER ISEED( 4 )
|
||
|
|
* DOUBLE PRECISION A( LDA, * ), X( * )
|
||
|
|
* ..
|
||
|
|
*
|
||
|
|
*
|
||
|
|
*> \par Purpose:
|
||
|
|
* =============
|
||
|
|
*>
|
||
|
|
*> \verbatim
|
||
|
|
*>
|
||
|
|
*> DLAROR pre- or post-multiplies an M by N matrix A by a random
|
||
|
|
*> orthogonal matrix U, overwriting A. A may optionally be initialized
|
||
|
|
*> to the identity matrix before multiplying by U. U is generated using
|
||
|
|
*> the method of G.W. Stewart (SIAM J. Numer. Anal. 17, 1980, 403-409).
|
||
|
|
*> \endverbatim
|
||
|
|
*
|
||
|
|
* Arguments:
|
||
|
|
* ==========
|
||
|
|
*
|
||
|
|
*> \param[in] SIDE
|
||
|
|
*> \verbatim
|
||
|
|
*> SIDE is CHARACTER*1
|
||
|
|
*> Specifies whether A is multiplied on the left or right by U.
|
||
|
|
*> = 'L': Multiply A on the left (premultiply) by U
|
||
|
|
*> = 'R': Multiply A on the right (postmultiply) by U'
|
||
|
|
*> = 'C' or 'T': Multiply A on the left by U and the right
|
||
|
|
*> by U' (Here, U' means U-transpose.)
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[in] INIT
|
||
|
|
*> \verbatim
|
||
|
|
*> INIT is CHARACTER*1
|
||
|
|
*> Specifies whether or not A should be initialized to the
|
||
|
|
*> identity matrix.
|
||
|
|
*> = 'I': Initialize A to (a section of) the identity matrix
|
||
|
|
*> before applying U.
|
||
|
|
*> = 'N': No initialization. Apply U to the input matrix A.
|
||
|
|
*>
|
||
|
|
*> INIT = 'I' may be used to generate square or rectangular
|
||
|
|
*> orthogonal matrices:
|
||
|
|
*>
|
||
|
|
*> For M = N and SIDE = 'L' or 'R', the rows will be orthogonal
|
||
|
|
*> to each other, as will the columns.
|
||
|
|
*>
|
||
|
|
*> If M < N, SIDE = 'R' produces a dense matrix whose rows are
|
||
|
|
*> orthogonal and whose columns are not, while SIDE = 'L'
|
||
|
|
*> produces a matrix whose rows are orthogonal, and whose first
|
||
|
|
*> M columns are orthogonal, and whose remaining columns are
|
||
|
|
*> zero.
|
||
|
|
*>
|
||
|
|
*> If M > N, SIDE = 'L' produces a dense matrix whose columns
|
||
|
|
*> are orthogonal and whose rows are not, while SIDE = 'R'
|
||
|
|
*> produces a matrix whose columns are orthogonal, and whose
|
||
|
|
*> first M rows are orthogonal, and whose remaining rows are
|
||
|
|
*> zero.
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[in] M
|
||
|
|
*> \verbatim
|
||
|
|
*> M is INTEGER
|
||
|
|
*> The number of rows of A.
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[in] N
|
||
|
|
*> \verbatim
|
||
|
|
*> N is INTEGER
|
||
|
|
*> The number of columns of A.
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[in,out] A
|
||
|
|
*> \verbatim
|
||
|
|
*> A is DOUBLE PRECISION array, dimension (LDA, N)
|
||
|
|
*> On entry, the array A.
|
||
|
|
*> On exit, overwritten by U A ( if SIDE = 'L' ),
|
||
|
|
*> or by A U ( if SIDE = 'R' ),
|
||
|
|
*> or by U A U' ( if SIDE = 'C' or 'T').
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[in] LDA
|
||
|
|
*> \verbatim
|
||
|
|
*> LDA is INTEGER
|
||
|
|
*> The leading dimension of the array A. LDA >= max(1,M).
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[in,out] ISEED
|
||
|
|
*> \verbatim
|
||
|
|
*> ISEED is INTEGER array, dimension (4)
|
||
|
|
*> On entry ISEED specifies the seed of the random number
|
||
|
|
*> generator. The array elements should be between 0 and 4095;
|
||
|
|
*> if not they will be reduced mod 4096. Also, ISEED(4) must
|
||
|
|
*> be odd. The random number generator uses a linear
|
||
|
|
*> congruential sequence limited to small integers, and so
|
||
|
|
*> should produce machine independent random numbers. The
|
||
|
|
*> values of ISEED are changed on exit, and can be used in the
|
||
|
|
*> next call to DLAROR to continue the same random number
|
||
|
|
*> sequence.
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[out] X
|
||
|
|
*> \verbatim
|
||
|
|
*> X is DOUBLE PRECISION array, dimension (3*MAX( M, N ))
|
||
|
|
*> Workspace of length
|
||
|
|
*> 2*M + N if SIDE = 'L',
|
||
|
|
*> 2*N + M if SIDE = 'R',
|
||
|
|
*> 3*N if SIDE = 'C' or 'T'.
|
||
|
|
*> \endverbatim
|
||
|
|
*>
|
||
|
|
*> \param[out] INFO
|
||
|
|
*> \verbatim
|
||
|
|
*> INFO is INTEGER
|
||
|
|
*> An error flag. It is set to:
|
||
|
|
*> = 0: normal return
|
||
|
|
*> < 0: if INFO = -k, the k-th argument had an illegal value
|
||
|
|
*> = 1: if the random numbers generated by DLARND are bad.
|
||
|
|
*> \endverbatim
|
||
|
|
*
|
||
|
|
* Authors:
|
||
|
|
* ========
|
||
|
|
*
|
||
|
|
*> \author Univ. of Tennessee
|
||
|
|
*> \author Univ. of California Berkeley
|
||
|
|
*> \author Univ. of Colorado Denver
|
||
|
|
*> \author NAG Ltd.
|
||
|
|
*
|
||
|
|
*> \ingroup double_matgen
|
||
|
|
*
|
||
|
|
* =====================================================================
|
||
|
|
SUBROUTINE DLAROR( SIDE, INIT, M, N, A, LDA, ISEED, X, INFO )
|
||
|
|
*
|
||
|
|
* -- LAPACK auxiliary routine --
|
||
|
|
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
||
|
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
||
|
|
*
|
||
|
|
* .. Scalar Arguments ..
|
||
|
|
CHARACTER INIT, SIDE
|
||
|
|
INTEGER INFO, LDA, M, N
|
||
|
|
* ..
|
||
|
|
* .. Array Arguments ..
|
||
|
|
INTEGER ISEED( 4 )
|
||
|
|
DOUBLE PRECISION A( LDA, * ), X( * )
|
||
|
|
* ..
|
||
|
|
*
|
||
|
|
* =====================================================================
|
||
|
|
*
|
||
|
|
* .. Parameters ..
|
||
|
|
DOUBLE PRECISION ZERO, ONE, TOOSML
|
||
|
|
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0,
|
||
|
|
$ TOOSML = 1.0D-20 )
|
||
|
|
* ..
|
||
|
|
* .. Local Scalars ..
|
||
|
|
INTEGER IROW, ITYPE, IXFRM, J, JCOL, KBEG, NXFRM
|
||
|
|
DOUBLE PRECISION FACTOR, XNORM, XNORMS
|
||
|
|
* ..
|
||
|
|
* .. External Functions ..
|
||
|
|
LOGICAL LSAME
|
||
|
|
DOUBLE PRECISION DLARND, DNRM2
|
||
|
|
EXTERNAL LSAME, DLARND, DNRM2
|
||
|
|
* ..
|
||
|
|
* .. External Subroutines ..
|
||
|
|
EXTERNAL DGEMV, DGER, DLASET, DSCAL, XERBLA
|
||
|
|
* ..
|
||
|
|
* .. Intrinsic Functions ..
|
||
|
|
INTRINSIC ABS, SIGN
|
||
|
|
* ..
|
||
|
|
* .. Executable Statements ..
|
||
|
|
*
|
||
|
|
INFO = 0
|
||
|
|
IF( N.EQ.0 .OR. M.EQ.0 )
|
||
|
|
$ RETURN
|
||
|
|
*
|
||
|
|
ITYPE = 0
|
||
|
|
IF( LSAME( SIDE, 'L' ) ) THEN
|
||
|
|
ITYPE = 1
|
||
|
|
ELSE IF( LSAME( SIDE, 'R' ) ) THEN
|
||
|
|
ITYPE = 2
|
||
|
|
ELSE IF( LSAME( SIDE, 'C' ) .OR. LSAME( SIDE, 'T' ) ) THEN
|
||
|
|
ITYPE = 3
|
||
|
|
END IF
|
||
|
|
*
|
||
|
|
* Check for argument errors.
|
||
|
|
*
|
||
|
|
IF( ITYPE.EQ.0 ) THEN
|
||
|
|
INFO = -1
|
||
|
|
ELSE IF( M.LT.0 ) THEN
|
||
|
|
INFO = -3
|
||
|
|
ELSE IF( N.LT.0 .OR. ( ITYPE.EQ.3 .AND. N.NE.M ) ) THEN
|
||
|
|
INFO = -4
|
||
|
|
ELSE IF( LDA.LT.M ) THEN
|
||
|
|
INFO = -6
|
||
|
|
END IF
|
||
|
|
IF( INFO.NE.0 ) THEN
|
||
|
|
CALL XERBLA( 'DLAROR', -INFO )
|
||
|
|
RETURN
|
||
|
|
END IF
|
||
|
|
*
|
||
|
|
IF( ITYPE.EQ.1 ) THEN
|
||
|
|
NXFRM = M
|
||
|
|
ELSE
|
||
|
|
NXFRM = N
|
||
|
|
END IF
|
||
|
|
*
|
||
|
|
* Initialize A to the identity matrix if desired
|
||
|
|
*
|
||
|
|
IF( LSAME( INIT, 'I' ) )
|
||
|
|
$ CALL DLASET( 'Full', M, N, ZERO, ONE, A, LDA )
|
||
|
|
*
|
||
|
|
* If no rotation possible, multiply by random +/-1
|
||
|
|
*
|
||
|
|
* Compute rotation by computing Householder transformations
|
||
|
|
* H(2), H(3), ..., H(nhouse)
|
||
|
|
*
|
||
|
|
DO 10 J = 1, NXFRM
|
||
|
|
X( J ) = ZERO
|
||
|
|
10 CONTINUE
|
||
|
|
*
|
||
|
|
DO 30 IXFRM = 2, NXFRM
|
||
|
|
KBEG = NXFRM - IXFRM + 1
|
||
|
|
*
|
||
|
|
* Generate independent normal( 0, 1 ) random numbers
|
||
|
|
*
|
||
|
|
DO 20 J = KBEG, NXFRM
|
||
|
|
X( J ) = DLARND( 3, ISEED )
|
||
|
|
20 CONTINUE
|
||
|
|
*
|
||
|
|
* Generate a Householder transformation from the random vector X
|
||
|
|
*
|
||
|
|
XNORM = DNRM2( IXFRM, X( KBEG ), 1 )
|
||
|
|
XNORMS = SIGN( XNORM, X( KBEG ) )
|
||
|
|
X( KBEG+NXFRM ) = SIGN( ONE, -X( KBEG ) )
|
||
|
|
FACTOR = XNORMS*( XNORMS+X( KBEG ) )
|
||
|
|
IF( ABS( FACTOR ).LT.TOOSML ) THEN
|
||
|
|
INFO = 1
|
||
|
|
CALL XERBLA( 'DLAROR', INFO )
|
||
|
|
RETURN
|
||
|
|
ELSE
|
||
|
|
FACTOR = ONE / FACTOR
|
||
|
|
END IF
|
||
|
|
X( KBEG ) = X( KBEG ) + XNORMS
|
||
|
|
*
|
||
|
|
* Apply Householder transformation to A
|
||
|
|
*
|
||
|
|
IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
|
||
|
|
*
|
||
|
|
* Apply H(k) from the left.
|
||
|
|
*
|
||
|
|
CALL DGEMV( 'T', IXFRM, N, ONE, A( KBEG, 1 ), LDA,
|
||
|
|
$ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
|
||
|
|
CALL DGER( IXFRM, N, -FACTOR, X( KBEG ), 1, X( 2*NXFRM+1 ),
|
||
|
|
$ 1, A( KBEG, 1 ), LDA )
|
||
|
|
*
|
||
|
|
END IF
|
||
|
|
*
|
||
|
|
IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
|
||
|
|
*
|
||
|
|
* Apply H(k) from the right.
|
||
|
|
*
|
||
|
|
CALL DGEMV( 'N', M, IXFRM, ONE, A( 1, KBEG ), LDA,
|
||
|
|
$ X( KBEG ), 1, ZERO, X( 2*NXFRM+1 ), 1 )
|
||
|
|
CALL DGER( M, IXFRM, -FACTOR, X( 2*NXFRM+1 ), 1, X( KBEG ),
|
||
|
|
$ 1, A( 1, KBEG ), LDA )
|
||
|
|
*
|
||
|
|
END IF
|
||
|
|
30 CONTINUE
|
||
|
|
*
|
||
|
|
X( 2*NXFRM ) = SIGN( ONE, DLARND( 3, ISEED ) )
|
||
|
|
*
|
||
|
|
* Scale the matrix A by D.
|
||
|
|
*
|
||
|
|
IF( ITYPE.EQ.1 .OR. ITYPE.EQ.3 ) THEN
|
||
|
|
DO 40 IROW = 1, M
|
||
|
|
CALL DSCAL( N, X( NXFRM+IROW ), A( IROW, 1 ), LDA )
|
||
|
|
40 CONTINUE
|
||
|
|
END IF
|
||
|
|
*
|
||
|
|
IF( ITYPE.EQ.2 .OR. ITYPE.EQ.3 ) THEN
|
||
|
|
DO 50 JCOL = 1, N
|
||
|
|
CALL DSCAL( M, X( NXFRM+JCOL ), A( 1, JCOL ), 1 )
|
||
|
|
50 CONTINUE
|
||
|
|
END IF
|
||
|
|
RETURN
|
||
|
|
*
|
||
|
|
* End of DLAROR
|
||
|
|
*
|
||
|
|
END
|