![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DORGRQ( M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) 2: * 3: * -- LAPACK routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * .. Scalar Arguments .. 9: INTEGER INFO, K, LDA, LWORK, M, N 10: * .. 11: * .. Array Arguments .. 12: DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) 13: * .. 14: * 15: * Purpose 16: * ======= 17: * 18: * DORGRQ generates an M-by-N real matrix Q with orthonormal rows, 19: * which is defined as the last M rows of a product of K elementary 20: * reflectors of order N 21: * 22: * Q = H(1) H(2) . . . H(k) 23: * 24: * as returned by DGERQF. 25: * 26: * Arguments 27: * ========= 28: * 29: * M (input) INTEGER 30: * The number of rows of the matrix Q. M >= 0. 31: * 32: * N (input) INTEGER 33: * The number of columns of the matrix Q. N >= M. 34: * 35: * K (input) INTEGER 36: * The number of elementary reflectors whose product defines the 37: * matrix Q. M >= K >= 0. 38: * 39: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 40: * On entry, the (m-k+i)-th row must contain the vector which 41: * defines the elementary reflector H(i), for i = 1,2,...,k, as 42: * returned by DGERQF in the last k rows of its array argument 43: * A. 44: * On exit, the M-by-N matrix Q. 45: * 46: * LDA (input) INTEGER 47: * The first dimension of the array A. LDA >= max(1,M). 48: * 49: * TAU (input) DOUBLE PRECISION array, dimension (K) 50: * TAU(i) must contain the scalar factor of the elementary 51: * reflector H(i), as returned by DGERQF. 52: * 53: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 54: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 55: * 56: * LWORK (input) INTEGER 57: * The dimension of the array WORK. LWORK >= max(1,M). 58: * For optimum performance LWORK >= M*NB, where NB is the 59: * optimal blocksize. 60: * 61: * If LWORK = -1, then a workspace query is assumed; the routine 62: * only calculates the optimal size of the WORK array, returns 63: * this value as the first entry of the WORK array, and no error 64: * message related to LWORK is issued by XERBLA. 65: * 66: * INFO (output) INTEGER 67: * = 0: successful exit 68: * < 0: if INFO = -i, the i-th argument has an illegal value 69: * 70: * ===================================================================== 71: * 72: * .. Parameters .. 73: DOUBLE PRECISION ZERO 74: PARAMETER ( ZERO = 0.0D+0 ) 75: * .. 76: * .. Local Scalars .. 77: LOGICAL LQUERY 78: INTEGER I, IB, II, IINFO, IWS, J, KK, L, LDWORK, 79: $ LWKOPT, NB, NBMIN, NX 80: * .. 81: * .. External Subroutines .. 82: EXTERNAL DLARFB, DLARFT, DORGR2, XERBLA 83: * .. 84: * .. Intrinsic Functions .. 85: INTRINSIC MAX, MIN 86: * .. 87: * .. External Functions .. 88: INTEGER ILAENV 89: EXTERNAL ILAENV 90: * .. 91: * .. Executable Statements .. 92: * 93: * Test the input arguments 94: * 95: INFO = 0 96: LQUERY = ( LWORK.EQ.-1 ) 97: IF( M.LT.0 ) THEN 98: INFO = -1 99: ELSE IF( N.LT.M ) THEN 100: INFO = -2 101: ELSE IF( K.LT.0 .OR. K.GT.M ) THEN 102: INFO = -3 103: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 104: INFO = -5 105: END IF 106: * 107: IF( INFO.EQ.0 ) THEN 108: IF( M.LE.0 ) THEN 109: LWKOPT = 1 110: ELSE 111: NB = ILAENV( 1, 'DORGRQ', ' ', M, N, K, -1 ) 112: LWKOPT = M*NB 113: END IF 114: WORK( 1 ) = LWKOPT 115: * 116: IF( LWORK.LT.MAX( 1, M ) .AND. .NOT.LQUERY ) THEN 117: INFO = -8 118: END IF 119: END IF 120: * 121: IF( INFO.NE.0 ) THEN 122: CALL XERBLA( 'DORGRQ', -INFO ) 123: RETURN 124: ELSE IF( LQUERY ) THEN 125: RETURN 126: END IF 127: * 128: * Quick return if possible 129: * 130: IF( M.LE.0 ) THEN 131: RETURN 132: END IF 133: * 134: NBMIN = 2 135: NX = 0 136: IWS = M 137: IF( NB.GT.1 .AND. NB.LT.K ) THEN 138: * 139: * Determine when to cross over from blocked to unblocked code. 140: * 141: NX = MAX( 0, ILAENV( 3, 'DORGRQ', ' ', M, N, K, -1 ) ) 142: IF( NX.LT.K ) THEN 143: * 144: * Determine if workspace is large enough for blocked code. 145: * 146: LDWORK = M 147: IWS = LDWORK*NB 148: IF( LWORK.LT.IWS ) THEN 149: * 150: * Not enough workspace to use optimal NB: reduce NB and 151: * determine the minimum value of NB. 152: * 153: NB = LWORK / LDWORK 154: NBMIN = MAX( 2, ILAENV( 2, 'DORGRQ', ' ', M, N, K, -1 ) ) 155: END IF 156: END IF 157: END IF 158: * 159: IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN 160: * 161: * Use blocked code after the first block. 162: * The last kk rows are handled by the block method. 163: * 164: KK = MIN( K, ( ( K-NX+NB-1 ) / NB )*NB ) 165: * 166: * Set A(1:m-kk,n-kk+1:n) to zero. 167: * 168: DO 20 J = N - KK + 1, N 169: DO 10 I = 1, M - KK 170: A( I, J ) = ZERO 171: 10 CONTINUE 172: 20 CONTINUE 173: ELSE 174: KK = 0 175: END IF 176: * 177: * Use unblocked code for the first or only block. 178: * 179: CALL DORGR2( M-KK, N-KK, K-KK, A, LDA, TAU, WORK, IINFO ) 180: * 181: IF( KK.GT.0 ) THEN 182: * 183: * Use blocked code 184: * 185: DO 50 I = K - KK + 1, K, NB 186: IB = MIN( NB, K-I+1 ) 187: II = M - K + I 188: IF( II.GT.1 ) THEN 189: * 190: * Form the triangular factor of the block reflector 191: * H = H(i+ib-1) . . . H(i+1) H(i) 192: * 193: CALL DLARFT( 'Backward', 'Rowwise', N-K+I+IB-1, IB, 194: $ A( II, 1 ), LDA, TAU( I ), WORK, LDWORK ) 195: * 196: * Apply H' to A(1:m-k+i-1,1:n-k+i+ib-1) from the right 197: * 198: CALL DLARFB( 'Right', 'Transpose', 'Backward', 'Rowwise', 199: $ II-1, N-K+I+IB-1, IB, A( II, 1 ), LDA, WORK, 200: $ LDWORK, A, LDA, WORK( IB+1 ), LDWORK ) 201: END IF 202: * 203: * Apply H' to columns 1:n-k+i+ib-1 of current block 204: * 205: CALL DORGR2( IB, N-K+I+IB-1, IB, A( II, 1 ), LDA, TAU( I ), 206: $ WORK, IINFO ) 207: * 208: * Set columns n-k+i+ib:n of current block to zero 209: * 210: DO 40 L = N - K + I + IB, N 211: DO 30 J = II, II + IB - 1 212: A( J, L ) = ZERO 213: 30 CONTINUE 214: 40 CONTINUE 215: 50 CONTINUE 216: END IF 217: * 218: WORK( 1 ) = IWS 219: RETURN 220: * 221: * End of DORGRQ 222: * 223: END