001:       SUBROUTINE CGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
002:      $                   WORK, 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, 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: *  This routine is deprecated and has been replaced by routine CGELSY.
022: *
023: *  CGELSX computes the minimum-norm solution to a complex linear least
024: *  squares problem:
025: *      minimize || A * X - B ||
026: *  using a complete orthogonal factorization of A.  A is an M-by-N
027: *  matrix which may be rank-deficient.
028: *
029: *  Several right hand side vectors b and solution vectors x can be
030: *  handled in a single call; they are stored as the columns of the
031: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
032: *  matrix X.
033: *
034: *  The routine first computes a QR factorization with column pivoting:
035: *      A * P = Q * [ R11 R12 ]
036: *                  [  0  R22 ]
037: *  with R11 defined as the largest leading submatrix whose estimated
038: *  condition number is less than 1/RCOND.  The order of R11, RANK,
039: *  is the effective rank of A.
040: *
041: *  Then, R22 is considered to be negligible, and R12 is annihilated
042: *  by unitary transformations from the right, arriving at the
043: *  complete orthogonal factorization:
044: *     A * P = Q * [ T11 0 ] * Z
045: *                 [  0  0 ]
046: *  The minimum-norm solution is then
047: *     X = P * Z' [ inv(T11)*Q1'*B ]
048: *                [        0       ]
049: *  where Q1 consists of the first RANK columns of Q.
050: *
051: *  Arguments
052: *  =========
053: *
054: *  M       (input) INTEGER
055: *          The number of rows of the matrix A.  M >= 0.
056: *
057: *  N       (input) INTEGER
058: *          The number of columns of the matrix A.  N >= 0.
059: *
060: *  NRHS    (input) INTEGER
061: *          The number of right hand sides, i.e., the number of
062: *          columns of matrices B and X. NRHS >= 0.
063: *
064: *  A       (input/output) COMPLEX array, dimension (LDA,N)
065: *          On entry, the M-by-N matrix A.
066: *          On exit, A has been overwritten by details of its
067: *          complete orthogonal factorization.
068: *
069: *  LDA     (input) INTEGER
070: *          The leading dimension of the array A.  LDA >= max(1,M).
071: *
072: *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
073: *          On entry, the M-by-NRHS right hand side matrix B.
074: *          On exit, the N-by-NRHS solution matrix X.
075: *          If m >= n and RANK = n, the residual sum-of-squares for
076: *          the solution in the i-th column is given by the sum of
077: *          squares of elements N+1:M in that column.
078: *
079: *  LDB     (input) INTEGER
080: *          The leading dimension of the array B. LDB >= max(1,M,N).
081: *
082: *  JPVT    (input/output) INTEGER array, dimension (N)
083: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
084: *          initial column, otherwise it is a free column.  Before
085: *          the QR factorization of A, all initial columns are
086: *          permuted to the leading positions; only the remaining
087: *          free columns are moved as a result of column pivoting
088: *          during the factorization.
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) COMPLEX array, dimension
104: *                      (min(M,N) + max( N, 2*min(M,N)+NRHS )),
105: *
106: *  RWORK   (workspace) REAL array, dimension (2*N)
107: *
108: *  INFO    (output) INTEGER
109: *          = 0:  successful exit
110: *          < 0:  if INFO = -i, the i-th argument had an illegal value
111: *
112: *  =====================================================================
113: *
114: *     .. Parameters ..
115:       INTEGER            IMAX, IMIN
116:       PARAMETER          ( IMAX = 1, IMIN = 2 )
117:       REAL               ZERO, ONE, DONE, NTDONE
118:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, DONE = ZERO,
119:      $                   NTDONE = ONE )
120:       COMPLEX            CZERO, CONE
121:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
122:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
123: *     ..
124: *     .. Local Scalars ..
125:       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
126:       REAL               ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
127:      $                   SMLNUM
128:       COMPLEX            C1, C2, S1, S2, T1, T2
129: *     ..
130: *     .. External Subroutines ..
131:       EXTERNAL           CGEQPF, CLAIC1, CLASCL, CLASET, CLATZM, CTRSM,
132:      $                   CTZRQF, CUNM2R, SLABAD, XERBLA
133: *     ..
134: *     .. External Functions ..
135:       REAL               CLANGE, SLAMCH
136:       EXTERNAL           CLANGE, SLAMCH
137: *     ..
138: *     .. Intrinsic Functions ..
139:       INTRINSIC          ABS, CONJG, MAX, MIN
140: *     ..
141: *     .. Executable Statements ..
142: *
143:       MN = MIN( M, N )
144:       ISMIN = MN + 1
145:       ISMAX = 2*MN + 1
146: *
147: *     Test the input arguments.
148: *
149:       INFO = 0
150:       IF( M.LT.0 ) THEN
151:          INFO = -1
152:       ELSE IF( N.LT.0 ) THEN
153:          INFO = -2
154:       ELSE IF( NRHS.LT.0 ) THEN
155:          INFO = -3
156:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
157:          INFO = -5
158:       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
159:          INFO = -7
160:       END IF
161: *
162:       IF( INFO.NE.0 ) THEN
163:          CALL XERBLA( 'CGELSX', -INFO )
164:          RETURN
165:       END IF
166: *
167: *     Quick return if possible
168: *
169:       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
170:          RANK = 0
171:          RETURN
172:       END IF
173: *
174: *     Get machine parameters
175: *
176:       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
177:       BIGNUM = ONE / SMLNUM
178:       CALL SLABAD( SMLNUM, BIGNUM )
179: *
180: *     Scale A, B if max elements outside range [SMLNUM,BIGNUM]
181: *
182:       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
183:       IASCL = 0
184:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
185: *
186: *        Scale matrix norm up to SMLNUM
187: *
188:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
189:          IASCL = 1
190:       ELSE IF( ANRM.GT.BIGNUM ) THEN
191: *
192: *        Scale matrix norm down to BIGNUM
193: *
194:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
195:          IASCL = 2
196:       ELSE IF( ANRM.EQ.ZERO ) THEN
197: *
198: *        Matrix all zero. Return zero solution.
199: *
200:          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
201:          RANK = 0
202:          GO TO 100
203:       END IF
204: *
205:       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
206:       IBSCL = 0
207:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
208: *
209: *        Scale matrix norm up to SMLNUM
210: *
211:          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
212:          IBSCL = 1
213:       ELSE IF( BNRM.GT.BIGNUM ) THEN
214: *
215: *        Scale matrix norm down to BIGNUM
216: *
217:          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
218:          IBSCL = 2
219:       END IF
220: *
221: *     Compute QR factorization with column pivoting of A:
222: *        A * P = Q * R
223: *
224:       CALL CGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), RWORK,
225:      $             INFO )
226: *
227: *     complex workspace MN+N. Real workspace 2*N. Details of Householder
228: *     rotations stored in WORK(1:MN).
229: *
230: *     Determine RANK using incremental condition estimation
231: *
232:       WORK( ISMIN ) = CONE
233:       WORK( ISMAX ) = CONE
234:       SMAX = ABS( A( 1, 1 ) )
235:       SMIN = SMAX
236:       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
237:          RANK = 0
238:          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
239:          GO TO 100
240:       ELSE
241:          RANK = 1
242:       END IF
243: *
244:    10 CONTINUE
245:       IF( RANK.LT.MN ) THEN
246:          I = RANK + 1
247:          CALL CLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
248:      $                A( I, I ), SMINPR, S1, C1 )
249:          CALL CLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
250:      $                A( I, I ), SMAXPR, S2, C2 )
251: *
252:          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
253:             DO 20 I = 1, RANK
254:                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
255:                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
256:    20       CONTINUE
257:             WORK( ISMIN+RANK ) = C1
258:             WORK( ISMAX+RANK ) = C2
259:             SMIN = SMINPR
260:             SMAX = SMAXPR
261:             RANK = RANK + 1
262:             GO TO 10
263:          END IF
264:       END IF
265: *
266: *     Logically partition R = [ R11 R12 ]
267: *                             [  0  R22 ]
268: *     where R11 = R(1:RANK,1:RANK)
269: *
270: *     [R11,R12] = [ T11, 0 ] * Y
271: *
272:       IF( RANK.LT.N )
273:      $   CALL CTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
274: *
275: *     Details of Householder rotations stored in WORK(MN+1:2*MN)
276: *
277: *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
278: *
279:       CALL CUNM2R( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
280:      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), INFO )
281: *
282: *     workspace NRHS
283: *
284: *      B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
285: *
286:       CALL CTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
287:      $            NRHS, CONE, A, LDA, B, LDB )
288: *
289:       DO 40 I = RANK + 1, N
290:          DO 30 J = 1, NRHS
291:             B( I, J ) = CZERO
292:    30    CONTINUE
293:    40 CONTINUE
294: *
295: *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
296: *
297:       IF( RANK.LT.N ) THEN
298:          DO 50 I = 1, RANK
299:             CALL CLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
300:      $                   CONJG( WORK( MN+I ) ), B( I, 1 ),
301:      $                   B( RANK+1, 1 ), LDB, WORK( 2*MN+1 ) )
302:    50    CONTINUE
303:       END IF
304: *
305: *     workspace NRHS
306: *
307: *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
308: *
309:       DO 90 J = 1, NRHS
310:          DO 60 I = 1, N
311:             WORK( 2*MN+I ) = NTDONE
312:    60    CONTINUE
313:          DO 80 I = 1, N
314:             IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
315:                IF( JPVT( I ).NE.I ) THEN
316:                   K = I
317:                   T1 = B( K, J )
318:                   T2 = B( JPVT( K ), J )
319:    70             CONTINUE
320:                   B( JPVT( K ), J ) = T1
321:                   WORK( 2*MN+K ) = DONE
322:                   T1 = T2
323:                   K = JPVT( K )
324:                   T2 = B( JPVT( K ), J )
325:                   IF( JPVT( K ).NE.I )
326:      $               GO TO 70
327:                   B( I, J ) = T1
328:                   WORK( 2*MN+K ) = DONE
329:                END IF
330:             END IF
331:    80    CONTINUE
332:    90 CONTINUE
333: *
334: *     Undo scaling
335: *
336:       IF( IASCL.EQ.1 ) THEN
337:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
338:          CALL CLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
339:      $                INFO )
340:       ELSE IF( IASCL.EQ.2 ) THEN
341:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
342:          CALL CLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
343:      $                INFO )
344:       END IF
345:       IF( IBSCL.EQ.1 ) THEN
346:          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
347:       ELSE IF( IBSCL.EQ.2 ) THEN
348:          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
349:       END IF
350: *
351:   100 CONTINUE
352: *
353:       RETURN
354: *
355: *     End of CGELSX
356: *
357:       END
358: