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