213
214
215
216
217
218
219 CHARACTER JOBU, JOBVT
220 INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
221
222
223 REAL RWORK( * ), S( * )
224 COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
225 $ WORK( * )
226
227
228
229
230
231 COMPLEX CZERO, CONE
232 parameter( czero = ( 0.0e0, 0.0e0 ),
233 $ cone = ( 1.0e0, 0.0e0 ) )
234 REAL ZERO, ONE
235 parameter( zero = 0.0e0, one = 1.0e0 )
236
237
238 LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
239 $ WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS
240 INTEGER BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL,
241 $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
242 $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
243 $ NRVT, WRKBL
244 INTEGER LWORK_CGEQRF, LWORK_CUNGQR_N, LWORK_CUNGQR_M,
245 $ LWORK_CGEBRD, LWORK_CUNGBR_P, LWORK_CUNGBR_Q,
246 $ LWORK_CGELQF, LWORK_CUNGLQ_N, LWORK_CUNGLQ_M
247 REAL ANRM, BIGNUM, EPS, SMLNUM
248
249
250 REAL DUM( 1 )
251 COMPLEX CDUM( 1 )
252
253
258
259
260 LOGICAL LSAME
261 INTEGER ILAENV
262 REAL CLANGE, SLAMCH, SROUNDUP_LWORK
265
266
267 INTRINSIC max, min, sqrt
268
269
270
271
272
273 info = 0
274 minmn = min( m, n )
275 wntua =
lsame( jobu,
'A' )
276 wntus =
lsame( jobu,
'S' )
277 wntuas = wntua .OR. wntus
278 wntuo =
lsame( jobu,
'O' )
279 wntun =
lsame( jobu,
'N' )
280 wntva =
lsame( jobvt,
'A' )
281 wntvs =
lsame( jobvt,
'S' )
282 wntvas = wntva .OR. wntvs
283 wntvo =
lsame( jobvt,
'O' )
284 wntvn =
lsame( jobvt,
'N' )
285 lquery = ( lwork.EQ.-1 )
286
287 IF( .NOT.( wntua .OR. wntus .OR. wntuo .OR. wntun ) ) THEN
288 info = -1
289 ELSE IF( .NOT.( wntva .OR. wntvs .OR. wntvo .OR. wntvn ) .OR.
290 $ ( wntvo .AND. wntuo ) ) THEN
291 info = -2
292 ELSE IF( m.LT.0 ) THEN
293 info = -3
294 ELSE IF( n.LT.0 ) THEN
295 info = -4
296 ELSE IF( lda.LT.max( 1, m ) ) THEN
297 info = -6
298 ELSE IF( ldu.LT.1 .OR. ( wntuas .AND. ldu.LT.m ) ) THEN
299 info = -9
300 ELSE IF( ldvt.LT.1 .OR. ( wntva .AND. ldvt.LT.n ) .OR.
301 $ ( wntvs .AND. ldvt.LT.minmn ) ) THEN
302 info = -11
303 END IF
304
305
306
307
308
309
310
311
312
313 IF( info.EQ.0 ) THEN
314 minwrk = 1
315 maxwrk = 1
316 IF( m.GE.n .AND. minmn.GT.0 ) THEN
317
318
319
320 mnthr =
ilaenv( 6,
'CGESVD', jobu // jobvt, m, n, 0, 0 )
321
322 CALL cgeqrf( m, n, a, lda, cdum(1), cdum(1), -1, ierr )
323 lwork_cgeqrf = int( cdum(1) )
324
325 CALL cungqr( m, n, n, a, lda, cdum(1), cdum(1), -1,
326 $ ierr )
327 lwork_cungqr_n = int( cdum(1) )
328 CALL cungqr( m, m, n, a, lda, cdum(1), cdum(1), -1,
329 $ ierr )
330 lwork_cungqr_m = int( cdum(1) )
331
332 CALL cgebrd( n, n, a, lda, s, dum(1), cdum(1),
333 $ cdum(1), cdum(1), -1, ierr )
334 lwork_cgebrd = int( cdum(1) )
335
336 CALL cungbr(
'P', n, n, n, a, lda, cdum(1),
337 $ cdum(1), -1, ierr )
338 lwork_cungbr_p = int( cdum(1) )
339 CALL cungbr(
'Q', n, n, n, a, lda, cdum(1),
340 $ cdum(1), -1, ierr )
341 lwork_cungbr_q = int( cdum(1) )
342
343 mnthr =
ilaenv( 6,
'CGESVD', jobu // jobvt, m, n, 0, 0 )
344 IF( m.GE.mnthr ) THEN
345 IF( wntun ) THEN
346
347
348
349 maxwrk = n + lwork_cgeqrf
350 maxwrk = max( maxwrk, 2*n+lwork_cgebrd )
351 IF( wntvo .OR. wntvas )
352 $ maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
353 minwrk = 3*n
354 ELSE IF( wntuo .AND. wntvn ) THEN
355
356
357
358 wrkbl = n + lwork_cgeqrf
359 wrkbl = max( wrkbl, n+lwork_cungqr_n )
360 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
361 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
362 maxwrk = max( n*n+wrkbl, n*n+m*n )
363 minwrk = 2*n + m
364 ELSE IF( wntuo .AND. wntvas ) THEN
365
366
367
368
369 wrkbl = n + lwork_cgeqrf
370 wrkbl = max( wrkbl, n+lwork_cungqr_n )
371 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
372 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
373 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
374 maxwrk = max( n*n+wrkbl, n*n+m*n )
375 minwrk = 2*n + m
376 ELSE IF( wntus .AND. wntvn ) THEN
377
378
379
380 wrkbl = n + lwork_cgeqrf
381 wrkbl = max( wrkbl, n+lwork_cungqr_n )
382 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
383 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
384 maxwrk = n*n + wrkbl
385 minwrk = 2*n + m
386 ELSE IF( wntus .AND. wntvo ) THEN
387
388
389
390 wrkbl = n + lwork_cgeqrf
391 wrkbl = max( wrkbl, n+lwork_cungqr_n )
392 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
393 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
394 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
395 maxwrk = 2*n*n + wrkbl
396 minwrk = 2*n + m
397 ELSE IF( wntus .AND. wntvas ) THEN
398
399
400
401
402 wrkbl = n + lwork_cgeqrf
403 wrkbl = max( wrkbl, n+lwork_cungqr_n )
404 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
405 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
406 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
407 maxwrk = n*n + wrkbl
408 minwrk = 2*n + m
409 ELSE IF( wntua .AND. wntvn ) THEN
410
411
412
413 wrkbl = n + lwork_cgeqrf
414 wrkbl = max( wrkbl, n+lwork_cungqr_m )
415 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
416 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
417 maxwrk = n*n + wrkbl
418 minwrk = 2*n + m
419 ELSE IF( wntua .AND. wntvo ) THEN
420
421
422
423 wrkbl = n + lwork_cgeqrf
424 wrkbl = max( wrkbl, n+lwork_cungqr_m )
425 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
426 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
427 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
428 maxwrk = 2*n*n + wrkbl
429 minwrk = 2*n + m
430 ELSE IF( wntua .AND. wntvas ) THEN
431
432
433
434
435 wrkbl = n + lwork_cgeqrf
436 wrkbl = max( wrkbl, n+lwork_cungqr_m )
437 wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
438 wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
439 wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
440 maxwrk = n*n + wrkbl
441 minwrk = 2*n + m
442 END IF
443 ELSE
444
445
446
447 CALL cgebrd( m, n, a, lda, s, dum(1), cdum(1),
448 $ cdum(1), cdum(1), -1, ierr )
449 lwork_cgebrd = int( cdum(1) )
450 maxwrk = 2*n + lwork_cgebrd
451 IF( wntus .OR. wntuo ) THEN
452 CALL cungbr(
'Q', m, n, n, a, lda, cdum(1),
453 $ cdum(1), -1, ierr )
454 lwork_cungbr_q = int( cdum(1) )
455 maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
456 END IF
457 IF( wntua ) THEN
458 CALL cungbr(
'Q', m, m, n, a, lda, cdum(1),
459 $ cdum(1), -1, ierr )
460 lwork_cungbr_q = int( cdum(1) )
461 maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
462 END IF
463 IF( .NOT.wntvn ) THEN
464 maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
465 END IF
466 minwrk = 2*n + m
467 END IF
468 ELSE IF( minmn.GT.0 ) THEN
469
470
471
472 mnthr =
ilaenv( 6,
'CGESVD', jobu // jobvt, m, n, 0, 0 )
473
474 CALL cgelqf( m, n, a, lda, cdum(1), cdum(1), -1, ierr )
475 lwork_cgelqf = int( cdum(1) )
476
477 CALL cunglq( n, n, m, cdum(1), n, cdum(1), cdum(1), -1,
478 $ ierr )
479 lwork_cunglq_n = int( cdum(1) )
480 CALL cunglq( m, n, m, a, lda, cdum(1), cdum(1), -1,
481 $ ierr )
482 lwork_cunglq_m = int( cdum(1) )
483
484 CALL cgebrd( m, m, a, lda, s, dum(1), cdum(1),
485 $ cdum(1), cdum(1), -1, ierr )
486 lwork_cgebrd = int( cdum(1) )
487
488 CALL cungbr(
'P', m, m, m, a, n, cdum(1),
489 $ cdum(1), -1, ierr )
490 lwork_cungbr_p = int( cdum(1) )
491
492 CALL cungbr(
'Q', m, m, m, a, n, cdum(1),
493 $ cdum(1), -1, ierr )
494 lwork_cungbr_q = int( cdum(1) )
495 IF( n.GE.mnthr ) THEN
496 IF( wntvn ) THEN
497
498
499
500 maxwrk = m + lwork_cgelqf
501 maxwrk = max( maxwrk, 2*m+lwork_cgebrd )
502 IF( wntuo .OR. wntuas )
503 $ maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
504 minwrk = 3*m
505 ELSE IF( wntvo .AND. wntun ) THEN
506
507
508
509 wrkbl = m + lwork_cgelqf
510 wrkbl = max( wrkbl, m+lwork_cunglq_m )
511 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
512 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
513 maxwrk = max( m*m+wrkbl, m*m+m*n )
514 minwrk = 2*m + n
515 ELSE IF( wntvo .AND. wntuas ) THEN
516
517
518
519
520 wrkbl = m + lwork_cgelqf
521 wrkbl = max( wrkbl, m+lwork_cunglq_m )
522 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
523 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
524 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
525 maxwrk = max( m*m+wrkbl, m*m+m*n )
526 minwrk = 2*m + n
527 ELSE IF( wntvs .AND. wntun ) THEN
528
529
530
531 wrkbl = m + lwork_cgelqf
532 wrkbl = max( wrkbl, m+lwork_cunglq_m )
533 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
534 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
535 maxwrk = m*m + wrkbl
536 minwrk = 2*m + n
537 ELSE IF( wntvs .AND. wntuo ) THEN
538
539
540
541 wrkbl = m + lwork_cgelqf
542 wrkbl = max( wrkbl, m+lwork_cunglq_m )
543 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
544 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
545 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
546 maxwrk = 2*m*m + wrkbl
547 minwrk = 2*m + n
548 ELSE IF( wntvs .AND. wntuas ) THEN
549
550
551
552
553 wrkbl = m + lwork_cgelqf
554 wrkbl = max( wrkbl, m+lwork_cunglq_m )
555 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
556 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
557 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
558 maxwrk = m*m + wrkbl
559 minwrk = 2*m + n
560 ELSE IF( wntva .AND. wntun ) THEN
561
562
563
564 wrkbl = m + lwork_cgelqf
565 wrkbl = max( wrkbl, m+lwork_cunglq_n )
566 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
567 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
568 maxwrk = m*m + wrkbl
569 minwrk = 2*m + n
570 ELSE IF( wntva .AND. wntuo ) THEN
571
572
573
574 wrkbl = m + lwork_cgelqf
575 wrkbl = max( wrkbl, m+lwork_cunglq_n )
576 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
577 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
578 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
579 maxwrk = 2*m*m + wrkbl
580 minwrk = 2*m + n
581 ELSE IF( wntva .AND. wntuas ) THEN
582
583
584
585
586 wrkbl = m + lwork_cgelqf
587 wrkbl = max( wrkbl, m+lwork_cunglq_n )
588 wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
589 wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
590 wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
591 maxwrk = m*m + wrkbl
592 minwrk = 2*m + n
593 END IF
594 ELSE
595
596
597
598 CALL cgebrd( m, n, a, lda, s, dum(1), cdum(1),
599 $ cdum(1), cdum(1), -1, ierr )
600 lwork_cgebrd = int( cdum(1) )
601 maxwrk = 2*m + lwork_cgebrd
602 IF( wntvs .OR. wntvo ) THEN
603
604 CALL cungbr(
'P', m, n, m, a, n, cdum(1),
605 $ cdum(1), -1, ierr )
606 lwork_cungbr_p = int( cdum(1) )
607 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
608 END IF
609 IF( wntva ) THEN
610 CALL cungbr(
'P', n, n, m, a, n, cdum(1),
611 $ cdum(1), -1, ierr )
612 lwork_cungbr_p = int( cdum(1) )
613 maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
614 END IF
615 IF( .NOT.wntun ) THEN
616 maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
617 END IF
618 minwrk = 2*m + n
619 END IF
620 END IF
621 maxwrk = max( minwrk, maxwrk )
623
624 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
625 info = -13
626 END IF
627 END IF
628
629 IF( info.NE.0 ) THEN
630 CALL xerbla(
'CGESVD', -info )
631 RETURN
632 ELSE IF( lquery ) THEN
633 RETURN
634 END IF
635
636
637
638 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
639 RETURN
640 END IF
641
642
643
645 smlnum = sqrt(
slamch(
'S' ) ) / eps
646 bignum = one / smlnum
647
648
649
650 anrm =
clange(
'M', m, n, a, lda, dum )
651 iscl = 0
652 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
653 iscl = 1
654 CALL clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
655 ELSE IF( anrm.GT.bignum ) THEN
656 iscl = 1
657 CALL clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
658 END IF
659
660 IF( m.GE.n ) THEN
661
662
663
664
665
666 IF( m.GE.mnthr ) THEN
667
668 IF( wntun ) THEN
669
670
671
672
673 itau = 1
674 iwork = itau + n
675
676
677
678
679
680 CALL cgeqrf( m, n, a, lda, work( itau ),
681 $ work( iwork ),
682 $ lwork-iwork+1, ierr )
683
684
685
686 IF( n .GT. 1 ) THEN
687 CALL claset(
'L', n-1, n-1, czero, czero, a( 2,
688 $ 1 ),
689 $ lda )
690 END IF
691 ie = 1
692 itauq = 1
693 itaup = itauq + n
694 iwork = itaup + n
695
696
697
698
699
700 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
701 $ work( itauq ),
702 $ work( itaup ), work( iwork ), lwork-iwork+1,
703 $ ierr )
704 ncvt = 0
705 IF( wntvo .OR. wntvas ) THEN
706
707
708
709
710
711 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
712 $ work( iwork ), lwork-iwork+1, ierr )
713 ncvt = n
714 END IF
715 irwork = ie + n
716
717
718
719
720
721
722 CALL cbdsqr(
'U', n, ncvt, 0, 0, s, rwork( ie ), a,
723 $ lda,
724 $ cdum, 1, cdum, 1, rwork( irwork ), info )
725
726
727
728 IF( wntvas )
729 $
CALL clacpy(
'F', n, n, a, lda, vt, ldvt )
730
731 ELSE IF( wntuo .AND. wntvn ) THEN
732
733
734
735
736
737 IF( lwork.GE.n*n+3*n ) THEN
738
739
740
741 ir = 1
742 IF( lwork.GE.max( wrkbl, lda*n )+lda*n ) THEN
743
744
745
746 ldwrku = lda
747 ldwrkr = lda
748 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n ) THEN
749
750
751
752 ldwrku = lda
753 ldwrkr = n
754 ELSE
755
756
757
758 ldwrku = ( lwork-n*n ) / n
759 ldwrkr = n
760 END IF
761 itau = ir + ldwrkr*n
762 iwork = itau + n
763
764
765
766
767
768 CALL cgeqrf( m, n, a, lda, work( itau ),
769 $ work( iwork ), lwork-iwork+1, ierr )
770
771
772
773 CALL clacpy(
'U', n, n, a, lda, work( ir ),
774 $ ldwrkr )
775 CALL claset(
'L', n-1, n-1, czero, czero,
776 $ work( ir+1 ), ldwrkr )
777
778
779
780
781
782 CALL cungqr( m, n, n, a, lda, work( itau ),
783 $ work( iwork ), lwork-iwork+1, ierr )
784 ie = 1
785 itauq = itau
786 itaup = itauq + n
787 iwork = itaup + n
788
789
790
791
792
793 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
794 $ rwork( ie ),
795 $ work( itauq ), work( itaup ),
796 $ work( iwork ), lwork-iwork+1, ierr )
797
798
799
800
801
802 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
803 $ work( itauq ), work( iwork ),
804 $ lwork-iwork+1, ierr )
805 irwork = ie + n
806
807
808
809
810
811
812 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ), cdum,
813 $ 1,
814 $ work( ir ), ldwrkr, cdum, 1,
815 $ rwork( irwork ), info )
816 iu = itauq
817
818
819
820
821
822
823 DO 10 i = 1, m, ldwrku
824 chunk = min( m-i+1, ldwrku )
825 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i,
826 $ 1 ),
827 $ lda, work( ir ), ldwrkr, czero,
828 $ work( iu ), ldwrku )
829 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
830 $ a( i, 1 ), lda )
831 10 CONTINUE
832
833 ELSE
834
835
836
837 ie = 1
838 itauq = 1
839 itaup = itauq + n
840 iwork = itaup + n
841
842
843
844
845
846 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
847 $ work( itauq ), work( itaup ),
848 $ work( iwork ), lwork-iwork+1, ierr )
849
850
851
852
853
854 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
855 $ work( iwork ), lwork-iwork+1, ierr )
856 irwork = ie + n
857
858
859
860
861
862
863 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ), cdum,
864 $ 1,
865 $ a, lda, cdum, 1, rwork( irwork ), info )
866
867 END IF
868
869 ELSE IF( wntuo .AND. wntvas ) THEN
870
871
872
873
874
875 IF( lwork.GE.n*n+3*n ) THEN
876
877
878
879 ir = 1
880 IF( lwork.GE.max( wrkbl, lda*n )+lda*n ) THEN
881
882
883
884 ldwrku = lda
885 ldwrkr = lda
886 ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n ) THEN
887
888
889
890 ldwrku = lda
891 ldwrkr = n
892 ELSE
893
894
895
896 ldwrku = ( lwork-n*n ) / n
897 ldwrkr = n
898 END IF
899 itau = ir + ldwrkr*n
900 iwork = itau + n
901
902
903
904
905
906 CALL cgeqrf( m, n, a, lda, work( itau ),
907 $ work( iwork ), lwork-iwork+1, ierr )
908
909
910
911 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
912 IF( n.GT.1 )
913 $
CALL claset(
'L', n-1, n-1, czero, czero,
914 $ vt( 2, 1 ), ldvt )
915
916
917
918
919
920 CALL cungqr( m, n, n, a, lda, work( itau ),
921 $ work( iwork ), lwork-iwork+1, ierr )
922 ie = 1
923 itauq = itau
924 itaup = itauq + n
925 iwork = itaup + n
926
927
928
929
930
931 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
932 $ work( itauq ), work( itaup ),
933 $ work( iwork ), lwork-iwork+1, ierr )
934 CALL clacpy(
'L', n, n, vt, ldvt, work( ir ),
935 $ ldwrkr )
936
937
938
939
940
941 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
942 $ work( itauq ), work( iwork ),
943 $ lwork-iwork+1, ierr )
944
945
946
947
948
949 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
950 $ work( iwork ), lwork-iwork+1, ierr )
951 irwork = ie + n
952
953
954
955
956
957
958
959 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ), vt,
960 $ ldvt, work( ir ), ldwrkr, cdum, 1,
961 $ rwork( irwork ), info )
962 iu = itauq
963
964
965
966
967
968
969 DO 20 i = 1, m, ldwrku
970 chunk = min( m-i+1, ldwrku )
971 CALL cgemm(
'N',
'N', chunk, n, n, cone, a( i,
972 $ 1 ),
973 $ lda, work( ir ), ldwrkr, czero,
974 $ work( iu ), ldwrku )
975 CALL clacpy(
'F', chunk, n, work( iu ), ldwrku,
976 $ a( i, 1 ), lda )
977 20 CONTINUE
978
979 ELSE
980
981
982
983 itau = 1
984 iwork = itau + n
985
986
987
988
989
990 CALL cgeqrf( m, n, a, lda, work( itau ),
991 $ work( iwork ), lwork-iwork+1, ierr )
992
993
994
995 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
996 IF( n.GT.1 )
997 $
CALL claset(
'L', n-1, n-1, czero, czero,
998 $ vt( 2, 1 ), ldvt )
999
1000
1001
1002
1003
1004 CALL cungqr( m, n, n, a, lda, work( itau ),
1005 $ work( iwork ), lwork-iwork+1, ierr )
1006 ie = 1
1007 itauq = itau
1008 itaup = itauq + n
1009 iwork = itaup + n
1010
1011
1012
1013
1014
1015 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1016 $ work( itauq ), work( itaup ),
1017 $ work( iwork ), lwork-iwork+1, ierr )
1018
1019
1020
1021
1022
1023 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1024 $ work( itauq ), a, lda, work( iwork ),
1025 $ lwork-iwork+1, ierr )
1026
1027
1028
1029
1030
1031 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
1032 $ work( iwork ), lwork-iwork+1, ierr )
1033 irwork = ie + n
1034
1035
1036
1037
1038
1039
1040
1041 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), vt,
1042 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1043 $ info )
1044
1045 END IF
1046
1047 ELSE IF( wntus ) THEN
1048
1049 IF( wntvn ) THEN
1050
1051
1052
1053
1054
1055 IF( lwork.GE.n*n+3*n ) THEN
1056
1057
1058
1059 ir = 1
1060 IF( lwork.GE.wrkbl+lda*n ) THEN
1061
1062
1063
1064 ldwrkr = lda
1065 ELSE
1066
1067
1068
1069 ldwrkr = n
1070 END IF
1071 itau = ir + ldwrkr*n
1072 iwork = itau + n
1073
1074
1075
1076
1077
1078 CALL cgeqrf( m, n, a, lda, work( itau ),
1079 $ work( iwork ), lwork-iwork+1, ierr )
1080
1081
1082
1083 CALL clacpy(
'U', n, n, a, lda, work( ir ),
1084 $ ldwrkr )
1085 CALL claset(
'L', n-1, n-1, czero, czero,
1086 $ work( ir+1 ), ldwrkr )
1087
1088
1089
1090
1091
1092 CALL cungqr( m, n, n, a, lda, work( itau ),
1093 $ work( iwork ), lwork-iwork+1, ierr )
1094 ie = 1
1095 itauq = itau
1096 itaup = itauq + n
1097 iwork = itaup + n
1098
1099
1100
1101
1102
1103 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1104 $ rwork( ie ), work( itauq ),
1105 $ work( itaup ), work( iwork ),
1106 $ lwork-iwork+1, ierr )
1107
1108
1109
1110
1111
1112 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1113 $ work( itauq ), work( iwork ),
1114 $ lwork-iwork+1, ierr )
1115 irwork = ie + n
1116
1117
1118
1119
1120
1121
1122 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ),
1123 $ cdum,
1124 $ 1, work( ir ), ldwrkr, cdum, 1,
1125 $ rwork( irwork ), info )
1126
1127
1128
1129
1130
1131
1132 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1133 $ work( ir ), ldwrkr, czero, u, ldu )
1134
1135 ELSE
1136
1137
1138
1139 itau = 1
1140 iwork = itau + n
1141
1142
1143
1144
1145
1146 CALL cgeqrf( m, n, a, lda, work( itau ),
1147 $ work( iwork ), lwork-iwork+1, ierr )
1148 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1149
1150
1151
1152
1153
1154 CALL cungqr( m, n, n, u, ldu, work( itau ),
1155 $ work( iwork ), lwork-iwork+1, ierr )
1156 ie = 1
1157 itauq = itau
1158 itaup = itauq + n
1159 iwork = itaup + n
1160
1161
1162
1163 IF( n .GT. 1 ) THEN
1164 CALL claset(
'L', n-1, n-1, czero, czero,
1165 $ a( 2, 1 ), lda )
1166 END IF
1167
1168
1169
1170
1171
1172 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1173 $ work( itauq ), work( itaup ),
1174 $ work( iwork ), lwork-iwork+1, ierr )
1175
1176
1177
1178
1179
1180 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1181 $ work( itauq ), u, ldu, work( iwork ),
1182 $ lwork-iwork+1, ierr )
1183 irwork = ie + n
1184
1185
1186
1187
1188
1189
1190 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ),
1191 $ cdum,
1192 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1193 $ info )
1194
1195 END IF
1196
1197 ELSE IF( wntvo ) THEN
1198
1199
1200
1201
1202
1203 IF( lwork.GE.2*n*n+3*n ) THEN
1204
1205
1206
1207 iu = 1
1208 IF( lwork.GE.wrkbl+2*lda*n ) THEN
1209
1210
1211
1212 ldwrku = lda
1213 ir = iu + ldwrku*n
1214 ldwrkr = lda
1215 ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1216
1217
1218
1219 ldwrku = lda
1220 ir = iu + ldwrku*n
1221 ldwrkr = n
1222 ELSE
1223
1224
1225
1226 ldwrku = n
1227 ir = iu + ldwrku*n
1228 ldwrkr = n
1229 END IF
1230 itau = ir + ldwrkr*n
1231 iwork = itau + n
1232
1233
1234
1235
1236
1237 CALL cgeqrf( m, n, a, lda, work( itau ),
1238 $ work( iwork ), lwork-iwork+1, ierr )
1239
1240
1241
1242 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1243 $ ldwrku )
1244 CALL claset(
'L', n-1, n-1, czero, czero,
1245 $ work( iu+1 ), ldwrku )
1246
1247
1248
1249
1250
1251 CALL cungqr( m, n, n, a, lda, work( itau ),
1252 $ work( iwork ), lwork-iwork+1, ierr )
1253 ie = 1
1254 itauq = itau
1255 itaup = itauq + n
1256 iwork = itaup + n
1257
1258
1259
1260
1261
1262
1263
1264 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1265 $ rwork( ie ), work( itauq ),
1266 $ work( itaup ), work( iwork ),
1267 $ lwork-iwork+1, ierr )
1268 CALL clacpy(
'U', n, n, work( iu ), ldwrku,
1269 $ work( ir ), ldwrkr )
1270
1271
1272
1273
1274
1275 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1276 $ work( itauq ), work( iwork ),
1277 $ lwork-iwork+1, ierr )
1278
1279
1280
1281
1282
1283
1284 CALL cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1285 $ work( itaup ), work( iwork ),
1286 $ lwork-iwork+1, ierr )
1287 irwork = ie + n
1288
1289
1290
1291
1292
1293
1294
1295 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1296 $ work( ir ), ldwrkr, work( iu ),
1297 $ ldwrku, cdum, 1, rwork( irwork ),
1298 $ info )
1299
1300
1301
1302
1303
1304
1305 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1306 $ work( iu ), ldwrku, czero, u, ldu )
1307
1308
1309
1310
1311
1312 CALL clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1313 $ lda )
1314
1315 ELSE
1316
1317
1318
1319 itau = 1
1320 iwork = itau + n
1321
1322
1323
1324
1325
1326 CALL cgeqrf( m, n, a, lda, work( itau ),
1327 $ work( iwork ), lwork-iwork+1, ierr )
1328 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1329
1330
1331
1332
1333
1334 CALL cungqr( m, n, n, u, ldu, work( itau ),
1335 $ work( iwork ), lwork-iwork+1, ierr )
1336 ie = 1
1337 itauq = itau
1338 itaup = itauq + n
1339 iwork = itaup + n
1340
1341
1342
1343 IF( n .GT. 1 ) THEN
1344 CALL claset(
'L', n-1, n-1, czero, czero,
1345 $ a( 2, 1 ), lda )
1346 END IF
1347
1348
1349
1350
1351
1352 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1353 $ work( itauq ), work( itaup ),
1354 $ work( iwork ), lwork-iwork+1, ierr )
1355
1356
1357
1358
1359
1360 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1361 $ work( itauq ), u, ldu, work( iwork ),
1362 $ lwork-iwork+1, ierr )
1363
1364
1365
1366
1367
1368 CALL cungbr(
'P', n, n, n, a, lda,
1369 $ work( itaup ),
1370 $ work( iwork ), lwork-iwork+1, ierr )
1371 irwork = ie + n
1372
1373
1374
1375
1376
1377
1378
1379 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1380 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1381 $ info )
1382
1383 END IF
1384
1385 ELSE IF( wntvas ) THEN
1386
1387
1388
1389
1390
1391
1392 IF( lwork.GE.n*n+3*n ) THEN
1393
1394
1395
1396 iu = 1
1397 IF( lwork.GE.wrkbl+lda*n ) THEN
1398
1399
1400
1401 ldwrku = lda
1402 ELSE
1403
1404
1405
1406 ldwrku = n
1407 END IF
1408 itau = iu + ldwrku*n
1409 iwork = itau + n
1410
1411
1412
1413
1414
1415 CALL cgeqrf( m, n, a, lda, work( itau ),
1416 $ work( iwork ), lwork-iwork+1, ierr )
1417
1418
1419
1420 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1421 $ ldwrku )
1422 CALL claset(
'L', n-1, n-1, czero, czero,
1423 $ work( iu+1 ), ldwrku )
1424
1425
1426
1427
1428
1429 CALL cungqr( m, n, n, a, lda, work( itau ),
1430 $ work( iwork ), lwork-iwork+1, ierr )
1431 ie = 1
1432 itauq = itau
1433 itaup = itauq + n
1434 iwork = itaup + n
1435
1436
1437
1438
1439
1440 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1441 $ rwork( ie ), work( itauq ),
1442 $ work( itaup ), work( iwork ),
1443 $ lwork-iwork+1, ierr )
1444 CALL clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1445 $ ldvt )
1446
1447
1448
1449
1450
1451 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1452 $ work( itauq ), work( iwork ),
1453 $ lwork-iwork+1, ierr )
1454
1455
1456
1457
1458
1459
1460 CALL cungbr(
'P', n, n, n, vt, ldvt,
1461 $ work( itaup ),
1462 $ work( iwork ), lwork-iwork+1, ierr )
1463 irwork = ie + n
1464
1465
1466
1467
1468
1469
1470
1471 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1472 $ vt,
1473 $ ldvt, work( iu ), ldwrku, cdum, 1,
1474 $ rwork( irwork ), info )
1475
1476
1477
1478
1479
1480
1481 CALL cgemm(
'N',
'N', m, n, n, cone, a, lda,
1482 $ work( iu ), ldwrku, czero, u, ldu )
1483
1484 ELSE
1485
1486
1487
1488 itau = 1
1489 iwork = itau + n
1490
1491
1492
1493
1494
1495 CALL cgeqrf( m, n, a, lda, work( itau ),
1496 $ work( iwork ), lwork-iwork+1, ierr )
1497 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1498
1499
1500
1501
1502
1503 CALL cungqr( m, n, n, u, ldu, work( itau ),
1504 $ work( iwork ), lwork-iwork+1, ierr )
1505
1506
1507
1508 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
1509 IF( n.GT.1 )
1510 $
CALL claset(
'L', n-1, n-1, czero, czero,
1511 $ vt( 2, 1 ), ldvt )
1512 ie = 1
1513 itauq = itau
1514 itaup = itauq + n
1515 iwork = itaup + n
1516
1517
1518
1519
1520
1521 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1522 $ work( itauq ), work( itaup ),
1523 $ work( iwork ), lwork-iwork+1, ierr )
1524
1525
1526
1527
1528
1529
1530 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
1531 $ work( itauq ), u, ldu, work( iwork ),
1532 $ lwork-iwork+1, ierr )
1533
1534
1535
1536
1537
1538 CALL cungbr(
'P', n, n, n, vt, ldvt,
1539 $ work( itaup ),
1540 $ work( iwork ), lwork-iwork+1, ierr )
1541 irwork = ie + n
1542
1543
1544
1545
1546
1547
1548
1549 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ),
1550 $ vt,
1551 $ ldvt, u, ldu, cdum, 1,
1552 $ rwork( irwork ), info )
1553
1554 END IF
1555
1556 END IF
1557
1558 ELSE IF( wntua ) THEN
1559
1560 IF( wntvn ) THEN
1561
1562
1563
1564
1565
1566 IF( lwork.GE.n*n+max( n+m, 3*n ) ) THEN
1567
1568
1569
1570 ir = 1
1571 IF( lwork.GE.wrkbl+lda*n ) THEN
1572
1573
1574
1575 ldwrkr = lda
1576 ELSE
1577
1578
1579
1580 ldwrkr = n
1581 END IF
1582 itau = ir + ldwrkr*n
1583 iwork = itau + n
1584
1585
1586
1587
1588
1589 CALL cgeqrf( m, n, a, lda, work( itau ),
1590 $ work( iwork ), lwork-iwork+1, ierr )
1591 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1592
1593
1594
1595 CALL clacpy(
'U', n, n, a, lda, work( ir ),
1596 $ ldwrkr )
1597 CALL claset(
'L', n-1, n-1, czero, czero,
1598 $ work( ir+1 ), ldwrkr )
1599
1600
1601
1602
1603
1604 CALL cungqr( m, m, n, u, ldu, work( itau ),
1605 $ work( iwork ), lwork-iwork+1, ierr )
1606 ie = 1
1607 itauq = itau
1608 itaup = itauq + n
1609 iwork = itaup + n
1610
1611
1612
1613
1614
1615 CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1616 $ rwork( ie ), work( itauq ),
1617 $ work( itaup ), work( iwork ),
1618 $ lwork-iwork+1, ierr )
1619
1620
1621
1622
1623
1624 CALL cungbr(
'Q', n, n, n, work( ir ), ldwrkr,
1625 $ work( itauq ), work( iwork ),
1626 $ lwork-iwork+1, ierr )
1627 irwork = ie + n
1628
1629
1630
1631
1632
1633
1634 CALL cbdsqr(
'U', n, 0, n, 0, s, rwork( ie ),
1635 $ cdum,
1636 $ 1, work( ir ), ldwrkr, cdum, 1,
1637 $ rwork( irwork ), info )
1638
1639
1640
1641
1642
1643
1644 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1645 $ work( ir ), ldwrkr, czero, a, lda )
1646
1647
1648
1649 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1650
1651 ELSE
1652
1653
1654
1655 itau = 1
1656 iwork = itau + n
1657
1658
1659
1660
1661
1662 CALL cgeqrf( m, n, a, lda, work( itau ),
1663 $ work( iwork ), lwork-iwork+1, ierr )
1664 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1665
1666
1667
1668
1669
1670 CALL cungqr( m, m, n, u, ldu, work( itau ),
1671 $ work( iwork ), lwork-iwork+1, ierr )
1672 ie = 1
1673 itauq = itau
1674 itaup = itauq + n
1675 iwork = itaup + n
1676
1677
1678
1679 IF( n .GT. 1 ) THEN
1680 CALL claset(
'L', n-1, n-1, czero, czero,
1681 $ a( 2, 1 ), lda )
1682 END IF
1683
1684
1685
1686
1687
1688 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1689 $ work( itauq ), work( itaup ),
1690 $ work( iwork ), lwork-iwork+1, ierr )
1691
1692
1693
1694
1695
1696
1697 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1698 $ work( itauq ), u, ldu, work( iwork ),
1699 $ lwork-iwork+1, ierr )
1700 irwork = ie + n
1701
1702
1703
1704
1705
1706
1707 CALL cbdsqr(
'U', n, 0, m, 0, s, rwork( ie ),
1708 $ cdum,
1709 $ 1, u, ldu, cdum, 1, rwork( irwork ),
1710 $ info )
1711
1712 END IF
1713
1714 ELSE IF( wntvo ) THEN
1715
1716
1717
1718
1719
1720 IF( lwork.GE.2*n*n+max( n+m, 3*n ) ) THEN
1721
1722
1723
1724 iu = 1
1725 IF( lwork.GE.wrkbl+2*lda*n ) THEN
1726
1727
1728
1729 ldwrku = lda
1730 ir = iu + ldwrku*n
1731 ldwrkr = lda
1732 ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1733
1734
1735
1736 ldwrku = lda
1737 ir = iu + ldwrku*n
1738 ldwrkr = n
1739 ELSE
1740
1741
1742
1743 ldwrku = n
1744 ir = iu + ldwrku*n
1745 ldwrkr = n
1746 END IF
1747 itau = ir + ldwrkr*n
1748 iwork = itau + n
1749
1750
1751
1752
1753
1754 CALL cgeqrf( m, n, a, lda, work( itau ),
1755 $ work( iwork ), lwork-iwork+1, ierr )
1756 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1757
1758
1759
1760
1761
1762 CALL cungqr( m, m, n, u, ldu, work( itau ),
1763 $ work( iwork ), lwork-iwork+1, ierr )
1764
1765
1766
1767 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1768 $ ldwrku )
1769 CALL claset(
'L', n-1, n-1, czero, czero,
1770 $ work( iu+1 ), ldwrku )
1771 ie = 1
1772 itauq = itau
1773 itaup = itauq + n
1774 iwork = itaup + n
1775
1776
1777
1778
1779
1780
1781
1782 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1783 $ rwork( ie ), work( itauq ),
1784 $ work( itaup ), work( iwork ),
1785 $ lwork-iwork+1, ierr )
1786 CALL clacpy(
'U', n, n, work( iu ), ldwrku,
1787 $ work( ir ), ldwrkr )
1788
1789
1790
1791
1792
1793 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1794 $ work( itauq ), work( iwork ),
1795 $ lwork-iwork+1, ierr )
1796
1797
1798
1799
1800
1801
1802 CALL cungbr(
'P', n, n, n, work( ir ), ldwrkr,
1803 $ work( itaup ), work( iwork ),
1804 $ lwork-iwork+1, ierr )
1805 irwork = ie + n
1806
1807
1808
1809
1810
1811
1812
1813 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1814 $ work( ir ), ldwrkr, work( iu ),
1815 $ ldwrku, cdum, 1, rwork( irwork ),
1816 $ info )
1817
1818
1819
1820
1821
1822
1823 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
1824 $ work( iu ), ldwrku, czero, a, lda )
1825
1826
1827
1828 CALL clacpy(
'F', m, n, a, lda, u, ldu )
1829
1830
1831
1832 CALL clacpy(
'F', n, n, work( ir ), ldwrkr, a,
1833 $ lda )
1834
1835 ELSE
1836
1837
1838
1839 itau = 1
1840 iwork = itau + n
1841
1842
1843
1844
1845
1846 CALL cgeqrf( m, n, a, lda, work( itau ),
1847 $ work( iwork ), lwork-iwork+1, ierr )
1848 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1849
1850
1851
1852
1853
1854 CALL cungqr( m, m, n, u, ldu, work( itau ),
1855 $ work( iwork ), lwork-iwork+1, ierr )
1856 ie = 1
1857 itauq = itau
1858 itaup = itauq + n
1859 iwork = itaup + n
1860
1861
1862
1863 IF( n .GT. 1 ) THEN
1864 CALL claset(
'L', n-1, n-1, czero, czero,
1865 $ a( 2, 1 ), lda )
1866 END IF
1867
1868
1869
1870
1871
1872 CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1873 $ work( itauq ), work( itaup ),
1874 $ work( iwork ), lwork-iwork+1, ierr )
1875
1876
1877
1878
1879
1880
1881 CALL cunmbr(
'Q',
'R',
'N', m, n, n, a, lda,
1882 $ work( itauq ), u, ldu, work( iwork ),
1883 $ lwork-iwork+1, ierr )
1884
1885
1886
1887
1888
1889 CALL cungbr(
'P', n, n, n, a, lda,
1890 $ work( itaup ),
1891 $ work( iwork ), lwork-iwork+1, ierr )
1892 irwork = ie + n
1893
1894
1895
1896
1897
1898
1899
1900 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ), a,
1901 $ lda, u, ldu, cdum, 1, rwork( irwork ),
1902 $ info )
1903
1904 END IF
1905
1906 ELSE IF( wntvas ) THEN
1907
1908
1909
1910
1911
1912
1913 IF( lwork.GE.n*n+max( n+m, 3*n ) ) THEN
1914
1915
1916
1917 iu = 1
1918 IF( lwork.GE.wrkbl+lda*n ) THEN
1919
1920
1921
1922 ldwrku = lda
1923 ELSE
1924
1925
1926
1927 ldwrku = n
1928 END IF
1929 itau = iu + ldwrku*n
1930 iwork = itau + n
1931
1932
1933
1934
1935
1936 CALL cgeqrf( m, n, a, lda, work( itau ),
1937 $ work( iwork ), lwork-iwork+1, ierr )
1938 CALL clacpy(
'L', m, n, a, lda, u, ldu )
1939
1940
1941
1942
1943
1944 CALL cungqr( m, m, n, u, ldu, work( itau ),
1945 $ work( iwork ), lwork-iwork+1, ierr )
1946
1947
1948
1949 CALL clacpy(
'U', n, n, a, lda, work( iu ),
1950 $ ldwrku )
1951 CALL claset(
'L', n-1, n-1, czero, czero,
1952 $ work( iu+1 ), ldwrku )
1953 ie = 1
1954 itauq = itau
1955 itaup = itauq + n
1956 iwork = itaup + n
1957
1958
1959
1960
1961
1962 CALL cgebrd( n, n, work( iu ), ldwrku, s,
1963 $ rwork( ie ), work( itauq ),
1964 $ work( itaup ), work( iwork ),
1965 $ lwork-iwork+1, ierr )
1966 CALL clacpy(
'U', n, n, work( iu ), ldwrku, vt,
1967 $ ldvt )
1968
1969
1970
1971
1972
1973 CALL cungbr(
'Q', n, n, n, work( iu ), ldwrku,
1974 $ work( itauq ), work( iwork ),
1975 $ lwork-iwork+1, ierr )
1976
1977
1978
1979
1980
1981
1982 CALL cungbr(
'P', n, n, n, vt, ldvt,
1983 $ work( itaup ),
1984 $ work( iwork ), lwork-iwork+1, ierr )
1985 irwork = ie + n
1986
1987
1988
1989
1990
1991
1992
1993 CALL cbdsqr(
'U', n, n, n, 0, s, rwork( ie ),
1994 $ vt,
1995 $ ldvt, work( iu ), ldwrku, cdum, 1,
1996 $ rwork( irwork ), info )
1997
1998
1999
2000
2001
2002
2003 CALL cgemm(
'N',
'N', m, n, n, cone, u, ldu,
2004 $ work( iu ), ldwrku, czero, a, lda )
2005
2006
2007
2008 CALL clacpy(
'F', m, n, a, lda, u, ldu )
2009
2010 ELSE
2011
2012
2013
2014 itau = 1
2015 iwork = itau + n
2016
2017
2018
2019
2020
2021 CALL cgeqrf( m, n, a, lda, work( itau ),
2022 $ work( iwork ), lwork-iwork+1, ierr )
2023 CALL clacpy(
'L', m, n, a, lda, u, ldu )
2024
2025
2026
2027
2028
2029 CALL cungqr( m, m, n, u, ldu, work( itau ),
2030 $ work( iwork ), lwork-iwork+1, ierr )
2031
2032
2033
2034 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
2035 IF( n.GT.1 )
2036 $
CALL claset(
'L', n-1, n-1, czero, czero,
2037 $ vt( 2, 1 ), ldvt )
2038 ie = 1
2039 itauq = itau
2040 itaup = itauq + n
2041 iwork = itaup + n
2042
2043
2044
2045
2046
2047 CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
2048 $ work( itauq ), work( itaup ),
2049 $ work( iwork ), lwork-iwork+1, ierr )
2050
2051
2052
2053
2054
2055
2056 CALL cunmbr(
'Q',
'R',
'N', m, n, n, vt, ldvt,
2057 $ work( itauq ), u, ldu, work( iwork ),
2058 $ lwork-iwork+1, ierr )
2059
2060
2061
2062
2063
2064 CALL cungbr(
'P', n, n, n, vt, ldvt,
2065 $ work( itaup ),
2066 $ work( iwork ), lwork-iwork+1, ierr )
2067 irwork = ie + n
2068
2069
2070
2071
2072
2073
2074
2075 CALL cbdsqr(
'U', n, n, m, 0, s, rwork( ie ),
2076 $ vt,
2077 $ ldvt, u, ldu, cdum, 1,
2078 $ rwork( irwork ), info )
2079
2080 END IF
2081
2082 END IF
2083
2084 END IF
2085
2086 ELSE
2087
2088
2089
2090
2091
2092
2093 ie = 1
2094 itauq = 1
2095 itaup = itauq + n
2096 iwork = itaup + n
2097
2098
2099
2100
2101
2102 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2103 $ work( itaup ), work( iwork ), lwork-iwork+1,
2104 $ ierr )
2105 IF( wntuas ) THEN
2106
2107
2108
2109
2110
2111
2112 CALL clacpy(
'L', m, n, a, lda, u, ldu )
2113 IF( wntus )
2114 $ ncu = n
2115 IF( wntua )
2116 $ ncu = m
2117 CALL cungbr(
'Q', m, ncu, n, u, ldu, work( itauq ),
2118 $ work( iwork ), lwork-iwork+1, ierr )
2119 END IF
2120 IF( wntvas ) THEN
2121
2122
2123
2124
2125
2126
2127 CALL clacpy(
'U', n, n, a, lda, vt, ldvt )
2128 CALL cungbr(
'P', n, n, n, vt, ldvt, work( itaup ),
2129 $ work( iwork ), lwork-iwork+1, ierr )
2130 END IF
2131 IF( wntuo ) THEN
2132
2133
2134
2135
2136
2137
2138 CALL cungbr(
'Q', m, n, n, a, lda, work( itauq ),
2139 $ work( iwork ), lwork-iwork+1, ierr )
2140 END IF
2141 IF( wntvo ) THEN
2142
2143
2144
2145
2146
2147
2148 CALL cungbr(
'P', n, n, n, a, lda, work( itaup ),
2149 $ work( iwork ), lwork-iwork+1, ierr )
2150 END IF
2151 irwork = ie + n
2152 IF( wntuas .OR. wntuo )
2153 $ nru = m
2154 IF( wntun )
2155 $ nru = 0
2156 IF( wntvas .OR. wntvo )
2157 $ ncvt = n
2158 IF( wntvn )
2159 $ ncvt = 0
2160 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
2161
2162
2163
2164
2165
2166
2167
2168 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2169 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2170 $ info )
2171 ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
2172
2173
2174
2175
2176
2177
2178
2179 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2180 $ lda, u, ldu, cdum, 1, rwork( irwork ),
2181 $ info )
2182 ELSE
2183
2184
2185
2186
2187
2188
2189
2190 CALL cbdsqr(
'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2191 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2192 $ info )
2193 END IF
2194
2195 END IF
2196
2197 ELSE
2198
2199
2200
2201
2202
2203 IF( n.GE.mnthr ) THEN
2204
2205 IF( wntvn ) THEN
2206
2207
2208
2209
2210 itau = 1
2211 iwork = itau + m
2212
2213
2214
2215
2216
2217 CALL cgelqf( m, n, a, lda, work( itau ),
2218 $ work( iwork ),
2219 $ lwork-iwork+1, ierr )
2220
2221
2222
2223 CALL claset(
'U', m-1, m-1, czero, czero, a( 1, 2 ),
2224 $ lda )
2225 ie = 1
2226 itauq = 1
2227 itaup = itauq + m
2228 iwork = itaup + m
2229
2230
2231
2232
2233
2234 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2235 $ work( itauq ),
2236 $ work( itaup ), work( iwork ), lwork-iwork+1,
2237 $ ierr )
2238 IF( wntuo .OR. wntuas ) THEN
2239
2240
2241
2242
2243
2244 CALL cungbr(
'Q', m, m, m, a, lda, work( itauq ),
2245 $ work( iwork ), lwork-iwork+1, ierr )
2246 END IF
2247 irwork = ie + m
2248 nru = 0
2249 IF( wntuo .OR. wntuas )
2250 $ nru = m
2251
2252
2253
2254
2255
2256
2257 CALL cbdsqr(
'U', m, 0, nru, 0, s, rwork( ie ), cdum,
2258 $ 1,
2259 $ a, lda, cdum, 1, rwork( irwork ), info )
2260
2261
2262
2263 IF( wntuas )
2264 $
CALL clacpy(
'F', m, m, a, lda, u, ldu )
2265
2266 ELSE IF( wntvo .AND. wntun ) THEN
2267
2268
2269
2270
2271
2272 IF( lwork.GE.m*m+3*m ) THEN
2273
2274
2275
2276 ir = 1
2277 IF( lwork.GE.max( wrkbl, lda*n )+lda*m ) THEN
2278
2279
2280
2281 ldwrku = lda
2282 chunk = n
2283 ldwrkr = lda
2284 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m ) THEN
2285
2286
2287
2288 ldwrku = lda
2289 chunk = n
2290 ldwrkr = m
2291 ELSE
2292
2293
2294
2295 ldwrku = m
2296 chunk = ( lwork-m*m ) / m
2297 ldwrkr = m
2298 END IF
2299 itau = ir + ldwrkr*m
2300 iwork = itau + m
2301
2302
2303
2304
2305
2306 CALL cgelqf( m, n, a, lda, work( itau ),
2307 $ work( iwork ), lwork-iwork+1, ierr )
2308
2309
2310
2311 CALL clacpy(
'L', m, m, a, lda, work( ir ),
2312 $ ldwrkr )
2313 CALL claset(
'U', m-1, m-1, czero, czero,
2314 $ work( ir+ldwrkr ), ldwrkr )
2315
2316
2317
2318
2319
2320 CALL cunglq( m, n, m, a, lda, work( itau ),
2321 $ work( iwork ), lwork-iwork+1, ierr )
2322 ie = 1
2323 itauq = itau
2324 itaup = itauq + m
2325 iwork = itaup + m
2326
2327
2328
2329
2330
2331 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
2332 $ rwork( ie ),
2333 $ work( itauq ), work( itaup ),
2334 $ work( iwork ), lwork-iwork+1, ierr )
2335
2336
2337
2338
2339
2340 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2341 $ work( itaup ), work( iwork ),
2342 $ lwork-iwork+1, ierr )
2343 irwork = ie + m
2344
2345
2346
2347
2348
2349
2350 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2351 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2352 $ rwork( irwork ), info )
2353 iu = itauq
2354
2355
2356
2357
2358
2359
2360 DO 30 i = 1, n, chunk
2361 blk = min( n-i+1, chunk )
2362 CALL cgemm(
'N',
'N', m, blk, m, cone,
2363 $ work( ir ),
2364 $ ldwrkr, a( 1, i ), lda, czero,
2365 $ work( iu ), ldwrku )
2366 CALL clacpy(
'F', m, blk, work( iu ), ldwrku,
2367 $ a( 1, i ), lda )
2368 30 CONTINUE
2369
2370 ELSE
2371
2372
2373
2374 ie = 1
2375 itauq = 1
2376 itaup = itauq + m
2377 iwork = itaup + m
2378
2379
2380
2381
2382
2383 CALL cgebrd( m, n, a, lda, s, rwork( ie ),
2384 $ work( itauq ), work( itaup ),
2385 $ work( iwork ), lwork-iwork+1, ierr )
2386
2387
2388
2389
2390
2391 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
2392 $ work( iwork ), lwork-iwork+1, ierr )
2393 irwork = ie + m
2394
2395
2396
2397
2398
2399
2400 CALL cbdsqr(
'L', m, n, 0, 0, s, rwork( ie ), a,
2401 $ lda,
2402 $ cdum, 1, cdum, 1, rwork( irwork ), info )
2403
2404 END IF
2405
2406 ELSE IF( wntvo .AND. wntuas ) THEN
2407
2408
2409
2410
2411
2412 IF( lwork.GE.m*m+3*m ) THEN
2413
2414
2415
2416 ir = 1
2417 IF( lwork.GE.max( wrkbl, lda*n )+lda*m ) THEN
2418
2419
2420
2421 ldwrku = lda
2422 chunk = n
2423 ldwrkr = lda
2424 ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m ) THEN
2425
2426
2427
2428 ldwrku = lda
2429 chunk = n
2430 ldwrkr = m
2431 ELSE
2432
2433
2434
2435 ldwrku = m
2436 chunk = ( lwork-m*m ) / m
2437 ldwrkr = m
2438 END IF
2439 itau = ir + ldwrkr*m
2440 iwork = itau + m
2441
2442
2443
2444
2445
2446 CALL cgelqf( m, n, a, lda, work( itau ),
2447 $ work( iwork ), lwork-iwork+1, ierr )
2448
2449
2450
2451 CALL clacpy(
'L', m, m, a, lda, u, ldu )
2452 CALL claset(
'U', m-1, m-1, czero, czero, u( 1,
2453 $ 2 ),
2454 $ ldu )
2455
2456
2457
2458
2459
2460 CALL cunglq( m, n, m, a, lda, work( itau ),
2461 $ work( iwork ), lwork-iwork+1, ierr )
2462 ie = 1
2463 itauq = itau
2464 itaup = itauq + m
2465 iwork = itaup + m
2466
2467
2468
2469
2470
2471 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2472 $ work( itauq ), work( itaup ),
2473 $ work( iwork ), lwork-iwork+1, ierr )
2474 CALL clacpy(
'U', m, m, u, ldu, work( ir ),
2475 $ ldwrkr )
2476
2477
2478
2479
2480
2481 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2482 $ work( itaup ), work( iwork ),
2483 $ lwork-iwork+1, ierr )
2484
2485
2486
2487
2488
2489 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2490 $ work( iwork ), lwork-iwork+1, ierr )
2491 irwork = ie + m
2492
2493
2494
2495
2496
2497
2498
2499 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2500 $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2501 $ rwork( irwork ), info )
2502 iu = itauq
2503
2504
2505
2506
2507
2508
2509 DO 40 i = 1, n, chunk
2510 blk = min( n-i+1, chunk )
2511 CALL cgemm(
'N',
'N', m, blk, m, cone,
2512 $ work( ir ),
2513 $ ldwrkr, a( 1, i ), lda, czero,
2514 $ work( iu ), ldwrku )
2515 CALL clacpy(
'F', m, blk, work( iu ), ldwrku,
2516 $ a( 1, i ), lda )
2517 40 CONTINUE
2518
2519 ELSE
2520
2521
2522
2523 itau = 1
2524 iwork = itau + m
2525
2526
2527
2528
2529
2530 CALL cgelqf( m, n, a, lda, work( itau ),
2531 $ work( iwork ), lwork-iwork+1, ierr )
2532
2533
2534
2535 CALL clacpy(
'L', m, m, a, lda, u, ldu )
2536 CALL claset(
'U', m-1, m-1, czero, czero, u( 1,
2537 $ 2 ),
2538 $ ldu )
2539
2540
2541
2542
2543
2544 CALL cunglq( m, n, m, a, lda, work( itau ),
2545 $ work( iwork ), lwork-iwork+1, ierr )
2546 ie = 1
2547 itauq = itau
2548 itaup = itauq + m
2549 iwork = itaup + m
2550
2551
2552
2553
2554
2555 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2556 $ work( itauq ), work( itaup ),
2557 $ work( iwork ), lwork-iwork+1, ierr )
2558
2559
2560
2561
2562
2563 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
2564 $ work( itaup ), a, lda, work( iwork ),
2565 $ lwork-iwork+1, ierr )
2566
2567
2568
2569
2570
2571 CALL cungbr(
'Q', m, m, m, u, ldu, work( itauq ),
2572 $ work( iwork ), lwork-iwork+1, ierr )
2573 irwork = ie + m
2574
2575
2576
2577
2578
2579
2580
2581 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ), a,
2582 $ lda,
2583 $ u, ldu, cdum, 1, rwork( irwork ), info )
2584
2585 END IF
2586
2587 ELSE IF( wntvs ) THEN
2588
2589 IF( wntun ) THEN
2590
2591
2592
2593
2594
2595 IF( lwork.GE.m*m+3*m ) THEN
2596
2597
2598
2599 ir = 1
2600 IF( lwork.GE.wrkbl+lda*m ) THEN
2601
2602
2603
2604 ldwrkr = lda
2605 ELSE
2606
2607
2608
2609 ldwrkr = m
2610 END IF
2611 itau = ir + ldwrkr*m
2612 iwork = itau + m
2613
2614
2615
2616
2617
2618 CALL cgelqf( m, n, a, lda, work( itau ),
2619 $ work( iwork ), lwork-iwork+1, ierr )
2620
2621
2622
2623 CALL clacpy(
'L', m, m, a, lda, work( ir ),
2624 $ ldwrkr )
2625 CALL claset(
'U', m-1, m-1, czero, czero,
2626 $ work( ir+ldwrkr ), ldwrkr )
2627
2628
2629
2630
2631
2632 CALL cunglq( m, n, m, a, lda, work( itau ),
2633 $ work( iwork ), lwork-iwork+1, ierr )
2634 ie = 1
2635 itauq = itau
2636 itaup = itauq + m
2637 iwork = itaup + m
2638
2639
2640
2641
2642
2643 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
2644 $ rwork( ie ), work( itauq ),
2645 $ work( itaup ), work( iwork ),
2646 $ lwork-iwork+1, ierr )
2647
2648
2649
2650
2651
2652
2653 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
2654 $ work( itaup ), work( iwork ),
2655 $ lwork-iwork+1, ierr )
2656 irwork = ie + m
2657
2658
2659
2660
2661
2662
2663 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
2664 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2665 $ rwork( irwork ), info )
2666
2667
2668
2669
2670
2671
2672 CALL cgemm(
'N',
'N', m, n, m, cone, work( ir ),
2673 $ ldwrkr, a, lda, czero, vt, ldvt )
2674
2675 ELSE
2676
2677
2678
2679 itau = 1
2680 iwork = itau + m
2681
2682
2683
2684
2685
2686 CALL cgelqf( m, n, a, lda, work( itau ),
2687 $ work( iwork ), lwork-iwork+1, ierr )
2688
2689
2690
2691 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
2692
2693
2694
2695
2696
2697 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2698 $ work( iwork ), lwork-iwork+1, ierr )
2699 ie = 1
2700 itauq = itau
2701 itaup = itauq + m
2702 iwork = itaup + m
2703
2704
2705
2706 CALL claset(
'U', m-1, m-1, czero, czero,
2707 $ a( 1, 2 ), lda )
2708
2709
2710
2711
2712
2713 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2714 $ work( itauq ), work( itaup ),
2715 $ work( iwork ), lwork-iwork+1, ierr )
2716
2717
2718
2719
2720
2721 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2722 $ work( itaup ), vt, ldvt,
2723 $ work( iwork ), lwork-iwork+1, ierr )
2724 irwork = ie + m
2725
2726
2727
2728
2729
2730
2731 CALL cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ),
2732 $ vt,
2733 $ ldvt, cdum, 1, cdum, 1,
2734 $ rwork( irwork ), info )
2735
2736 END IF
2737
2738 ELSE IF( wntuo ) THEN
2739
2740
2741
2742
2743
2744 IF( lwork.GE.2*m*m+3*m ) THEN
2745
2746
2747
2748 iu = 1
2749 IF( lwork.GE.wrkbl+2*lda*m ) THEN
2750
2751
2752
2753 ldwrku = lda
2754 ir = iu + ldwrku*m
2755 ldwrkr = lda
2756 ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
2757
2758
2759
2760 ldwrku = lda
2761 ir = iu + ldwrku*m
2762 ldwrkr = m
2763 ELSE
2764
2765
2766
2767 ldwrku = m
2768 ir = iu + ldwrku*m
2769 ldwrkr = m
2770 END IF
2771 itau = ir + ldwrkr*m
2772 iwork = itau + m
2773
2774
2775
2776
2777
2778 CALL cgelqf( m, n, a, lda, work( itau ),
2779 $ work( iwork ), lwork-iwork+1, ierr )
2780
2781
2782
2783 CALL clacpy(
'L', m, m, a, lda, work( iu ),
2784 $ ldwrku )
2785 CALL claset(
'U', m-1, m-1, czero, czero,
2786 $ work( iu+ldwrku ), ldwrku )
2787
2788
2789
2790
2791
2792 CALL cunglq( m, n, m, a, lda, work( itau ),
2793 $ work( iwork ), lwork-iwork+1, ierr )
2794 ie = 1
2795 itauq = itau
2796 itaup = itauq + m
2797 iwork = itaup + m
2798
2799
2800
2801
2802
2803
2804
2805 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2806 $ rwork( ie ), work( itauq ),
2807 $ work( itaup ), work( iwork ),
2808 $ lwork-iwork+1, ierr )
2809 CALL clacpy(
'L', m, m, work( iu ), ldwrku,
2810 $ work( ir ), ldwrkr )
2811
2812
2813
2814
2815
2816
2817 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
2818 $ work( itaup ), work( iwork ),
2819 $ lwork-iwork+1, ierr )
2820
2821
2822
2823
2824
2825 CALL cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
2826 $ work( itauq ), work( iwork ),
2827 $ lwork-iwork+1, ierr )
2828 irwork = ie + m
2829
2830
2831
2832
2833
2834
2835
2836 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
2837 $ work( iu ), ldwrku, work( ir ),
2838 $ ldwrkr, cdum, 1, rwork( irwork ),
2839 $ info )
2840
2841
2842
2843
2844
2845
2846 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
2847 $ ldwrku, a, lda, czero, vt, ldvt )
2848
2849
2850
2851
2852
2853 CALL clacpy(
'F', m, m, work( ir ), ldwrkr, a,
2854 $ lda )
2855
2856 ELSE
2857
2858
2859
2860 itau = 1
2861 iwork = itau + m
2862
2863
2864
2865
2866
2867 CALL cgelqf( m, n, a, lda, work( itau ),
2868 $ work( iwork ), lwork-iwork+1, ierr )
2869 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
2870
2871
2872
2873
2874
2875 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2876 $ work( iwork ), lwork-iwork+1, ierr )
2877 ie = 1
2878 itauq = itau
2879 itaup = itauq + m
2880 iwork = itaup + m
2881
2882
2883
2884 CALL claset(
'U', m-1, m-1, czero, czero,
2885 $ a( 1, 2 ), lda )
2886
2887
2888
2889
2890
2891 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2892 $ work( itauq ), work( itaup ),
2893 $ work( iwork ), lwork-iwork+1, ierr )
2894
2895
2896
2897
2898
2899 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
2900 $ work( itaup ), vt, ldvt,
2901 $ work( iwork ), lwork-iwork+1, ierr )
2902
2903
2904
2905
2906
2907 CALL cungbr(
'Q', m, m, m, a, lda,
2908 $ work( itauq ),
2909 $ work( iwork ), lwork-iwork+1, ierr )
2910 irwork = ie + m
2911
2912
2913
2914
2915
2916
2917
2918 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ),
2919 $ vt,
2920 $ ldvt, a, lda, cdum, 1,
2921 $ rwork( irwork ), info )
2922
2923 END IF
2924
2925 ELSE IF( wntuas ) THEN
2926
2927
2928
2929
2930
2931
2932 IF( lwork.GE.m*m+3*m ) THEN
2933
2934
2935
2936 iu = 1
2937 IF( lwork.GE.wrkbl+lda*m ) THEN
2938
2939
2940
2941 ldwrku = lda
2942 ELSE
2943
2944
2945
2946 ldwrku = m
2947 END IF
2948 itau = iu + ldwrku*m
2949 iwork = itau + m
2950
2951
2952
2953
2954
2955 CALL cgelqf( m, n, a, lda, work( itau ),
2956 $ work( iwork ), lwork-iwork+1, ierr )
2957
2958
2959
2960 CALL clacpy(
'L', m, m, a, lda, work( iu ),
2961 $ ldwrku )
2962 CALL claset(
'U', m-1, m-1, czero, czero,
2963 $ work( iu+ldwrku ), ldwrku )
2964
2965
2966
2967
2968
2969 CALL cunglq( m, n, m, a, lda, work( itau ),
2970 $ work( iwork ), lwork-iwork+1, ierr )
2971 ie = 1
2972 itauq = itau
2973 itaup = itauq + m
2974 iwork = itaup + m
2975
2976
2977
2978
2979
2980 CALL cgebrd( m, m, work( iu ), ldwrku, s,
2981 $ rwork( ie ), work( itauq ),
2982 $ work( itaup ), work( iwork ),
2983 $ lwork-iwork+1, ierr )
2984 CALL clacpy(
'L', m, m, work( iu ), ldwrku, u,
2985 $ ldu )
2986
2987
2988
2989
2990
2991
2992 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
2993 $ work( itaup ), work( iwork ),
2994 $ lwork-iwork+1, ierr )
2995
2996
2997
2998
2999
3000 CALL cungbr(
'Q', m, m, m, u, ldu,
3001 $ work( itauq ),
3002 $ work( iwork ), lwork-iwork+1, ierr )
3003 irwork = ie + m
3004
3005
3006
3007
3008
3009
3010
3011 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3012 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3013 $ rwork( irwork ), info )
3014
3015
3016
3017
3018
3019
3020 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3021 $ ldwrku, a, lda, czero, vt, ldvt )
3022
3023 ELSE
3024
3025
3026
3027 itau = 1
3028 iwork = itau + m
3029
3030
3031
3032
3033
3034 CALL cgelqf( m, n, a, lda, work( itau ),
3035 $ work( iwork ), lwork-iwork+1, ierr )
3036 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3037
3038
3039
3040
3041
3042 CALL cunglq( m, n, m, vt, ldvt, work( itau ),
3043 $ work( iwork ), lwork-iwork+1, ierr )
3044
3045
3046
3047 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3048 CALL claset(
'U', m-1, m-1, czero, czero,
3049 $ u( 1, 2 ), ldu )
3050 ie = 1
3051 itauq = itau
3052 itaup = itauq + m
3053 iwork = itaup + m
3054
3055
3056
3057
3058
3059 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3060 $ work( itauq ), work( itaup ),
3061 $ work( iwork ), lwork-iwork+1, ierr )
3062
3063
3064
3065
3066
3067
3068 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3069 $ work( itaup ), vt, ldvt,
3070 $ work( iwork ), lwork-iwork+1, ierr )
3071
3072
3073
3074
3075
3076 CALL cungbr(
'Q', m, m, m, u, ldu,
3077 $ work( itauq ),
3078 $ work( iwork ), lwork-iwork+1, ierr )
3079 irwork = ie + m
3080
3081
3082
3083
3084
3085
3086
3087 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ),
3088 $ vt,
3089 $ ldvt, u, ldu, cdum, 1,
3090 $ rwork( irwork ), info )
3091
3092 END IF
3093
3094 END IF
3095
3096 ELSE IF( wntva ) THEN
3097
3098 IF( wntun ) THEN
3099
3100
3101
3102
3103
3104 IF( lwork.GE.m*m+max( n+m, 3*m ) ) THEN
3105
3106
3107
3108 ir = 1
3109 IF( lwork.GE.wrkbl+lda*m ) THEN
3110
3111
3112
3113 ldwrkr = lda
3114 ELSE
3115
3116
3117
3118 ldwrkr = m
3119 END IF
3120 itau = ir + ldwrkr*m
3121 iwork = itau + m
3122
3123
3124
3125
3126
3127 CALL cgelqf( m, n, a, lda, work( itau ),
3128 $ work( iwork ), lwork-iwork+1, ierr )
3129 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3130
3131
3132
3133 CALL clacpy(
'L', m, m, a, lda, work( ir ),
3134 $ ldwrkr )
3135 CALL claset(
'U', m-1, m-1, czero, czero,
3136 $ work( ir+ldwrkr ), ldwrkr )
3137
3138
3139
3140
3141
3142 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3143 $ work( iwork ), lwork-iwork+1, ierr )
3144 ie = 1
3145 itauq = itau
3146 itaup = itauq + m
3147 iwork = itaup + m
3148
3149
3150
3151
3152
3153 CALL cgebrd( m, m, work( ir ), ldwrkr, s,
3154 $ rwork( ie ), work( itauq ),
3155 $ work( itaup ), work( iwork ),
3156 $ lwork-iwork+1, ierr )
3157
3158
3159
3160
3161
3162
3163 CALL cungbr(
'P', m, m, m, work( ir ), ldwrkr,
3164 $ work( itaup ), work( iwork ),
3165 $ lwork-iwork+1, ierr )
3166 irwork = ie + m
3167
3168
3169
3170
3171
3172
3173 CALL cbdsqr(
'U', m, m, 0, 0, s, rwork( ie ),
3174 $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3175 $ rwork( irwork ), info )
3176
3177
3178
3179
3180
3181
3182 CALL cgemm(
'N',
'N', m, n, m, cone, work( ir ),
3183 $ ldwrkr, vt, ldvt, czero, a, lda )
3184
3185
3186
3187 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3188
3189 ELSE
3190
3191
3192
3193 itau = 1
3194 iwork = itau + m
3195
3196
3197
3198
3199
3200 CALL cgelqf( m, n, a, lda, work( itau ),
3201 $ work( iwork ), lwork-iwork+1, ierr )
3202 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3203
3204
3205
3206
3207
3208 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3209 $ work( iwork ), lwork-iwork+1, ierr )
3210 ie = 1
3211 itauq = itau
3212 itaup = itauq + m
3213 iwork = itaup + m
3214
3215
3216
3217 CALL claset(
'U', m-1, m-1, czero, czero,
3218 $ a( 1, 2 ), lda )
3219
3220
3221
3222
3223
3224 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3225 $ work( itauq ), work( itaup ),
3226 $ work( iwork ), lwork-iwork+1, ierr )
3227
3228
3229
3230
3231
3232
3233 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3234 $ work( itaup ), vt, ldvt,
3235 $ work( iwork ), lwork-iwork+1, ierr )
3236 irwork = ie + m
3237
3238
3239
3240
3241
3242
3243 CALL cbdsqr(
'U', m, n, 0, 0, s, rwork( ie ),
3244 $ vt,
3245 $ ldvt, cdum, 1, cdum, 1,
3246 $ rwork( irwork ), info )
3247
3248 END IF
3249
3250 ELSE IF( wntuo ) THEN
3251
3252
3253
3254
3255
3256 IF( lwork.GE.2*m*m+max( n+m, 3*m ) ) THEN
3257
3258
3259
3260 iu = 1
3261 IF( lwork.GE.wrkbl+2*lda*m ) THEN
3262
3263
3264
3265 ldwrku = lda
3266 ir = iu + ldwrku*m
3267 ldwrkr = lda
3268 ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
3269
3270
3271
3272 ldwrku = lda
3273 ir = iu + ldwrku*m
3274 ldwrkr = m
3275 ELSE
3276
3277
3278
3279 ldwrku = m
3280 ir = iu + ldwrku*m
3281 ldwrkr = m
3282 END IF
3283 itau = ir + ldwrkr*m
3284 iwork = itau + m
3285
3286
3287
3288
3289
3290 CALL cgelqf( m, n, a, lda, work( itau ),
3291 $ work( iwork ), lwork-iwork+1, ierr )
3292 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3293
3294
3295
3296
3297
3298 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3299 $ work( iwork ), lwork-iwork+1, ierr )
3300
3301
3302
3303 CALL clacpy(
'L', m, m, a, lda, work( iu ),
3304 $ ldwrku )
3305 CALL claset(
'U', m-1, m-1, czero, czero,
3306 $ work( iu+ldwrku ), ldwrku )
3307 ie = 1
3308 itauq = itau
3309 itaup = itauq + m
3310 iwork = itaup + m
3311
3312
3313
3314
3315
3316
3317
3318 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3319 $ rwork( ie ), work( itauq ),
3320 $ work( itaup ), work( iwork ),
3321 $ lwork-iwork+1, ierr )
3322 CALL clacpy(
'L', m, m, work( iu ), ldwrku,
3323 $ work( ir ), ldwrkr )
3324
3325
3326
3327
3328
3329
3330 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
3331 $ work( itaup ), work( iwork ),
3332 $ lwork-iwork+1, ierr )
3333
3334
3335
3336
3337
3338 CALL cungbr(
'Q', m, m, m, work( ir ), ldwrkr,
3339 $ work( itauq ), work( iwork ),
3340 $ lwork-iwork+1, ierr )
3341 irwork = ie + m
3342
3343
3344
3345
3346
3347
3348
3349 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3350 $ work( iu ), ldwrku, work( ir ),
3351 $ ldwrkr, cdum, 1, rwork( irwork ),
3352 $ info )
3353
3354
3355
3356
3357
3358
3359 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3360 $ ldwrku, vt, ldvt, czero, a, lda )
3361
3362
3363
3364 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3365
3366
3367
3368 CALL clacpy(
'F', m, m, work( ir ), ldwrkr, a,
3369 $ lda )
3370
3371 ELSE
3372
3373
3374
3375 itau = 1
3376 iwork = itau + m
3377
3378
3379
3380
3381
3382 CALL cgelqf( m, n, a, lda, work( itau ),
3383 $ work( iwork ), lwork-iwork+1, ierr )
3384 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3385
3386
3387
3388
3389
3390 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3391 $ work( iwork ), lwork-iwork+1, ierr )
3392 ie = 1
3393 itauq = itau
3394 itaup = itauq + m
3395 iwork = itaup + m
3396
3397
3398
3399 CALL claset(
'U', m-1, m-1, czero, czero,
3400 $ a( 1, 2 ), lda )
3401
3402
3403
3404
3405
3406 CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3407 $ work( itauq ), work( itaup ),
3408 $ work( iwork ), lwork-iwork+1, ierr )
3409
3410
3411
3412
3413
3414
3415 CALL cunmbr(
'P',
'L',
'C', m, n, m, a, lda,
3416 $ work( itaup ), vt, ldvt,
3417 $ work( iwork ), lwork-iwork+1, ierr )
3418
3419
3420
3421
3422
3423 CALL cungbr(
'Q', m, m, m, a, lda,
3424 $ work( itauq ),
3425 $ work( iwork ), lwork-iwork+1, ierr )
3426 irwork = ie + m
3427
3428
3429
3430
3431
3432
3433
3434 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ),
3435 $ vt,
3436 $ ldvt, a, lda, cdum, 1,
3437 $ rwork( irwork ), info )
3438
3439 END IF
3440
3441 ELSE IF( wntuas ) THEN
3442
3443
3444
3445
3446
3447
3448 IF( lwork.GE.m*m+max( n+m, 3*m ) ) THEN
3449
3450
3451
3452 iu = 1
3453 IF( lwork.GE.wrkbl+lda*m ) THEN
3454
3455
3456
3457 ldwrku = lda
3458 ELSE
3459
3460
3461
3462 ldwrku = m
3463 END IF
3464 itau = iu + ldwrku*m
3465 iwork = itau + m
3466
3467
3468
3469
3470
3471 CALL cgelqf( m, n, a, lda, work( itau ),
3472 $ work( iwork ), lwork-iwork+1, ierr )
3473 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3474
3475
3476
3477
3478
3479 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3480 $ work( iwork ), lwork-iwork+1, ierr )
3481
3482
3483
3484 CALL clacpy(
'L', m, m, a, lda, work( iu ),
3485 $ ldwrku )
3486 CALL claset(
'U', m-1, m-1, czero, czero,
3487 $ work( iu+ldwrku ), ldwrku )
3488 ie = 1
3489 itauq = itau
3490 itaup = itauq + m
3491 iwork = itaup + m
3492
3493
3494
3495
3496
3497 CALL cgebrd( m, m, work( iu ), ldwrku, s,
3498 $ rwork( ie ), work( itauq ),
3499 $ work( itaup ), work( iwork ),
3500 $ lwork-iwork+1, ierr )
3501 CALL clacpy(
'L', m, m, work( iu ), ldwrku, u,
3502 $ ldu )
3503
3504
3505
3506
3507
3508 CALL cungbr(
'P', m, m, m, work( iu ), ldwrku,
3509 $ work( itaup ), work( iwork ),
3510 $ lwork-iwork+1, ierr )
3511
3512
3513
3514
3515
3516 CALL cungbr(
'Q', m, m, m, u, ldu,
3517 $ work( itauq ),
3518 $ work( iwork ), lwork-iwork+1, ierr )
3519 irwork = ie + m
3520
3521
3522
3523
3524
3525
3526
3527 CALL cbdsqr(
'U', m, m, m, 0, s, rwork( ie ),
3528 $ work( iu ), ldwrku, u, ldu, cdum, 1,
3529 $ rwork( irwork ), info )
3530
3531
3532
3533
3534
3535
3536 CALL cgemm(
'N',
'N', m, n, m, cone, work( iu ),
3537 $ ldwrku, vt, ldvt, czero, a, lda )
3538
3539
3540
3541 CALL clacpy(
'F', m, n, a, lda, vt, ldvt )
3542
3543 ELSE
3544
3545
3546
3547 itau = 1
3548 iwork = itau + m
3549
3550
3551
3552
3553
3554 CALL cgelqf( m, n, a, lda, work( itau ),
3555 $ work( iwork ), lwork-iwork+1, ierr )
3556 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3557
3558
3559
3560
3561
3562 CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3563 $ work( iwork ), lwork-iwork+1, ierr )
3564
3565
3566
3567 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3568 CALL claset(
'U', m-1, m-1, czero, czero,
3569 $ u( 1, 2 ), ldu )
3570 ie = 1
3571 itauq = itau
3572 itaup = itauq + m
3573 iwork = itaup + m
3574
3575
3576
3577
3578
3579 CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3580 $ work( itauq ), work( itaup ),
3581 $ work( iwork ), lwork-iwork+1, ierr )
3582
3583
3584
3585
3586
3587
3588 CALL cunmbr(
'P',
'L',
'C', m, n, m, u, ldu,
3589 $ work( itaup ), vt, ldvt,
3590 $ work( iwork ), lwork-iwork+1, ierr )
3591
3592
3593
3594
3595
3596 CALL cungbr(
'Q', m, m, m, u, ldu,
3597 $ work( itauq ),
3598 $ work( iwork ), lwork-iwork+1, ierr )
3599 irwork = ie + m
3600
3601
3602
3603
3604
3605
3606
3607 CALL cbdsqr(
'U', m, n, m, 0, s, rwork( ie ),
3608 $ vt,
3609 $ ldvt, u, ldu, cdum, 1,
3610 $ rwork( irwork ), info )
3611
3612 END IF
3613
3614 END IF
3615
3616 END IF
3617
3618 ELSE
3619
3620
3621
3622
3623
3624
3625 ie = 1
3626 itauq = 1
3627 itaup = itauq + m
3628 iwork = itaup + m
3629
3630
3631
3632
3633
3634 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3635 $ work( itaup ), work( iwork ), lwork-iwork+1,
3636 $ ierr )
3637 IF( wntuas ) THEN
3638
3639
3640
3641
3642
3643
3644 CALL clacpy(
'L', m, m, a, lda, u, ldu )
3645 CALL cungbr(
'Q', m, m, n, u, ldu, work( itauq ),
3646 $ work( iwork ), lwork-iwork+1, ierr )
3647 END IF
3648 IF( wntvas ) THEN
3649
3650
3651
3652
3653
3654
3655 CALL clacpy(
'U', m, n, a, lda, vt, ldvt )
3656 IF( wntva )
3657 $ nrvt = n
3658 IF( wntvs )
3659 $ nrvt = m
3660 CALL cungbr(
'P', nrvt, n, m, vt, ldvt, work( itaup ),
3661 $ work( iwork ), lwork-iwork+1, ierr )
3662 END IF
3663 IF( wntuo ) THEN
3664
3665
3666
3667
3668
3669
3670 CALL cungbr(
'Q', m, m, n, a, lda, work( itauq ),
3671 $ work( iwork ), lwork-iwork+1, ierr )
3672 END IF
3673 IF( wntvo ) THEN
3674
3675
3676
3677
3678
3679
3680 CALL cungbr(
'P', m, n, m, a, lda, work( itaup ),
3681 $ work( iwork ), lwork-iwork+1, ierr )
3682 END IF
3683 irwork = ie + m
3684 IF( wntuas .OR. wntuo )
3685 $ nru = m
3686 IF( wntun )
3687 $ nru = 0
3688 IF( wntvas .OR. wntvo )
3689 $ ncvt = n
3690 IF( wntvn )
3691 $ ncvt = 0
3692 IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
3693
3694
3695
3696
3697
3698
3699
3700 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3701 $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3702 $ info )
3703 ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
3704
3705
3706
3707
3708
3709
3710
3711 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3712 $ lda, u, ldu, cdum, 1, rwork( irwork ),
3713 $ info )
3714 ELSE
3715
3716
3717
3718
3719
3720
3721
3722 CALL cbdsqr(
'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3723 $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3724 $ info )
3725 END IF
3726
3727 END IF
3728
3729 END IF
3730
3731
3732
3733 IF( iscl.EQ.1 ) THEN
3734 IF( anrm.GT.bignum )
3735 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3736 $ ierr )
3737 IF( info.NE.0 .AND. anrm.GT.bignum )
3738 $
CALL slascl(
'G', 0, 0, bignum, anrm, minmn-1, 1,
3739 $ rwork( ie ), minmn, ierr )
3740 IF( anrm.LT.smlnum )
3741 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3742 $ ierr )
3743 IF( info.NE.0 .AND. anrm.LT.smlnum )
3744 $
CALL slascl(
'G', 0, 0, smlnum, anrm, minmn-1, 1,
3745 $ rwork( ie ), minmn, ierr )
3746 END IF
3747
3748
3749
3751
3752 RETURN
3753
3754
3755
subroutine xerbla(srname, info)
subroutine cbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, rwork, info)
CBDSQR
subroutine cgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
CGEBRD
subroutine cgelqf(m, n, a, lda, tau, work, lwork, info)
CGELQF
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
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 slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL 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 cungbr(vect, m, n, k, a, lda, tau, work, lwork, info)
CUNGBR
subroutine cunglq(m, n, k, a, lda, tau, work, lwork, info)
CUNGLQ
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
subroutine cunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMBR