001:       SUBROUTINE CGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
002:      $                   WORK, LWORK, RWORK, INFO )
003: *
004: *  -- LAPACK driver routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
010:       REAL               RCOND
011: *     ..
012: *     .. Array Arguments ..
013:       INTEGER            JPVT( * )
014:       REAL               RWORK( * )
015:       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CGELSY computes the minimum-norm solution to a complex linear least
022: *  squares problem:
023: *      minimize || A * X - B ||
024: *  using a complete orthogonal factorization of A.  A is an M-by-N
025: *  matrix which may be rank-deficient.
026: *
027: *  Several right hand side vectors b and solution vectors x can be
028: *  handled in a single call; they are stored as the columns of the
029: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
030: *  matrix X.
031: *
032: *  The routine first computes a QR factorization with column pivoting:
033: *      A * P = Q * [ R11 R12 ]
034: *                  [  0  R22 ]
035: *  with R11 defined as the largest leading submatrix whose estimated
036: *  condition number is less than 1/RCOND.  The order of R11, RANK,
037: *  is the effective rank of A.
038: *
039: *  Then, R22 is considered to be negligible, and R12 is annihilated
040: *  by unitary transformations from the right, arriving at the
041: *  complete orthogonal factorization:
042: *     A * P = Q * [ T11 0 ] * Z
043: *                 [  0  0 ]
044: *  The minimum-norm solution is then
045: *     X = P * Z' [ inv(T11)*Q1'*B ]
046: *                [        0       ]
047: *  where Q1 consists of the first RANK columns of Q.
048: *
049: *  This routine is basically identical to the original xGELSX except
050: *  three differences:
051: *    o The permutation of matrix B (the right hand side) is faster and
052: *      more simple.
053: *    o The call to the subroutine xGEQPF has been substituted by the
054: *      the call to the subroutine xGEQP3. This subroutine is a Blas-3
055: *      version of the QR factorization with column pivoting.
056: *    o Matrix B (the right hand side) is updated with Blas-3.
057: *
058: *  Arguments
059: *  =========
060: *
061: *  M       (input) INTEGER
062: *          The number of rows of the matrix A.  M >= 0.
063: *
064: *  N       (input) INTEGER
065: *          The number of columns of the matrix A.  N >= 0.
066: *
067: *  NRHS    (input) INTEGER
068: *          The number of right hand sides, i.e., the number of
069: *          columns of matrices B and X. NRHS >= 0.
070: *
071: *  A       (input/output) COMPLEX array, dimension (LDA,N)
072: *          On entry, the M-by-N matrix A.
073: *          On exit, A has been overwritten by details of its
074: *          complete orthogonal factorization.
075: *
076: *  LDA     (input) INTEGER
077: *          The leading dimension of the array A.  LDA >= max(1,M).
078: *
079: *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
080: *          On entry, the M-by-NRHS right hand side matrix B.
081: *          On exit, the N-by-NRHS solution matrix X.
082: *
083: *  LDB     (input) INTEGER
084: *          The leading dimension of the array B. LDB >= max(1,M,N).
085: *
086: *  JPVT    (input/output) INTEGER array, dimension (N)
087: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
088: *          to the front of AP, otherwise column i is a free column.
089: *          On exit, if JPVT(i) = k, then the i-th column of A*P
090: *          was the k-th column of A.
091: *
092: *  RCOND   (input) REAL
093: *          RCOND is used to determine the effective rank of A, which
094: *          is defined as the order of the largest leading triangular
095: *          submatrix R11 in the QR factorization with pivoting of A,
096: *          whose estimated condition number < 1/RCOND.
097: *
098: *  RANK    (output) INTEGER
099: *          The effective rank of A, i.e., the order of the submatrix
100: *          R11.  This is the same as the order of the submatrix T11
101: *          in the complete orthogonal factorization of A.
102: *
103: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
104: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
105: *
106: *  LWORK   (input) INTEGER
107: *          The dimension of the array WORK.
108: *          The unblocked strategy requires that:
109: *            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
110: *          where MN = min(M,N).
111: *          The block algorithm requires that:
112: *            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
113: *          where NB is an upper bound on the blocksize returned
114: *          by ILAENV for the routines CGEQP3, CTZRZF, CTZRQF, CUNMQR,
115: *          and CUNMRZ.
116: *
117: *          If LWORK = -1, then a workspace query is assumed; the routine
118: *          only calculates the optimal size of the WORK array, returns
119: *          this value as the first entry of the WORK array, and no error
120: *          message related to LWORK is issued by XERBLA.
121: *
122: *  RWORK   (workspace) REAL array, dimension (2*N)
123: *
124: *  INFO    (output) INTEGER
125: *          = 0: successful exit
126: *          < 0: if INFO = -i, the i-th argument had an illegal value
127: *
128: *  Further Details
129: *  ===============
130: *
131: *  Based on contributions by
132: *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
133: *    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
134: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
135: *
136: *  =====================================================================
137: *
138: *     .. Parameters ..
139:       INTEGER            IMAX, IMIN
140:       PARAMETER          ( IMAX = 1, IMIN = 2 )
141:       REAL               ZERO, ONE
142:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
143:       COMPLEX            CZERO, CONE
144:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
145:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
146: *     ..
147: *     .. Local Scalars ..
148:       LOGICAL            LQUERY
149:       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
150:      $                   NB, NB1, NB2, NB3, NB4
151:       REAL               ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
152:      $                   SMLNUM, WSIZE
153:       COMPLEX            C1, C2, S1, S2
154: *     ..
155: *     .. External Subroutines ..
156:       EXTERNAL           CCOPY, CGEQP3, CLAIC1, CLASCL, CLASET, CTRSM,
157:      $                   CTZRZF, CUNMQR, CUNMRZ, SLABAD, XERBLA
158: *     ..
159: *     .. External Functions ..
160:       INTEGER            ILAENV
161:       REAL               CLANGE, SLAMCH
162:       EXTERNAL           CLANGE, ILAENV, SLAMCH
163: *     ..
164: *     .. Intrinsic Functions ..
165:       INTRINSIC          ABS, MAX, MIN, REAL, CMPLX
166: *     ..
167: *     .. Executable Statements ..
168: *
169:       MN = MIN( M, N )
170:       ISMIN = MN + 1
171:       ISMAX = 2*MN + 1
172: *
173: *     Test the input arguments.
174: *
175:       INFO = 0
176:       NB1 = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
177:       NB2 = ILAENV( 1, 'CGERQF', ' ', M, N, -1, -1 )
178:       NB3 = ILAENV( 1, 'CUNMQR', ' ', M, N, NRHS, -1 )
179:       NB4 = ILAENV( 1, 'CUNMRQ', ' ', M, N, NRHS, -1 )
180:       NB = MAX( NB1, NB2, NB3, NB4 )
181:       LWKOPT = MAX( 1, MN+2*N+NB*(N+1), 2*MN+NB*NRHS )
182:       WORK( 1 ) = CMPLX( LWKOPT )
183:       LQUERY = ( LWORK.EQ.-1 )
184:       IF( M.LT.0 ) THEN
185:          INFO = -1
186:       ELSE IF( N.LT.0 ) THEN
187:          INFO = -2
188:       ELSE IF( NRHS.LT.0 ) THEN
189:          INFO = -3
190:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
191:          INFO = -5
192:       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
193:          INFO = -7
194:       ELSE IF( LWORK.LT.( MN+MAX( 2*MN, N+1, MN+NRHS ) ) .AND.
195:      $   .NOT.LQUERY ) THEN
196:          INFO = -12
197:       END IF
198: *
199:       IF( INFO.NE.0 ) THEN
200:          CALL XERBLA( 'CGELSY', -INFO )
201:          RETURN
202:       ELSE IF( LQUERY ) THEN
203:          RETURN
204:       END IF
205: *
206: *     Quick return if possible
207: *
208:       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
209:          RANK = 0
210:          RETURN
211:       END IF
212: *
213: *     Get machine parameters
214: *
215:       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
216:       BIGNUM = ONE / SMLNUM
217:       CALL SLABAD( SMLNUM, BIGNUM )
218: *
219: *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
220: *
221:       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
222:       IASCL = 0
223:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
224: *
225: *        Scale matrix norm up to SMLNUM
226: *
227:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
228:          IASCL = 1
229:       ELSE IF( ANRM.GT.BIGNUM ) THEN
230: *
231: *        Scale matrix norm down to BIGNUM
232: *
233:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
234:          IASCL = 2
235:       ELSE IF( ANRM.EQ.ZERO ) THEN
236: *
237: *        Matrix all zero. Return zero solution.
238: *
239:          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
240:          RANK = 0
241:          GO TO 70
242:       END IF
243: *
244:       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
245:       IBSCL = 0
246:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
247: *
248: *        Scale matrix norm up to SMLNUM
249: *
250:          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
251:          IBSCL = 1
252:       ELSE IF( BNRM.GT.BIGNUM ) THEN
253: *
254: *        Scale matrix norm down to BIGNUM
255: *
256:          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
257:          IBSCL = 2
258:       END IF
259: *
260: *     Compute QR factorization with column pivoting of A:
261: *        A * P = Q * R
262: *
263:       CALL CGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
264:      $             LWORK-MN, RWORK, INFO )
265:       WSIZE = MN + REAL( WORK( MN+1 ) )
266: *
267: *     complex workspace: MN+NB*(N+1). real workspace 2*N.
268: *     Details of Householder rotations stored in WORK(1:MN).
269: *
270: *     Determine RANK using incremental condition estimation
271: *
272:       WORK( ISMIN ) = CONE
273:       WORK( ISMAX ) = CONE
274:       SMAX = ABS( A( 1, 1 ) )
275:       SMIN = SMAX
276:       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
277:          RANK = 0
278:          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
279:          GO TO 70
280:       ELSE
281:          RANK = 1
282:       END IF
283: *
284:    10 CONTINUE
285:       IF( RANK.LT.MN ) THEN
286:          I = RANK + 1
287:          CALL CLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
288:      $                A( I, I ), SMINPR, S1, C1 )
289:          CALL CLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
290:      $                A( I, I ), SMAXPR, S2, C2 )
291: *
292:          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
293:             DO 20 I = 1, RANK
294:                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
295:                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
296:    20       CONTINUE
297:             WORK( ISMIN+RANK ) = C1
298:             WORK( ISMAX+RANK ) = C2
299:             SMIN = SMINPR
300:             SMAX = SMAXPR
301:             RANK = RANK + 1
302:             GO TO 10
303:          END IF
304:       END IF
305: *
306: *     complex workspace: 3*MN.
307: *
308: *     Logically partition R = [ R11 R12 ]
309: *                             [  0  R22 ]
310: *     where R11 = R(1:RANK,1:RANK)
311: *
312: *     [R11,R12] = [ T11, 0 ] * Y
313: *
314:       IF( RANK.LT.N )
315:      $   CALL CTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
316:      $                LWORK-2*MN, INFO )
317: *
318: *     complex workspace: 2*MN.
319: *     Details of Householder rotations stored in WORK(MN+1:2*MN)
320: *
321: *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
322: *
323:       CALL CUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
324:      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
325:       WSIZE = MAX( WSIZE, 2*MN+REAL( WORK( 2*MN+1 ) ) )
326: *
327: *     complex workspace: 2*MN+NB*NRHS.
328: *
329: *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
330: *
331:       CALL CTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
332:      $            NRHS, CONE, A, LDA, B, LDB )
333: *
334:       DO 40 J = 1, NRHS
335:          DO 30 I = RANK + 1, N
336:             B( I, J ) = CZERO
337:    30    CONTINUE
338:    40 CONTINUE
339: *
340: *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
341: *
342:       IF( RANK.LT.N ) THEN
343:          CALL CUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK,
344:      $                N-RANK, A, LDA, WORK( MN+1 ), B, LDB,
345:      $                WORK( 2*MN+1 ), LWORK-2*MN, INFO )
346:       END IF
347: *
348: *     complex workspace: 2*MN+NRHS.
349: *
350: *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
351: *
352:       DO 60 J = 1, NRHS
353:          DO 50 I = 1, N
354:             WORK( JPVT( I ) ) = B( I, J )
355:    50    CONTINUE
356:          CALL CCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
357:    60 CONTINUE
358: *
359: *     complex workspace: N.
360: *
361: *     Undo scaling
362: *
363:       IF( IASCL.EQ.1 ) THEN
364:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
365:          CALL CLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
366:      $                INFO )
367:       ELSE IF( IASCL.EQ.2 ) THEN
368:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
369:          CALL CLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
370:      $                INFO )
371:       END IF
372:       IF( IBSCL.EQ.1 ) THEN
373:          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
374:       ELSE IF( IBSCL.EQ.2 ) THEN
375:          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
376:       END IF
377: *
378:    70 CONTINUE
379:       WORK( 1 ) = CMPLX( LWKOPT )
380: *
381:       RETURN
382: *
383: *     End of CGELSY
384: *
385:       END
386: