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