LAPACK: SRC/zgelsy.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00002      $                   WORK, LWORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
00011       DOUBLE PRECISION   RCOND
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            JPVT( * )
00015       DOUBLE PRECISION   RWORK( * )
00016       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  ZGELSY computes the minimum-norm solution to a complex linear least
00023 *  squares problem:
00024 *      minimize || A * X - B ||
00025 *  using a complete orthogonal factorization of A.  A is an M-by-N
00026 *  matrix which may be rank-deficient.
00027 *
00028 *  Several right hand side vectors b and solution vectors x can be
00029 *  handled in a single call; they are stored as the columns of the
00030 *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
00031 *  matrix X.
00032 *
00033 *  The routine first computes a QR factorization with column pivoting:
00034 *      A * P = Q * [ R11 R12 ]
00035 *                  [  0  R22 ]
00036 *  with R11 defined as the largest leading submatrix whose estimated
00037 *  condition number is less than 1/RCOND.  The order of R11, RANK,
00038 *  is the effective rank of A.
00039 *
00040 *  Then, R22 is considered to be negligible, and R12 is annihilated
00041 *  by unitary transformations from the right, arriving at the
00042 *  complete orthogonal factorization:
00043 *     A * P = Q * [ T11 0 ] * Z
00044 *                 [  0  0 ]
00045 *  The minimum-norm solution is then
00046 *     X = P * Z' [ inv(T11)*Q1'*B ]
00047 *                [        0       ]
00048 *  where Q1 consists of the first RANK columns of Q.
00049 *
00050 *  This routine is basically identical to the original xGELSX except
00051 *  three differences:
00052 *    o The permutation of matrix B (the right hand side) is faster and
00053 *      more simple.
00054 *    o The call to the subroutine xGEQPF has been substituted by the
00055 *      the call to the subroutine xGEQP3. This subroutine is a Blas-3
00056 *      version of the QR factorization with column pivoting.
00057 *    o Matrix B (the right hand side) is updated with Blas-3.
00058 *
00059 *  Arguments
00060 *  =========
00061 *
00062 *  M       (input) INTEGER
00063 *          The number of rows of the matrix A.  M >= 0.
00064 *
00065 *  N       (input) INTEGER
00066 *          The number of columns of the matrix A.  N >= 0.
00067 *
00068 *  NRHS    (input) INTEGER
00069 *          The number of right hand sides, i.e., the number of
00070 *          columns of matrices B and X. NRHS >= 0.
00071 *
00072 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00073 *          On entry, the M-by-N matrix A.
00074 *          On exit, A has been overwritten by details of its
00075 *          complete orthogonal factorization.
00076 *
00077 *  LDA     (input) INTEGER
00078 *          The leading dimension of the array A.  LDA >= max(1,M).
00079 *
00080 *  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
00081 *          On entry, the M-by-NRHS right hand side matrix B.
00082 *          On exit, the N-by-NRHS solution matrix X.
00083 *
00084 *  LDB     (input) INTEGER
00085 *          The leading dimension of the array B. LDB >= max(1,M,N).
00086 *
00087 *  JPVT    (input/output) INTEGER array, dimension (N)
00088 *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00089 *          to the front of AP, otherwise column i is a free column.
00090 *          On exit, if JPVT(i) = k, then the i-th column of A*P
00091 *          was the k-th column of A.
00092 *
00093 *  RCOND   (input) DOUBLE PRECISION
00094 *          RCOND is used to determine the effective rank of A, which
00095 *          is defined as the order of the largest leading triangular
00096 *          submatrix R11 in the QR factorization with pivoting of A,
00097 *          whose estimated condition number < 1/RCOND.
00098 *
00099 *  RANK    (output) INTEGER
00100 *          The effective rank of A, i.e., the order of the submatrix
00101 *          R11.  This is the same as the order of the submatrix T11
00102 *          in the complete orthogonal factorization of A.
00103 *
00104 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
00105 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00106 *
00107 *  LWORK   (input) INTEGER
00108 *          The dimension of the array WORK.
00109 *          The unblocked strategy requires that:
00110 *            LWORK >= MN + MAX( 2*MN, N+1, MN+NRHS )
00111 *          where MN = min(M,N).
00112 *          The block algorithm requires that:
00113 *            LWORK >= MN + MAX( 2*MN, NB*(N+1), MN+MN*NB, MN+NB*NRHS )
00114 *          where NB is an upper bound on the blocksize returned
00115 *          by ILAENV for the routines ZGEQP3, ZTZRZF, CTZRQF, ZUNMQR,
00116 *          and ZUNMRZ.
00117 *
00118 *          If LWORK = -1, then a workspace query is assumed; the routine
00119 *          only calculates the optimal size of the WORK array, returns
00120 *          this value as the first entry of the WORK array, and no error
00121 *          message related to LWORK is issued by XERBLA.
00122 *
00123 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
00124 *
00125 *  INFO    (output) INTEGER
00126 *          = 0: successful exit
00127 *          < 0: if INFO = -i, the i-th argument had an illegal value
00128 *
00129 *  Further Details
00130 *  ===============
00131 *
00132 *  Based on contributions by
00133 *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
00134 *    E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00135 *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00136 *
00137 *  =====================================================================
00138 *
00139 *     .. Parameters ..
00140       INTEGER            IMAX, IMIN
00141       PARAMETER          ( IMAX = 1, IMIN = 2 )
00142       DOUBLE PRECISION   ZERO, ONE
00143       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00144       COMPLEX*16         CZERO, CONE
00145       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00146      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00147 *     ..
00148 *     .. Local Scalars ..
00149       LOGICAL            LQUERY
00150       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKOPT, MN,
00151      $                   NB, NB1, NB2, NB3, NB4
00152       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
00153      $                   SMLNUM, WSIZE
00154       COMPLEX*16         C1, C2, S1, S2
00155 *     ..
00156 *     .. External Subroutines ..
00157       EXTERNAL           DLABAD, XERBLA, ZCOPY, ZGEQP3, ZLAIC1, ZLASCL,
00158      $                   ZLASET, ZTRSM, ZTZRZF, ZUNMQR, ZUNMRZ
00159 *     ..
00160 *     .. External Functions ..
00161       INTEGER            ILAENV
00162       DOUBLE PRECISION   DLAMCH, ZLANGE
00163       EXTERNAL           ILAENV, DLAMCH, ZLANGE
00164 *     ..
00165 *     .. Intrinsic Functions ..
00166       INTRINSIC          ABS, DBLE, DCMPLX, MAX, MIN
00167 *     ..
00168 *     .. Executable Statements ..
00169 *
00170       MN = MIN( M, N )
00171       ISMIN = MN + 1
00172       ISMAX = 2*MN + 1
00173 *
00174 *     Test the input arguments.
00175 *
00176       INFO = 0
00177       NB1 = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
00178       NB2 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
00179       NB3 = ILAENV( 1, 'ZUNMQR', ' ', M, N, NRHS, -1 )
00180       NB4 = ILAENV( 1, 'ZUNMRQ', ' ', M, N, NRHS, -1 )
00181       NB = MAX( NB1, NB2, NB3, NB4 )
00182       LWKOPT = MAX( 1, MN+2*N+NB*( N+1 ), 2*MN+NB*NRHS )
00183       WORK( 1 ) = DCMPLX( LWKOPT )
00184       LQUERY = ( LWORK.EQ.-1 )
00185       IF( M.LT.0 ) THEN
00186          INFO = -1
00187       ELSE IF( N.LT.0 ) THEN
00188          INFO = -2
00189       ELSE IF( NRHS.LT.0 ) THEN
00190          INFO = -3
00191       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00192          INFO = -5
00193       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
00194          INFO = -7
00195       ELSE IF( LWORK.LT.( MN+MAX( 2*MN, N+1, MN+NRHS ) ) .AND. .NOT.
00196      $         LQUERY ) THEN
00197          INFO = -12
00198       END IF
00199 *
00200       IF( INFO.NE.0 ) THEN
00201          CALL XERBLA( 'ZGELSY', -INFO )
00202          RETURN
00203       ELSE IF( LQUERY ) THEN
00204          RETURN
00205       END IF
00206 *
00207 *     Quick return if possible
00208 *
00209       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
00210          RANK = 0
00211          RETURN
00212       END IF
00213 *
00214 *     Get machine parameters
00215 *
00216       SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
00217       BIGNUM = ONE / SMLNUM
00218       CALL DLABAD( SMLNUM, BIGNUM )
00219 *
00220 *     Scale A, B if max entries outside range [SMLNUM,BIGNUM]
00221 *
00222       ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
00223       IASCL = 0
00224       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00225 *
00226 *        Scale matrix norm up to SMLNUM
00227 *
00228          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00229          IASCL = 1
00230       ELSE IF( ANRM.GT.BIGNUM ) THEN
00231 *
00232 *        Scale matrix norm down to BIGNUM
00233 *
00234          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00235          IASCL = 2
00236       ELSE IF( ANRM.EQ.ZERO ) THEN
00237 *
00238 *        Matrix all zero. Return zero solution.
00239 *
00240          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00241          RANK = 0
00242          GO TO 70
00243       END IF
00244 *
00245       BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
00246       IBSCL = 0
00247       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00248 *
00249 *        Scale matrix norm up to SMLNUM
00250 *
00251          CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00252          IBSCL = 1
00253       ELSE IF( BNRM.GT.BIGNUM ) THEN
00254 *
00255 *        Scale matrix norm down to BIGNUM
00256 *
00257          CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00258          IBSCL = 2
00259       END IF
00260 *
00261 *     Compute QR factorization with column pivoting of A:
00262 *        A * P = Q * R
00263 *
00264       CALL ZGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
00265      $             LWORK-MN, RWORK, INFO )
00266       WSIZE = MN + DBLE( WORK( MN+1 ) )
00267 *
00268 *     complex workspace: MN+NB*(N+1). real workspace 2*N.
00269 *     Details of Householder rotations stored in WORK(1:MN).
00270 *
00271 *     Determine RANK using incremental condition estimation
00272 *
00273       WORK( ISMIN ) = CONE
00274       WORK( ISMAX ) = CONE
00275       SMAX = ABS( A( 1, 1 ) )
00276       SMIN = SMAX
00277       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
00278          RANK = 0
00279          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00280          GO TO 70
00281       ELSE
00282          RANK = 1
00283       END IF
00284 *
00285    10 CONTINUE
00286       IF( RANK.LT.MN ) THEN
00287          I = RANK + 1
00288          CALL ZLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
00289      $                A( I, I ), SMINPR, S1, C1 )
00290          CALL ZLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
00291      $                A( I, I ), SMAXPR, S2, C2 )
00292 *
00293          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
00294             DO 20 I = 1, RANK
00295                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
00296                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
00297    20       CONTINUE
00298             WORK( ISMIN+RANK ) = C1
00299             WORK( ISMAX+RANK ) = C2
00300             SMIN = SMINPR
00301             SMAX = SMAXPR
00302             RANK = RANK + 1
00303             GO TO 10
00304          END IF
00305       END IF
00306 *
00307 *     complex workspace: 3*MN.
00308 *
00309 *     Logically partition R = [ R11 R12 ]
00310 *                             [  0  R22 ]
00311 *     where R11 = R(1:RANK,1:RANK)
00312 *
00313 *     [R11,R12] = [ T11, 0 ] * Y
00314 *
00315       IF( RANK.LT.N )
00316      $   CALL ZTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
00317      $                LWORK-2*MN, INFO )
00318 *
00319 *     complex workspace: 2*MN.
00320 *     Details of Householder rotations stored in WORK(MN+1:2*MN)
00321 *
00322 *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
00323 *
00324       CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
00325      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
00326       WSIZE = MAX( WSIZE, 2*MN+DBLE( WORK( 2*MN+1 ) ) )
00327 *
00328 *     complex workspace: 2*MN+NB*NRHS.
00329 *
00330 *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
00331 *
00332       CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
00333      $            NRHS, CONE, A, LDA, B, LDB )
00334 *
00335       DO 40 J = 1, NRHS
00336          DO 30 I = RANK + 1, N
00337             B( I, J ) = CZERO
00338    30    CONTINUE
00339    40 CONTINUE
00340 *
00341 *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
00342 *
00343       IF( RANK.LT.N ) THEN
00344          CALL ZUNMRZ( 'Left', 'Conjugate transpose', N, NRHS, RANK,
00345      $                N-RANK, A, LDA, WORK( MN+1 ), B, LDB,
00346      $                WORK( 2*MN+1 ), LWORK-2*MN, INFO )
00347       END IF
00348 *
00349 *     complex workspace: 2*MN+NRHS.
00350 *
00351 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
00352 *
00353       DO 60 J = 1, NRHS
00354          DO 50 I = 1, N
00355             WORK( JPVT( I ) ) = B( I, J )
00356    50    CONTINUE
00357          CALL ZCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
00358    60 CONTINUE
00359 *
00360 *     complex workspace: N.
00361 *
00362 *     Undo scaling
00363 *
00364       IF( IASCL.EQ.1 ) THEN
00365          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00366          CALL ZLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
00367      $                INFO )
00368       ELSE IF( IASCL.EQ.2 ) THEN
00369          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00370          CALL ZLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
00371      $                INFO )
00372       END IF
00373       IF( IBSCL.EQ.1 ) THEN
00374          CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00375       ELSE IF( IBSCL.EQ.2 ) THEN
00376          CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00377       END IF
00378 *
00379    70 CONTINUE
00380       WORK( 1 ) = DCMPLX( LWKOPT )
00381 *
00382       RETURN
00383 *
00384 *     End of ZGELSY
00385 *
00386       END