LAPACK: SRC/cgelsx.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00002      $                   WORK, 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, M, N, NRHS, RANK
00011       REAL               RCOND
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            JPVT( * )
00015       REAL               RWORK( * )
00016       COMPLEX            A( LDA, * ), B( LDB, * ), WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  This routine is deprecated and has been replaced by routine CGELSY.
00023 *
00024 *  CGELSX computes the minimum-norm solution to a complex linear least
00025 *  squares problem:
00026 *      minimize || A * X - B ||
00027 *  using a complete orthogonal factorization of A.  A is an M-by-N
00028 *  matrix which may be rank-deficient.
00029 *
00030 *  Several right hand side vectors b and solution vectors x can be
00031 *  handled in a single call; they are stored as the columns of the
00032 *  M-by-NRHS right hand side matrix B and the N-by-NRHS solution
00033 *  matrix X.
00034 *
00035 *  The routine first computes a QR factorization with column pivoting:
00036 *      A * P = Q * [ R11 R12 ]
00037 *                  [  0  R22 ]
00038 *  with R11 defined as the largest leading submatrix whose estimated
00039 *  condition number is less than 1/RCOND.  The order of R11, RANK,
00040 *  is the effective rank of A.
00041 *
00042 *  Then, R22 is considered to be negligible, and R12 is annihilated
00043 *  by unitary transformations from the right, arriving at the
00044 *  complete orthogonal factorization:
00045 *     A * P = Q * [ T11 0 ] * Z
00046 *                 [  0  0 ]
00047 *  The minimum-norm solution is then
00048 *     X = P * Z' [ inv(T11)*Q1'*B ]
00049 *                [        0       ]
00050 *  where Q1 consists of the first RANK columns of Q.
00051 *
00052 *  Arguments
00053 *  =========
00054 *
00055 *  M       (input) INTEGER
00056 *          The number of rows of the matrix A.  M >= 0.
00057 *
00058 *  N       (input) INTEGER
00059 *          The number of columns of the matrix A.  N >= 0.
00060 *
00061 *  NRHS    (input) INTEGER
00062 *          The number of right hand sides, i.e., the number of
00063 *          columns of matrices B and X. NRHS >= 0.
00064 *
00065 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00066 *          On entry, the M-by-N matrix A.
00067 *          On exit, A has been overwritten by details of its
00068 *          complete orthogonal factorization.
00069 *
00070 *  LDA     (input) INTEGER
00071 *          The leading dimension of the array A.  LDA >= max(1,M).
00072 *
00073 *  B       (input/output) COMPLEX array, dimension (LDB,NRHS)
00074 *          On entry, the M-by-NRHS right hand side matrix B.
00075 *          On exit, the N-by-NRHS solution matrix X.
00076 *          If m >= n and RANK = n, the residual sum-of-squares for
00077 *          the solution in the i-th column is given by the sum of
00078 *          squares of elements N+1:M in that column.
00079 *
00080 *  LDB     (input) INTEGER
00081 *          The leading dimension of the array B. LDB >= max(1,M,N).
00082 *
00083 *  JPVT    (input/output) INTEGER array, dimension (N)
00084 *          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
00085 *          initial column, otherwise it is a free column.  Before
00086 *          the QR factorization of A, all initial columns are
00087 *          permuted to the leading positions; only the remaining
00088 *          free columns are moved as a result of column pivoting
00089 *          during the factorization.
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) REAL
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) COMPLEX array, dimension
00105 *                      (min(M,N) + max( N, 2*min(M,N)+NRHS )),
00106 *
00107 *  RWORK   (workspace) REAL array, dimension (2*N)
00108 *
00109 *  INFO    (output) INTEGER
00110 *          = 0:  successful exit
00111 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00112 *
00113 *  =====================================================================
00114 *
00115 *     .. Parameters ..
00116       INTEGER            IMAX, IMIN
00117       PARAMETER          ( IMAX = 1, IMIN = 2 )
00118       REAL               ZERO, ONE, DONE, NTDONE
00119       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, DONE = ZERO,
00120      $                   NTDONE = ONE )
00121       COMPLEX            CZERO, CONE
00122       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00123      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00124 *     ..
00125 *     .. Local Scalars ..
00126       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
00127       REAL               ANRM, BIGNUM, BNRM, SMAX, SMAXPR, SMIN, SMINPR,
00128      $                   SMLNUM
00129       COMPLEX            C1, C2, S1, S2, T1, T2
00130 *     ..
00131 *     .. External Subroutines ..
00132       EXTERNAL           CGEQPF, CLAIC1, CLASCL, CLASET, CLATZM, CTRSM,
00133      $                   CTZRQF, CUNM2R, SLABAD, XERBLA
00134 *     ..
00135 *     .. External Functions ..
00136       REAL               CLANGE, SLAMCH
00137       EXTERNAL           CLANGE, SLAMCH
00138 *     ..
00139 *     .. Intrinsic Functions ..
00140       INTRINSIC          ABS, CONJG, MAX, MIN
00141 *     ..
00142 *     .. Executable Statements ..
00143 *
00144       MN = MIN( M, N )
00145       ISMIN = MN + 1
00146       ISMAX = 2*MN + 1
00147 *
00148 *     Test the input arguments.
00149 *
00150       INFO = 0
00151       IF( M.LT.0 ) THEN
00152          INFO = -1
00153       ELSE IF( N.LT.0 ) THEN
00154          INFO = -2
00155       ELSE IF( NRHS.LT.0 ) THEN
00156          INFO = -3
00157       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00158          INFO = -5
00159       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
00160          INFO = -7
00161       END IF
00162 *
00163       IF( INFO.NE.0 ) THEN
00164          CALL XERBLA( 'CGELSX', -INFO )
00165          RETURN
00166       END IF
00167 *
00168 *     Quick return if possible
00169 *
00170       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
00171          RANK = 0
00172          RETURN
00173       END IF
00174 *
00175 *     Get machine parameters
00176 *
00177       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
00178       BIGNUM = ONE / SMLNUM
00179       CALL SLABAD( SMLNUM, BIGNUM )
00180 *
00181 *     Scale A, B if max elements outside range [SMLNUM,BIGNUM]
00182 *
00183       ANRM = CLANGE( 'M', M, N, A, LDA, RWORK )
00184       IASCL = 0
00185       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00186 *
00187 *        Scale matrix norm up to SMLNUM
00188 *
00189          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00190          IASCL = 1
00191       ELSE IF( ANRM.GT.BIGNUM ) THEN
00192 *
00193 *        Scale matrix norm down to BIGNUM
00194 *
00195          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00196          IASCL = 2
00197       ELSE IF( ANRM.EQ.ZERO ) THEN
00198 *
00199 *        Matrix all zero. Return zero solution.
00200 *
00201          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00202          RANK = 0
00203          GO TO 100
00204       END IF
00205 *
00206       BNRM = CLANGE( 'M', M, NRHS, B, LDB, RWORK )
00207       IBSCL = 0
00208       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00209 *
00210 *        Scale matrix norm up to SMLNUM
00211 *
00212          CALL CLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00213          IBSCL = 1
00214       ELSE IF( BNRM.GT.BIGNUM ) THEN
00215 *
00216 *        Scale matrix norm down to BIGNUM
00217 *
00218          CALL CLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00219          IBSCL = 2
00220       END IF
00221 *
00222 *     Compute QR factorization with column pivoting of A:
00223 *        A * P = Q * R
00224 *
00225       CALL CGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), RWORK,
00226      $             INFO )
00227 *
00228 *     complex workspace MN+N. Real workspace 2*N. Details of Householder
00229 *     rotations stored in WORK(1:MN).
00230 *
00231 *     Determine RANK using incremental condition estimation
00232 *
00233       WORK( ISMIN ) = CONE
00234       WORK( ISMAX ) = CONE
00235       SMAX = ABS( A( 1, 1 ) )
00236       SMIN = SMAX
00237       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
00238          RANK = 0
00239          CALL CLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00240          GO TO 100
00241       ELSE
00242          RANK = 1
00243       END IF
00244 *
00245    10 CONTINUE
00246       IF( RANK.LT.MN ) THEN
00247          I = RANK + 1
00248          CALL CLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
00249      $                A( I, I ), SMINPR, S1, C1 )
00250          CALL CLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
00251      $                A( I, I ), SMAXPR, S2, C2 )
00252 *
00253          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
00254             DO 20 I = 1, RANK
00255                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
00256                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
00257    20       CONTINUE
00258             WORK( ISMIN+RANK ) = C1
00259             WORK( ISMAX+RANK ) = C2
00260             SMIN = SMINPR
00261             SMAX = SMAXPR
00262             RANK = RANK + 1
00263             GO TO 10
00264          END IF
00265       END IF
00266 *
00267 *     Logically partition R = [ R11 R12 ]
00268 *                             [  0  R22 ]
00269 *     where R11 = R(1:RANK,1:RANK)
00270 *
00271 *     [R11,R12] = [ T11, 0 ] * Y
00272 *
00273       IF( RANK.LT.N )
00274      $   CALL CTZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
00275 *
00276 *     Details of Householder rotations stored in WORK(MN+1:2*MN)
00277 *
00278 *     B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
00279 *
00280       CALL CUNM2R( 'Left', 'Conjugate transpose', M, NRHS, MN, A, LDA,
00281      $             WORK( 1 ), B, LDB, WORK( 2*MN+1 ), INFO )
00282 *
00283 *     workspace NRHS
00284 *
00285 *      B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
00286 *
00287       CALL CTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
00288      $            NRHS, CONE, A, LDA, B, LDB )
00289 *
00290       DO 40 I = RANK + 1, N
00291          DO 30 J = 1, NRHS
00292             B( I, J ) = CZERO
00293    30    CONTINUE
00294    40 CONTINUE
00295 *
00296 *     B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS)
00297 *
00298       IF( RANK.LT.N ) THEN
00299          DO 50 I = 1, RANK
00300             CALL CLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
00301      $                   CONJG( WORK( MN+I ) ), B( I, 1 ),
00302      $                   B( RANK+1, 1 ), LDB, WORK( 2*MN+1 ) )
00303    50    CONTINUE
00304       END IF
00305 *
00306 *     workspace NRHS
00307 *
00308 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
00309 *
00310       DO 90 J = 1, NRHS
00311          DO 60 I = 1, N
00312             WORK( 2*MN+I ) = NTDONE
00313    60    CONTINUE
00314          DO 80 I = 1, N
00315             IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
00316                IF( JPVT( I ).NE.I ) THEN
00317                   K = I
00318                   T1 = B( K, J )
00319                   T2 = B( JPVT( K ), J )
00320    70             CONTINUE
00321                   B( JPVT( K ), J ) = T1
00322                   WORK( 2*MN+K ) = DONE
00323                   T1 = T2
00324                   K = JPVT( K )
00325                   T2 = B( JPVT( K ), J )
00326                   IF( JPVT( K ).NE.I )
00327      $               GO TO 70
00328                   B( I, J ) = T1
00329                   WORK( 2*MN+K ) = DONE
00330                END IF
00331             END IF
00332    80    CONTINUE
00333    90 CONTINUE
00334 *
00335 *     Undo scaling
00336 *
00337       IF( IASCL.EQ.1 ) THEN
00338          CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00339          CALL CLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
00340      $                INFO )
00341       ELSE IF( IASCL.EQ.2 ) THEN
00342          CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00343          CALL CLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
00344      $                INFO )
00345       END IF
00346       IF( IBSCL.EQ.1 ) THEN
00347          CALL CLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00348       ELSE IF( IBSCL.EQ.2 ) THEN
00349          CALL CLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00350       END IF
00351 *
00352   100 CONTINUE
00353 *
00354       RETURN
00355 *
00356 *     End of CGELSX
00357 *
00358       END