1*> \brief \b CLATSQR 2* 3* Definition: 4* =========== 5* 6* SUBROUTINE CLATSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, 7* LWORK, INFO) 8* 9* .. Scalar Arguments .. 10* INTEGER INFO, LDA, M, N, MB, NB, LDT, LWORK 11* .. 12* .. Array Arguments .. 13* COMPLEX A( LDA, * ), T( LDT, * ), WORK( * ) 14* .. 15* 16* 17*> \par Purpose: 18* ============= 19*> 20*> \verbatim 21*> 22*> CLATSQR computes a blocked Tall-Skinny QR factorization of 23*> a complex M-by-N matrix A for M >= N: 24*> 25*> A = Q * ( R ), 26*> ( 0 ) 27*> 28*> where: 29*> 30*> Q is a M-by-M orthogonal matrix, stored on exit in an implicit 31*> form in the elements below the diagonal of the array A and in 32*> the elements of the array T; 33*> 34*> R is an upper-triangular N-by-N matrix, stored on exit in 35*> the elements on and above the diagonal of the array A. 36*> 37*> 0 is a (M-N)-by-N zero matrix, and is not stored. 38*> 39*> \endverbatim 40* 41* Arguments: 42* ========== 43* 44*> \param[in] M 45*> \verbatim 46*> M is INTEGER 47*> The number of rows of the matrix A. M >= 0. 48*> \endverbatim 49*> 50*> \param[in] N 51*> \verbatim 52*> N is INTEGER 53*> The number of columns of the matrix A. M >= N >= 0. 54*> \endverbatim 55*> 56*> \param[in] MB 57*> \verbatim 58*> MB is INTEGER 59*> The row block size to be used in the blocked QR. 60*> MB > N. 61*> \endverbatim 62*> 63*> \param[in] NB 64*> \verbatim 65*> NB is INTEGER 66*> The column block size to be used in the blocked QR. 67*> N >= NB >= 1. 68*> \endverbatim 69*> 70*> \param[in,out] A 71*> \verbatim 72*> A is COMPLEX array, dimension (LDA,N) 73*> On entry, the M-by-N matrix A. 74*> On exit, the elements on and above the diagonal 75*> of the array contain the N-by-N upper triangular matrix R; 76*> the elements below the diagonal represent Q by the columns 77*> of blocked V (see Further Details). 78*> \endverbatim 79*> 80*> \param[in] LDA 81*> \verbatim 82*> LDA is INTEGER 83*> The leading dimension of the array A. LDA >= max(1,M). 84*> \endverbatim 85*> 86*> \param[out] T 87*> \verbatim 88*> T is COMPLEX array, 89*> dimension (LDT, N * Number_of_row_blocks) 90*> where Number_of_row_blocks = CEIL((M-N)/(MB-N)) 91*> The blocked upper triangular block reflectors stored in compact form 92*> as a sequence of upper triangular blocks. 93*> See Further Details below. 94*> \endverbatim 95*> 96*> \param[in] LDT 97*> \verbatim 98*> LDT is INTEGER 99*> The leading dimension of the array T. LDT >= NB. 100*> \endverbatim 101*> 102*> \param[out] WORK 103*> \verbatim 104*> (workspace) COMPLEX array, dimension (MAX(1,LWORK)) 105*> \endverbatim 106*> 107*> \param[in] LWORK 108*> \verbatim 109*> The dimension of the array WORK. LWORK >= NB*N. 110*> If LWORK = -1, then a workspace query is assumed; the routine 111*> only calculates the optimal size of the WORK array, returns 112*> this value as the first entry of the WORK array, and no error 113*> message related to LWORK is issued by XERBLA. 114*> \endverbatim 115*> 116*> \param[out] INFO 117*> \verbatim 118*> INFO is INTEGER 119*> = 0: successful exit 120*> < 0: if INFO = -i, the i-th argument had an illegal value 121*> \endverbatim 122* 123* Authors: 124* ======== 125* 126*> \author Univ. of Tennessee 127*> \author Univ. of California Berkeley 128*> \author Univ. of Colorado Denver 129*> \author NAG Ltd. 130* 131*> \par Further Details: 132* ===================== 133*> 134*> \verbatim 135*> Tall-Skinny QR (TSQR) performs QR by a sequence of orthogonal transformations, 136*> representing Q as a product of other orthogonal matrices 137*> Q = Q(1) * Q(2) * . . . * Q(k) 138*> where each Q(i) zeros out subdiagonal entries of a block of MB rows of A: 139*> Q(1) zeros out the subdiagonal entries of rows 1:MB of A 140*> Q(2) zeros out the bottom MB-N rows of rows [1:N,MB+1:2*MB-N] of A 141*> Q(3) zeros out the bottom MB-N rows of rows [1:N,2*MB-N+1:3*MB-2*N] of A 142*> . . . 143*> 144*> Q(1) is computed by GEQRT, which represents Q(1) by Householder vectors 145*> stored under the diagonal of rows 1:MB of A, and by upper triangular 146*> block reflectors, stored in array T(1:LDT,1:N). 147*> For more information see Further Details in GEQRT. 148*> 149*> Q(i) for i>1 is computed by TPQRT, which represents Q(i) by Householder vectors 150*> stored in rows [(i-1)*(MB-N)+N+1:i*(MB-N)+N] of A, and by upper triangular 151*> block reflectors, stored in array T(1:LDT,(i-1)*N+1:i*N). 152*> The last Q(k) may use fewer rows. 153*> For more information see Further Details in TPQRT. 154*> 155*> For more details of the overall algorithm, see the description of 156*> Sequential TSQR in Section 2.2 of [1]. 157*> 158*> [1] “Communication-Optimal Parallel and Sequential QR and LU Factorizations,” 159*> J. Demmel, L. Grigori, M. Hoemmen, J. Langou, 160*> SIAM J. Sci. Comput, vol. 34, no. 1, 2012 161*> \endverbatim 162*> 163* ===================================================================== 164 SUBROUTINE CLATSQR( M, N, MB, NB, A, LDA, T, LDT, WORK, 165 $ LWORK, INFO) 166* 167* -- LAPACK computational routine -- 168* -- LAPACK is a software package provided by Univ. of Tennessee, -- 169* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd. -- 170* 171* .. Scalar Arguments .. 172 INTEGER INFO, LDA, M, N, MB, NB, LDT, LWORK 173* .. 174* .. Array Arguments .. 175 COMPLEX A( LDA, * ), WORK( * ), T(LDT, *) 176* .. 177* 178* ===================================================================== 179* 180* .. 181* .. Local Scalars .. 182 LOGICAL LQUERY 183 INTEGER I, II, KK, CTR 184* .. 185* .. EXTERNAL FUNCTIONS .. 186 LOGICAL LSAME 187 EXTERNAL LSAME 188* .. EXTERNAL SUBROUTINES .. 189 EXTERNAL CGEQRT, CTPQRT, XERBLA 190* .. INTRINSIC FUNCTIONS .. 191 INTRINSIC MAX, MIN, MOD 192* .. 193* .. EXECUTABLE STATEMENTS .. 194* 195* TEST THE INPUT ARGUMENTS 196* 197 INFO = 0 198* 199 LQUERY = ( LWORK.EQ.-1 ) 200* 201 IF( M.LT.0 ) THEN 202 INFO = -1 203 ELSE IF( N.LT.0 .OR. M.LT.N ) THEN 204 INFO = -2 205 ELSE IF( MB.LE.N ) THEN 206 INFO = -3 207 ELSE IF( NB.LT.1 .OR. ( NB.GT.N .AND. N.GT.0 )) THEN 208 INFO = -4 209 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 210 INFO = -5 211 ELSE IF( LDT.LT.NB ) THEN 212 INFO = -8 213 ELSE IF( LWORK.LT.(N*NB) .AND. (.NOT.LQUERY) ) THEN 214 INFO = -10 215 END IF 216 IF( INFO.EQ.0) THEN 217 WORK(1) = NB*N 218 END IF 219 IF( INFO.NE.0 ) THEN 220 CALL XERBLA( 'CLATSQR', -INFO ) 221 RETURN 222 ELSE IF (LQUERY) THEN 223 RETURN 224 END IF 225* 226* Quick return if possible 227* 228 IF( MIN(M,N).EQ.0 ) THEN 229 RETURN 230 END IF 231* 232* The QR Decomposition 233* 234 IF ((MB.LE.N).OR.(MB.GE.M)) THEN 235 CALL CGEQRT( M, N, NB, A, LDA, T, LDT, WORK, INFO) 236 RETURN 237 END IF 238 KK = MOD((M-N),(MB-N)) 239 II=M-KK+1 240* 241* Compute the QR factorization of the first block A(1:MB,1:N) 242* 243 CALL CGEQRT( MB, N, NB, A(1,1), LDA, T, LDT, WORK, INFO ) 244 CTR = 1 245* 246 DO I = MB+1, II-MB+N , (MB-N) 247* 248* Compute the QR factorization of the current block A(I:I+MB-N,1:N) 249* 250 CALL CTPQRT( MB-N, N, 0, NB, A(1,1), LDA, A( I, 1 ), LDA, 251 $ T(1,CTR * N + 1), 252 $ LDT, WORK, INFO ) 253 CTR = CTR + 1 254 END DO 255* 256* Compute the QR factorization of the last block A(II:M,1:N) 257* 258 IF (II.LE.M) THEN 259 CALL CTPQRT( KK, N, 0, NB, A(1,1), LDA, A( II, 1 ), LDA, 260 $ T(1, CTR * N + 1), LDT, 261 $ WORK, INFO ) 262 END IF 263* 264 work( 1 ) = N*NB 265 RETURN 266* 267* End of CLATSQR 268* 269 END 270