162
163
164
165
166
167
168 CHARACTER TRANS
169 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
170
171
172 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
173
174
175
176
177
178
179 REAL ZERO, ONE
180 parameter( zero = 0.0e0, one = 1.0e0 )
181 COMPLEX CZERO
182 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
183
184
185 LOGICAL LQUERY, TRAN
186 INTEGER I, IASCL, IBSCL, J, MAXMN, BROW,
187 $ SCLLEN, TSZO, TSZM, LWO, LWM, LW1, LW2,
188 $ WSIZEO, WSIZEM, INFO2
189 REAL ANRM, BIGNUM, BNRM, SMLNUM, DUM( 1 )
190 COMPLEX TQ( 5 ), WORKQ( 1 )
191
192
193 LOGICAL LSAME
194 REAL SLAMCH, CLANGE, SROUNDUP_LWORK
196
197
200
201
202 INTRINSIC max, min, int
203
204
205
206
207
208 info = 0
209 maxmn = max( m, n )
210 tran =
lsame( trans,
'C' )
211
212 lquery = ( lwork.EQ.-1 .OR. lwork.EQ.-2 )
213 IF( .NOT.(
lsame( trans,
'N' ) .OR.
214 $
lsame( trans,
'C' ) ) )
THEN
215 info = -1
216 ELSE IF( m.LT.0 ) THEN
217 info = -2
218 ELSE IF( n.LT.0 ) THEN
219 info = -3
220 ELSE IF( nrhs.LT.0 ) THEN
221 info = -4
222 ELSE IF( lda.LT.max( 1, m ) ) THEN
223 info = -6
224 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
225 info = -8
226 END IF
227
228 IF( info.EQ.0 ) THEN
229
230
231
232 IF( m.GE.n ) THEN
233 CALL cgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
234 tszo = int( tq( 1 ) )
235 lwo = int( workq( 1 ) )
236 CALL cgemqr(
'L', trans, m, nrhs, n, a, lda, tq,
237 $ tszo, b, ldb, workq, -1, info2 )
238 lwo = max( lwo, int( workq( 1 ) ) )
239 CALL cgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
240 tszm = int( tq( 1 ) )
241 lwm = int( workq( 1 ) )
242 CALL cgemqr(
'L', trans, m, nrhs, n, a, lda, tq,
243 $ tszm, b, ldb, workq, -1, info2 )
244 lwm = max( lwm, int( workq( 1 ) ) )
245 wsizeo = tszo + lwo
246 wsizem = tszm + lwm
247 ELSE
248 CALL cgelq( m, n, a, lda, tq, -1, workq, -1, info2 )
249 tszo = int( tq( 1 ) )
250 lwo = int( workq( 1 ) )
251 CALL cgemlq(
'L', trans, n, nrhs, m, a, lda, tq,
252 $ tszo, b, ldb, workq, -1, info2 )
253 lwo = max( lwo, int( workq( 1 ) ) )
254 CALL cgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
255 tszm = int( tq( 1 ) )
256 lwm = int( workq( 1 ) )
257 CALL cgemlq(
'L', trans, n, nrhs, m, a, lda, tq,
258 $ tszm, b, ldb, workq, -1, info2 )
259 lwm = max( lwm, int( workq( 1 ) ) )
260 wsizeo = tszo + lwo
261 wsizem = tszm + lwm
262 END IF
263
264 IF( ( lwork.LT.wsizem ).AND.( .NOT.lquery ) ) THEN
265 info = -10
266 END IF
267
269
270 END IF
271
272 IF( info.NE.0 ) THEN
273 CALL xerbla(
'CGETSLS', -info )
274 RETURN
275 END IF
276 IF( lquery ) THEN
278 RETURN
279 END IF
280 IF( lwork.LT.wsizeo ) THEN
281 lw1 = tszm
282 lw2 = lwm
283 ELSE
284 lw1 = tszo
285 lw2 = lwo
286 END IF
287
288
289
290 IF( min( m, n, nrhs ).EQ.0 ) THEN
291 CALL claset(
'FULL', max( m, n ), nrhs, czero, czero,
292 $ b, ldb )
293 RETURN
294 END IF
295
296
297
299 bignum = one / smlnum
300
301
302
303 anrm =
clange(
'M', m, n, a, lda, dum )
304 iascl = 0
305 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
306
307
308
309 CALL clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
310 iascl = 1
311 ELSE IF( anrm.GT.bignum ) THEN
312
313
314
315 CALL clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
316 iascl = 2
317 ELSE IF( anrm.EQ.zero ) THEN
318
319
320
321 CALL claset(
'F', maxmn, nrhs, czero, czero, b, ldb )
322 GO TO 50
323 END IF
324
325 brow = m
326 IF ( tran ) THEN
327 brow = n
328 END IF
329 bnrm =
clange(
'M', brow, nrhs, b, ldb, dum )
330 ibscl = 0
331 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
332
333
334
335 CALL clascl(
'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
336 $ info )
337 ibscl = 1
338 ELSE IF( bnrm.GT.bignum ) THEN
339
340
341
342 CALL clascl(
'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
343 $ info )
344 ibscl = 2
345 END IF
346
347 IF ( m.GE.n ) THEN
348
349
350
351 CALL cgeqr( m, n, a, lda, work( lw2+1 ), lw1,
352 $ work( 1 ), lw2, info )
353 IF ( .NOT.tran ) THEN
354
355
356
357
358
359 CALL cgemqr(
'L' ,
'C', m, nrhs, n, a, lda,
360 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
361 $ info )
362
363
364
365 CALL ctrtrs(
'U',
'N',
'N', n, nrhs,
366 $ a, lda, b, ldb, info )
367 IF( info.GT.0 ) THEN
368 RETURN
369 END IF
370 scllen = n
371 ELSE
372
373
374
375
376
377 CALL ctrtrs(
'U',
'C',
'N', n, nrhs,
378 $ a, lda, b, ldb, info )
379
380 IF( info.GT.0 ) THEN
381 RETURN
382 END IF
383
384
385
386 DO 20 j = 1, nrhs
387 DO 10 i = n + 1, m
388 b( i, j ) = czero
389 10 CONTINUE
390 20 CONTINUE
391
392
393
394 CALL cgemqr(
'L',
'N', m, nrhs, n, a, lda,
395 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
396 $ info )
397
398 scllen = m
399
400 END IF
401
402 ELSE
403
404
405
406 CALL cgelq( m, n, a, lda, work( lw2+1 ), lw1,
407 $ work( 1 ), lw2, info )
408
409
410
411 IF( .NOT.tran ) THEN
412
413
414
415
416
417 CALL ctrtrs(
'L',
'N',
'N', m, nrhs,
418 $ a, lda, b, ldb, info )
419
420 IF( info.GT.0 ) THEN
421 RETURN
422 END IF
423
424
425
426 DO 40 j = 1, nrhs
427 DO 30 i = m + 1, n
428 b( i, j ) = czero
429 30 CONTINUE
430 40 CONTINUE
431
432
433
434 CALL cgemlq(
'L',
'C', n, nrhs, m, a, lda,
435 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
436 $ info )
437
438
439
440 scllen = n
441
442 ELSE
443
444
445
446
447
448 CALL cgemlq(
'L',
'N', n, nrhs, m, a, lda,
449 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
450 $ info )
451
452
453
454
455
456 CALL ctrtrs(
'L',
'C',
'N', m, nrhs,
457 $ a, lda, b, ldb, info )
458
459 IF( info.GT.0 ) THEN
460 RETURN
461 END IF
462
463 scllen = m
464
465 END IF
466
467 END IF
468
469
470
471 IF( iascl.EQ.1 ) THEN
472 CALL clascl(
'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
473 $ info )
474 ELSE IF( iascl.EQ.2 ) THEN
475 CALL clascl(
'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
476 $ info )
477 END IF
478 IF( ibscl.EQ.1 ) THEN
479 CALL clascl(
'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
480 $ info )
481 ELSE IF( ibscl.EQ.2 ) THEN
482 CALL clascl(
'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
483 $ info )
484 END IF
485
486 50 CONTINUE
488 RETURN
489
490
491
subroutine xerbla(srname, info)
subroutine cgelq(m, n, a, lda, t, tsize, work, lwork, info)
CGELQ
subroutine cgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
CGEMLQ
subroutine cgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
CGEMQR
subroutine cgeqr(m, n, a, lda, t, tsize, work, lwork, info)
CGEQR
real function slamch(cmach)
SLAMCH
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine ctrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
CTRTRS