zgels.f

计算矩阵的经典开源库．全世界都在用它．相信你也不能例外．

         CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
         IASCL = 1
      ELSE IF( ANRM.GT.BIGNUM ) THEN
*
*        Scale matrix norm down to BIGNUM
*
         OPCNT( GELS ) = OPCNT( GELS ) + DBLE( 6*M*N )
         CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
         IASCL = 2
      ELSE IF( ANRM.EQ.ZERO ) THEN
*
*        Matrix all zero. Return zero solution.
*
         CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
         GO TO 50
      END IF
*
      BROW = M
      IF( TPSD )
     $   BROW = N
      BNRM = ZLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
      IBSCL = 0
      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
*
*        Scale matrix norm up to SMLNUM
*
         OPCNT( GELS ) = OPCNT( GELS ) + DBLE( 6*BROW*NRHS )
         CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
     $                INFO )
         IBSCL = 1
      ELSE IF( BNRM.GT.BIGNUM ) THEN
*
*        Scale matrix norm down to BIGNUM
*
         OPCNT( GELS ) = OPCNT( GELS ) + DBLE( 6*BROW*NRHS )
         CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
     $                INFO )
         IBSCL = 2
      END IF
*
      IF( M.GE.N ) THEN
*
*        compute QR factorization of A
*
         NB = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
         OPCNT( GEQRF ) = OPCNT( GEQRF ) +
     $                    DOPLA( 'ZGEQRF', M, N, 0, 0, NB )
         T1 = DSECND( )
         CALL ZGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
     $                INFO )
         T2 = DSECND( )
         TIMNG( GEQRF ) = TIMNG( GEQRF ) + ( T2-T1 )
*
*        workspace at least N, optimally N*NB
*
         IF( .NOT.TPSD ) THEN
*
*           Least-Squares Problem min || A * X - B ||
*
*           B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
*
            NB = ILAENV( 1, 'ZUNMQR', 'LC', M, NRHS, N, -1 )
            OPCNT( UNMQR ) = OPCNT( UNMQR ) +
     $                       DOPLA( 'ZUNMQR', M, NRHS, N, 0, NB )
            T1 = DSECND( )
            CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, N, A,
     $                   LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
            T2 = DSECND( )
            TIMNG( UNMQR ) = TIMNG( UNMQR ) + ( T2-T1 )
*
*           workspace at least NRHS, optimally NRHS*NB
*
*           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
*
            OPCNT( TRSM ) = OPCNT( TRSM ) +
     $                      DOPBL3( 'ZTRSM ', N, NRHS, 0 )
            T1 = DSECND( )
            CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
     $                  NRHS, CONE, A, LDA, B, LDB )
            T2 = DSECND( )
            TIMNG( TRSM ) = TIMNG( TRSM ) + ( T2-T1 )
*
            SCLLEN = N
*
         ELSE
*
*           Overdetermined system of equations A' * X = B
*
*           B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS)
*
            OPCNT( TRSM ) = OPCNT( TRSM ) +
     $                      DOPBL3( 'ZTRSM ', N, NRHS, 0 )
            T1 = DSECND( )
            CALL ZTRSM( 'Left', 'Upper', 'Conjugate transpose',
     $                  'Non-unit', N, NRHS, CONE, A, LDA, B, LDB )
            T2 = DSECND( )
            TIMNG( TRSM ) = TIMNG( TRSM ) + ( T2-T1 )
*
*           B(N+1:M,1:NRHS) = ZERO
*
            DO 20 J = 1, NRHS
               DO 10 I = N + 1, M
                  B( I, J ) = CZERO
   10          CONTINUE
   20       CONTINUE
*
*           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
*
            NB = ILAENV( 1, 'ZUNMQR', 'LN', M, NRHS, N, -1 )
            OPCNT( UNMQR ) = OPCNT( UNMQR ) +
     $                       DOPLA( 'ZUNMQR', M, NRHS, N, 0, NB )
            T1 = DSECND( )
            CALL ZUNMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
     $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
            T2 = DSECND( )
            TIMNG( UNMQR ) = TIMNG( UNMQR ) + ( T2-T1 )
