![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DORMRZ( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, 2: $ WORK, LWORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * January 2007 8: * 9: * .. Scalar Arguments .. 10: CHARACTER SIDE, TRANS 11: INTEGER INFO, K, L, LDA, LDC, LWORK, M, N 12: * .. 13: * .. Array Arguments .. 14: DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DORMRZ overwrites the general real M-by-N matrix C with 21: * 22: * SIDE = 'L' SIDE = 'R' 23: * TRANS = 'N': Q * C C * Q 24: * TRANS = 'T': Q**T * C C * Q**T 25: * 26: * where Q is a real orthogonal matrix defined as the product of k 27: * elementary reflectors 28: * 29: * Q = H(1) H(2) . . . H(k) 30: * 31: * as returned by DTZRZF. Q is of order M if SIDE = 'L' and of order N 32: * if SIDE = 'R'. 33: * 34: * Arguments 35: * ========= 36: * 37: * SIDE (input) CHARACTER*1 38: * = 'L': apply Q or Q**T from the Left; 39: * = 'R': apply Q or Q**T from the Right. 40: * 41: * TRANS (input) CHARACTER*1 42: * = 'N': No transpose, apply Q; 43: * = 'T': Transpose, apply Q**T. 44: * 45: * M (input) INTEGER 46: * The number of rows of the matrix C. M >= 0. 47: * 48: * N (input) INTEGER 49: * The number of columns of the matrix C. N >= 0. 50: * 51: * K (input) INTEGER 52: * The number of elementary reflectors whose product defines 53: * the matrix Q. 54: * If SIDE = 'L', M >= K >= 0; 55: * if SIDE = 'R', N >= K >= 0. 56: * 57: * L (input) INTEGER 58: * The number of columns of the matrix A containing 59: * the meaningful part of the Householder reflectors. 60: * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0. 61: * 62: * A (input) DOUBLE PRECISION array, dimension 63: * (LDA,M) if SIDE = 'L', 64: * (LDA,N) if SIDE = 'R' 65: * The i-th row must contain the vector which defines the 66: * elementary reflector H(i), for i = 1,2,...,k, as returned by 67: * DTZRZF in the last k rows of its array argument A. 68: * A is modified by the routine but restored on exit. 69: * 70: * LDA (input) INTEGER 71: * The leading dimension of the array A. LDA >= max(1,K). 72: * 73: * TAU (input) DOUBLE PRECISION array, dimension (K) 74: * TAU(i) must contain the scalar factor of the elementary 75: * reflector H(i), as returned by DTZRZF. 76: * 77: * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 78: * On entry, the M-by-N matrix C. 79: * On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q. 80: * 81: * LDC (input) INTEGER 82: * The leading dimension of the array C. LDC >= max(1,M). 83: * 84: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 85: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 86: * 87: * LWORK (input) INTEGER 88: * The dimension of the array WORK. 89: * If SIDE = 'L', LWORK >= max(1,N); 90: * if SIDE = 'R', LWORK >= max(1,M). 91: * For optimum performance LWORK >= N*NB if SIDE = 'L', and 92: * LWORK >= M*NB if SIDE = 'R', where NB is the optimal 93: * blocksize. 94: * 95: * If LWORK = -1, then a workspace query is assumed; the routine 96: * only calculates the optimal size of the WORK array, returns 97: * this value as the first entry of the WORK array, and no error 98: * message related to LWORK is issued by XERBLA. 99: * 100: * INFO (output) INTEGER 101: * = 0: successful exit 102: * < 0: if INFO = -i, the i-th argument had an illegal value 103: * 104: * Further Details 105: * =============== 106: * 107: * Based on contributions by 108: * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA 109: * 110: * ===================================================================== 111: * 112: * .. Parameters .. 113: INTEGER NBMAX, LDT 114: PARAMETER ( NBMAX = 64, LDT = NBMAX+1 ) 115: * .. 116: * .. Local Scalars .. 117: LOGICAL LEFT, LQUERY, NOTRAN 118: CHARACTER TRANST 119: INTEGER I, I1, I2, I3, IB, IC, IINFO, IWS, JA, JC, 120: $ LDWORK, LWKOPT, MI, NB, NBMIN, NI, NQ, NW 121: * .. 122: * .. Local Arrays .. 123: DOUBLE PRECISION T( LDT, NBMAX ) 124: * .. 125: * .. External Functions .. 126: LOGICAL LSAME 127: INTEGER ILAENV 128: EXTERNAL LSAME, ILAENV 129: * .. 130: * .. External Subroutines .. 131: EXTERNAL DLARZB, DLARZT, DORMR3, XERBLA 132: * .. 133: * .. Intrinsic Functions .. 134: INTRINSIC MAX, MIN 135: * .. 136: * .. Executable Statements .. 137: * 138: * Test the input arguments 139: * 140: INFO = 0 141: LEFT = LSAME( SIDE, 'L' ) 142: NOTRAN = LSAME( TRANS, 'N' ) 143: LQUERY = ( LWORK.EQ.-1 ) 144: * 145: * NQ is the order of Q and NW is the minimum dimension of WORK 146: * 147: IF( LEFT ) THEN 148: NQ = M 149: NW = MAX( 1, N ) 150: ELSE 151: NQ = N 152: NW = MAX( 1, M ) 153: END IF 154: IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN 155: INFO = -1 156: ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 157: INFO = -2 158: ELSE IF( M.LT.0 ) THEN 159: INFO = -3 160: ELSE IF( N.LT.0 ) THEN 161: INFO = -4 162: ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN 163: INFO = -5 164: ELSE IF( L.LT.0 .OR. ( LEFT .AND. ( L.GT.M ) ) .OR. 165: $ ( .NOT.LEFT .AND. ( L.GT.N ) ) ) THEN 166: INFO = -6 167: ELSE IF( LDA.LT.MAX( 1, K ) ) THEN 168: INFO = -8 169: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 170: INFO = -11 171: END IF 172: * 173: IF( INFO.EQ.0 ) THEN 174: IF( M.EQ.0 .OR. N.EQ.0 ) THEN 175: LWKOPT = 1 176: ELSE 177: * 178: * Determine the block size. NB may be at most NBMAX, where 179: * NBMAX is used to define the local array T. 180: * 181: NB = MIN( NBMAX, ILAENV( 1, 'DORMRQ', SIDE // TRANS, M, N, 182: $ K, -1 ) ) 183: LWKOPT = NW*NB 184: END IF 185: WORK( 1 ) = LWKOPT 186: * 187: IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN 188: INFO = -13 189: END IF 190: END IF 191: * 192: IF( INFO.NE.0 ) THEN 193: CALL XERBLA( 'DORMRZ', -INFO ) 194: RETURN 195: ELSE IF( LQUERY ) THEN 196: RETURN 197: END IF 198: * 199: * Quick return if possible 200: * 201: IF( M.EQ.0 .OR. N.EQ.0 ) THEN 202: WORK( 1 ) = 1 203: RETURN 204: END IF 205: * 206: NBMIN = 2 207: LDWORK = NW 208: IF( NB.GT.1 .AND. NB.LT.K ) THEN 209: IWS = NW*NB 210: IF( LWORK.LT.IWS ) THEN 211: NB = LWORK / LDWORK 212: NBMIN = MAX( 2, ILAENV( 2, 'DORMRQ', SIDE // TRANS, M, N, K, 213: $ -1 ) ) 214: END IF 215: ELSE 216: IWS = NW 217: END IF 218: * 219: IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN 220: * 221: * Use unblocked code 222: * 223: CALL DORMR3( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, 224: $ WORK, IINFO ) 225: ELSE 226: * 227: * Use blocked code 228: * 229: IF( ( LEFT .AND. .NOT.NOTRAN ) .OR. 230: $ ( .NOT.LEFT .AND. NOTRAN ) ) THEN 231: I1 = 1 232: I2 = K 233: I3 = NB 234: ELSE 235: I1 = ( ( K-1 ) / NB )*NB + 1 236: I2 = 1 237: I3 = -NB 238: END IF 239: * 240: IF( LEFT ) THEN 241: NI = N 242: JC = 1 243: JA = M - L + 1 244: ELSE 245: MI = M 246: IC = 1 247: JA = N - L + 1 248: END IF 249: * 250: IF( NOTRAN ) THEN 251: TRANST = 'T' 252: ELSE 253: TRANST = 'N' 254: END IF 255: * 256: DO 10 I = I1, I2, I3 257: IB = MIN( NB, K-I+1 ) 258: * 259: * Form the triangular factor of the block reflector 260: * H = H(i+ib-1) . . . H(i+1) H(i) 261: * 262: CALL DLARZT( 'Backward', 'Rowwise', L, IB, A( I, JA ), LDA, 263: $ TAU( I ), T, LDT ) 264: * 265: IF( LEFT ) THEN 266: * 267: * H or H' is applied to C(i:m,1:n) 268: * 269: MI = M - I + 1 270: IC = I 271: ELSE 272: * 273: * H or H' is applied to C(1:m,i:n) 274: * 275: NI = N - I + 1 276: JC = I 277: END IF 278: * 279: * Apply H or H' 280: * 281: CALL DLARZB( SIDE, TRANST, 'Backward', 'Rowwise', MI, NI, 282: $ IB, L, A( I, JA ), LDA, T, LDT, C( IC, JC ), 283: $ LDC, WORK, LDWORK ) 284: 10 CONTINUE 285: * 286: END IF 287: * 288: WORK( 1 ) = LWKOPT 289: * 290: RETURN 291: * 292: * End of DORMRZ 293: * 294: END