```001:       SUBROUTINE CGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
002:      \$                  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:       CHARACTER          TRANS
010:       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  CGELS solves overdetermined or underdetermined complex linear systems
020: *  involving an M-by-N matrix A, or its conjugate-transpose, using a QR
021: *  or LQ factorization of A.  It is assumed that A has full rank.
022: *
023: *  The following options are provided:
024: *
025: *  1. If TRANS = 'N' and m >= n:  find the least squares solution of
026: *     an overdetermined system, i.e., solve the least squares problem
027: *                  minimize || B - A*X ||.
028: *
029: *  2. If TRANS = 'N' and m < n:  find the minimum norm solution of
030: *     an underdetermined system A * X = B.
031: *
032: *  3. If TRANS = 'C' and m >= n:  find the minimum norm solution of
033: *     an undetermined system A**H * X = B.
034: *
035: *  4. If TRANS = 'C' and m < n:  find the least squares solution of
036: *     an overdetermined system, i.e., solve the least squares problem
037: *                  minimize || B - A**H * X ||.
038: *
039: *  Several right hand side vectors b and solution vectors x can be
040: *  handled in a single call; they are stored as the columns of the
041: *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
042: *  matrix X.
043: *
044: *  Arguments
045: *  =========
046: *
047: *  TRANS   (input) CHARACTER*1
048: *          = 'N': the linear system involves A;
049: *          = 'C': the linear system involves A**H.
050: *
051: *  M       (input) INTEGER
052: *          The number of rows of the matrix A.  M >= 0.
053: *
054: *  N       (input) INTEGER
055: *          The number of columns of the matrix A.  N >= 0.
056: *
057: *  NRHS    (input) INTEGER
058: *          The number of right hand sides, i.e., the number of
059: *          columns of the matrices B and X. NRHS >= 0.
060: *
061: *  A       (input/output) COMPLEX array, dimension (LDA,N)
062: *          On entry, the M-by-N matrix A.
063: *            if M >= N, A is overwritten by details of its QR
064: *                       factorization as returned by CGEQRF;
065: *            if M <  N, A is overwritten by details of its LQ
066: *                       factorization as returned by CGELQF.
067: *
068: *  LDA     (input) INTEGER
069: *          The leading dimension of the array A.  LDA >= max(1,M).
070: *
071: *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
072: *          On entry, the matrix B of right hand side vectors, stored
073: *          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
074: *          if TRANS = 'C'.
075: *          On exit, if INFO = 0, B is overwritten by the solution
076: *          vectors, stored columnwise:
077: *          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
078: *          squares solution vectors; the residual sum of squares for the
079: *          solution in each column is given by the sum of squares of the
080: *          modulus of elements N+1 to M in that column;
081: *          if TRANS = 'N' and m < n, rows 1 to N of B contain the
082: *          minimum norm solution vectors;
083: *          if TRANS = 'C' and m >= n, rows 1 to M of B contain the
084: *          minimum norm solution vectors;
085: *          if TRANS = 'C' and m < n, rows 1 to M of B contain the
086: *          least squares solution vectors; the residual sum of squares
087: *          for the solution in each column is given by the sum of
088: *          squares of the modulus of elements M+1 to N in that column.
089: *
090: *  LDB     (input) INTEGER
091: *          The leading dimension of the array B. LDB >= MAX(1,M,N).
092: *
093: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
094: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
095: *
096: *  LWORK   (input) INTEGER
097: *          The dimension of the array WORK.
098: *          LWORK >= max( 1, MN + max( MN, NRHS ) ).
099: *          For optimal performance,
100: *          LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
101: *          where MN = min(M,N) and NB is the optimum block size.
102: *
103: *          If LWORK = -1, then a workspace query is assumed; the routine
104: *          only calculates the optimal size of the WORK array, returns
105: *          this value as the first entry of the WORK array, and no error
106: *          message related to LWORK is issued by XERBLA.
107: *
108: *  INFO    (output) INTEGER
109: *          = 0:  successful exit
110: *          < 0:  if INFO = -i, the i-th argument had an illegal value
111: *          > 0:  if INFO =  i, the i-th diagonal element of the
112: *                triangular factor of A is zero, so that A does not have
113: *                full rank; the least squares solution could not be
114: *                computed.
115: *
116: *  =====================================================================
117: *
118: *     .. Parameters ..
119:       REAL               ZERO, ONE
120:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
121:       COMPLEX            CZERO
122:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
123: *     ..
124: *     .. Local Scalars ..
125:       LOGICAL            LQUERY, TPSD
126:       INTEGER            BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
127:       REAL               ANRM, BIGNUM, BNRM, SMLNUM
128: *     ..
129: *     .. Local Arrays ..
130:       REAL               RWORK( 1 )
131: *     ..
132: *     .. External Functions ..
133:       LOGICAL            LSAME
134:       INTEGER            ILAENV
135:       REAL               CLANGE, SLAMCH
136:       EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
137: *     ..
138: *     .. External Subroutines ..
139:       EXTERNAL           CGELQF, CGEQRF, CLASCL, CLASET, CTRTRS, CUNMLQ,
141: *     ..
142: *     .. Intrinsic Functions ..
143:       INTRINSIC          MAX, MIN, REAL
144: *     ..
145: *     .. Executable Statements ..
146: *
147: *     Test the input arguments.
148: *
149:       INFO = 0
150:       MN = MIN( M, N )
151:       LQUERY = ( LWORK.EQ.-1 )
152:       IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'C' ) ) ) THEN
153:          INFO = -1
154:       ELSE IF( M.LT.0 ) THEN
155:          INFO = -2
156:       ELSE IF( N.LT.0 ) THEN
157:          INFO = -3
158:       ELSE IF( NRHS.LT.0 ) THEN
159:          INFO = -4
160:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
161:          INFO = -6
162:       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
163:          INFO = -8
164:       ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND.
165:      \$   .NOT.LQUERY ) THEN
166:          INFO = -10
167:       END IF
168: *
169: *     Figure out optimal block size
170: *
171:       IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
172: *
173:          TPSD = .TRUE.
174:          IF( LSAME( TRANS, 'N' ) )
175:      \$      TPSD = .FALSE.
176: *
177:          IF( M.GE.N ) THEN
178:             NB = ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
179:             IF( TPSD ) THEN
180:                NB = MAX( NB, ILAENV( 1, 'CUNMQR', 'LN', M, NRHS, N,
181:      \$              -1 ) )
182:             ELSE
183:                NB = MAX( NB, ILAENV( 1, 'CUNMQR', 'LC', M, NRHS, N,
184:      \$              -1 ) )
185:             END IF
186:          ELSE
187:             NB = ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
188:             IF( TPSD ) THEN
189:                NB = MAX( NB, ILAENV( 1, 'CUNMLQ', 'LC', N, NRHS, M,
190:      \$              -1 ) )
191:             ELSE
192:                NB = MAX( NB, ILAENV( 1, 'CUNMLQ', 'LN', N, NRHS, M,
193:      \$              -1 ) )
194:             END IF
195:          END IF
196: *
197:          WSIZE = MAX( 1, MN + MAX( MN, NRHS )*NB )
198:          WORK( 1 ) = REAL( WSIZE )
199: *
200:       END IF
201: *
202:       IF( INFO.NE.0 ) THEN
203:          CALL XERBLA( 'CGELS ', -INFO )
204:          RETURN
205:       ELSE IF( LQUERY ) THEN
206:          RETURN
207:       END IF
208: *
209: *     Quick return if possible
210: *
211:       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
212:          CALL CLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
213:          RETURN
214:       END IF
215: *
216: *     Get machine parameters
217: *
218:       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
219:       BIGNUM = ONE / SMLNUM
220:       CALL SLABAD( SMLNUM, BIGNUM )
221: *
222: *     Scale A, B if max element outside range [SMLNUM,BIGNUM]
223: *
224:       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
225:       IASCL = 0
226:       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
227: *
228: *        Scale matrix norm up to SMLNUM
229: *
230:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
231:          IASCL = 1
232:       ELSE IF( ANRM.GT.BIGNUM ) THEN
233: *
234: *        Scale matrix norm down to BIGNUM
235: *
236:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
237:          IASCL = 2
238:       ELSE IF( ANRM.EQ.ZERO ) THEN
239: *
240: *        Matrix all zero. Return zero solution.
241: *
242:          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
243:          GO TO 50
244:       END IF
245: *
246:       BROW = M
247:       IF( TPSD )
248:      \$   BROW = N
249:       BNRM = CLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
250:       IBSCL = 0
251:       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
252: *
253: *        Scale matrix norm up to SMLNUM
254: *
255:          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
256:      \$                INFO )
257:          IBSCL = 1
258:       ELSE IF( BNRM.GT.BIGNUM ) THEN
259: *
260: *        Scale matrix norm down to BIGNUM
261: *
262:          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
263:      \$                INFO )
264:          IBSCL = 2
265:       END IF
266: *
267:       IF( M.GE.N ) THEN
268: *
269: *        compute QR factorization of A
270: *
271:          CALL CGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
272:      \$                INFO )
273: *
274: *        workspace at least N, optimally N*NB
275: *
276:          IF( .NOT.TPSD ) THEN
277: *
278: *           Least-Squares Problem min || A * X - B ||
279: *
280: *           B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
281: *
282:             CALL CUNMQR( 'Left', 'Conjugate transpose', M, NRHS, N, A,
283:      \$                   LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
284:      \$                   INFO )
285: *
286: *           workspace at least NRHS, optimally NRHS*NB
287: *
288: *           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
289: *
290:             CALL CTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
291:      \$                   A, LDA, B, LDB, INFO )
292: *
293:             IF( INFO.GT.0 ) THEN
294:                RETURN
295:             END IF
296: *
297:             SCLLEN = N
298: *
299:          ELSE
300: *
301: *           Overdetermined system of equations A' * X = B
302: *
303: *           B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS)
304: *
305:             CALL CTRTRS( 'Upper', 'Conjugate transpose','Non-unit',
306:      \$                   N, NRHS, A, LDA, B, LDB, INFO )
307: *
308:             IF( INFO.GT.0 ) THEN
309:                RETURN
310:             END IF
311: *
312: *           B(N+1:M,1:NRHS) = ZERO
313: *
314:             DO 20 J = 1, NRHS
315:                DO 10 I = N + 1, M
316:                   B( I, J ) = CZERO
317:    10          CONTINUE
318:    20       CONTINUE
319: *
320: *           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
321: *
322:             CALL CUNMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
323:      \$                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
324:      \$                   INFO )
325: *
326: *           workspace at least NRHS, optimally NRHS*NB
327: *
328:             SCLLEN = M
329: *
330:          END IF
331: *
332:       ELSE
333: *
334: *        Compute LQ factorization of A
335: *
336:          CALL CGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
337:      \$                INFO )
338: *
339: *        workspace at least M, optimally M*NB.
340: *
341:          IF( .NOT.TPSD ) THEN
342: *
343: *           underdetermined system of equations A * X = B
344: *
345: *           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
346: *
347:             CALL CTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
348:      \$                   A, LDA, B, LDB, INFO )
349: *
350:             IF( INFO.GT.0 ) THEN
351:                RETURN
352:             END IF
353: *
354: *           B(M+1:N,1:NRHS) = 0
355: *
356:             DO 40 J = 1, NRHS
357:                DO 30 I = M + 1, N
358:                   B( I, J ) = CZERO
359:    30          CONTINUE
360:    40       CONTINUE
361: *
362: *           B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS)
363: *
364:             CALL CUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A,
365:      \$                   LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
366:      \$                   INFO )
367: *
368: *           workspace at least NRHS, optimally NRHS*NB
369: *
370:             SCLLEN = N
371: *
372:          ELSE
373: *
374: *           overdetermined system min || A' * X - B ||
375: *
376: *           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
377: *
378:             CALL CUNMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
379:      \$                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
380:      \$                   INFO )
381: *
382: *           workspace at least NRHS, optimally NRHS*NB
383: *
384: *           B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS)
385: *
386:             CALL CTRTRS( 'Lower', 'Conjugate transpose', 'Non-unit',
387:      \$                   M, NRHS, A, LDA, B, LDB, INFO )
388: *
389:             IF( INFO.GT.0 ) THEN
390:                RETURN
391:             END IF
392: *
393:             SCLLEN = M
394: *
395:          END IF
396: *
397:       END IF
398: *
399: *     Undo scaling
400: *
401:       IF( IASCL.EQ.1 ) THEN
402:          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
403:      \$                INFO )
404:       ELSE IF( IASCL.EQ.2 ) THEN
405:          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
406:      \$                INFO )
407:       END IF
408:       IF( IBSCL.EQ.1 ) THEN
409:          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
410:      \$                INFO )
411:       ELSE IF( IBSCL.EQ.2 ) THEN
412:          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
413:      \$                INFO )
414:       END IF
415: *
416:    50 CONTINUE
417:       WORK( 1 ) = REAL( WSIZE )
418: *
419:       RETURN
420: *
421: *     End of CGELS
422: *
423:       END
424: ```