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