*
*           workspace at least NRHS, optimally NRHS*NB
*
            SCLLEN = M
*
         END IF
*
      ELSE
*
*        Compute LQ factorization of A
*
         NB = ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
         OPCNT( GELQF ) = OPCNT( GELQF ) +
     $                    DOPLA( 'ZGELQF', M, N, 0, 0, NB )
         T1 = DSECND( )
         CALL ZGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
     $                INFO )
         T2 = DSECND( )
         TIMNG( GELQF ) = TIMNG( GELQF ) + ( T2-T1 )
*
*        workspace at least M, optimally M*NB.
*
         IF( .NOT.TPSD ) THEN
*
*           underdetermined system of equations A * X = B
*
*           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
*
            OPCNT( TRSM ) = OPCNT( TRSM ) +
     $                      DOPBL3( 'ZTRSM ', M, NRHS, 0 )
            T1 = DSECND( )
            CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', M,
     $                  NRHS, CONE, A, LDA, B, LDB )
            T2 = DSECND( )
            TIMNG( TRSM ) = TIMNG( TRSM ) + ( T2-T1 )
*
*           B(M+1:N,1:NRHS) = 0
*
            DO 40 J = 1, NRHS
               DO 30 I = M + 1, N
                  B( I, J ) = CZERO
   30          CONTINUE
   40       CONTINUE
*
*           B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS)
*
            NB = ILAENV( 1, 'ZUNMLQ', 'LC', N, NRHS, M, -1 )
            OPCNT( UNMLQ ) = OPCNT( UNMLQ ) +
     $                       DOPLA( 'ZUNMLQ', N, NRHS, M, 0, NB )
            T1 = DSECND( )
            CALL ZUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A,
     $                   LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
            T2 = DSECND( )
            TIMNG( UNMLQ ) = TIMNG( UNMLQ ) + ( T2-T1 )
*
*           workspace at least NRHS, optimally NRHS*NB
*
            SCLLEN = N
*
         ELSE
*
*           overdetermined system min || A' * X - B ||
*
*           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
*
            NB = ILAENV( 1, 'ZUNMLQ', 'LN', N, NRHS, M, -1 )
            OPCNT( UNMLQ ) = OPCNT( UNMLQ ) +
     $                       DOPLA( 'ZUNMLQ', N, NRHS, M, 0, NB )
            T1 = DSECND( )
            CALL ZUNMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
     $                   WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
     $                   INFO )
            T2 = DSECND( )
            TIMNG( UNMLQ ) = TIMNG( UNMLQ ) + ( T2-T1 )
*
*           workspace at least NRHS, optimally NRHS*NB
*
*           B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS)
*
            OPCNT( TRSM ) = OPCNT( TRSM ) +
     $                      DOPBL3( 'ZTRSM ', M, NRHS, 0 )
            T1 = DSECND( )
            CALL ZTRSM( 'Left', 'Lower', 'Conjugate transpose',
     $                  'Non-unit', M, NRHS, CONE, A, LDA, B, LDB )
            T2 = DSECND( )
            TIMNG( TRSM ) = TIMNG( TRSM ) + ( T2-T1 )
*
            SCLLEN = M
*
         END IF
*
      END IF
*
*     Undo scaling
*
      IF( IASCL.EQ.1 ) THEN
         OPCNT( GELS ) = OPCNT( GELS ) + DBLE( 6*SCLLEN*NRHS )
         CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      ELSE IF( IASCL.EQ.2 ) THEN
         OPCNT( GELS ) = OPCNT( GELS ) + DBLE( 6*SCLLEN*NRHS )
         CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      END IF
      IF( IBSCL.EQ.1 ) THEN
         OPCNT( GELS ) = OPCNT( GELS ) + DBLE( 6*SCLLEN*NRHS )
         CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      ELSE IF( IBSCL.EQ.2 ) THEN
         OPCNT( GELS ) = OPCNT( GELS ) + DBLE( 6*SCLLEN*NRHS )
         CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
     $                INFO )
      END IF
*
   50 CONTINUE
      WORK( 1 ) = DBLE( WSIZE )
*
      RETURN
*
*     End of ZGELS
*
      END