LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgelsx.f
Go to the documentation of this file.
00001 *> \brief <b> SGELSX solves overdetermined or underdetermined systems for GE matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGELSX + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelsx.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelsx.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelsx.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00022 *                          WORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
00026 *       REAL               RCOND
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       INTEGER            JPVT( * )
00030 *       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> This routine is deprecated and has been replaced by routine SGELSY.
00040 *>
00041 *> SGELSX computes the minimum-norm solution to a real linear least
00042 *> squares problem:
00043 *>     minimize || A * X - B ||
00044 *> using a complete orthogonal factorization of A.  A is an M-by-N
00045 *> matrix which may be rank-deficient.
00046 *>
00047 *> Several right hand side vectors b and solution vectors x can be 
00048 *> handled in a single call; they are stored as the columns of the
00049 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
00050 *> matrix X.
00051 *>
00052 *> The routine first computes a QR factorization with column pivoting:
00053 *>     A * P = Q * [ R11 R12 ]
00054 *>                 [  0  R22 ]
00055 *> with R11 defined as the largest leading submatrix whose estimated
00056 *> condition number is less than 1/RCOND.  The order of R11, RANK,
00057 *> is the effective rank of A.
00058 *>
00059 *> Then, R22 is considered to be negligible, and R12 is annihilated
00060 *> by orthogonal transformations from the right, arriving at the
00061 *> complete orthogonal factorization:
00062 *>    A * P = Q * [ T11 0 ] * Z
00063 *>                [  0  0 ]
00064 *> The minimum-norm solution is then
00065 *>    X = P * Z**T [ inv(T11)*Q1**T*B ]
00066 *>                 [        0         ]
00067 *> where Q1 consists of the first RANK columns of Q.
00068 *> \endverbatim
00069 *
00070 *  Arguments:
00071 *  ==========
00072 *
00073 *> \param[in] M
00074 *> \verbatim
00075 *>          M is INTEGER
00076 *>          The number of rows of the matrix A.  M >= 0.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] N
00080 *> \verbatim
00081 *>          N is INTEGER
00082 *>          The number of columns of the matrix A.  N >= 0.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] NRHS
00086 *> \verbatim
00087 *>          NRHS is INTEGER
00088 *>          The number of right hand sides, i.e., the number of
00089 *>          columns of matrices B and X. NRHS >= 0.
00090 *> \endverbatim
00091 *>
00092 *> \param[in,out] A
00093 *> \verbatim
00094 *>          A is REAL array, dimension (LDA,N)
00095 *>          On entry, the M-by-N matrix A.
00096 *>          On exit, A has been overwritten by details of its
00097 *>          complete orthogonal factorization.
00098 *> \endverbatim
00099 *>
00100 *> \param[in] LDA
00101 *> \verbatim
00102 *>          LDA is INTEGER
00103 *>          The leading dimension of the array A.  LDA >= max(1,M).
00104 *> \endverbatim
00105 *>
00106 *> \param[in,out] B
00107 *> \verbatim
00108 *>          B is REAL array, dimension (LDB,NRHS)
00109 *>          On entry, the M-by-NRHS right hand side matrix B.
00110 *>          On exit, the N-by-NRHS solution matrix X.
00111 *>          If m >= n and RANK = n, the residual sum-of-squares for
00112 *>          the solution in the i-th column is given by the sum of
00113 *>          squares of elements N+1:M in that column.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] LDB
00117 *> \verbatim
00118 *>          LDB is INTEGER
00119 *>          The leading dimension of the array B. LDB >= max(1,M,N).
00120 *> \endverbatim
00121 *>
00122 *> \param[in,out] JPVT
00123 *> \verbatim
00124 *>          JPVT is INTEGER array, dimension (N)
00125 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is an
00126 *>          initial column, otherwise it is a free column.  Before
00127 *>          the QR factorization of A, all initial columns are
00128 *>          permuted to the leading positions; only the remaining
00129 *>          free columns are moved as a result of column pivoting
00130 *>          during the factorization.
00131 *>          On exit, if JPVT(i) = k, then the i-th column of A*P
00132 *>          was the k-th column of A.
00133 *> \endverbatim
00134 *>
00135 *> \param[in] RCOND
00136 *> \verbatim
00137 *>          RCOND is REAL
00138 *>          RCOND is used to determine the effective rank of A, which
00139 *>          is defined as the order of the largest leading triangular
00140 *>          submatrix R11 in the QR factorization with pivoting of A,
00141 *>          whose estimated condition number < 1/RCOND.
00142 *> \endverbatim
00143 *>
00144 *> \param[out] RANK
00145 *> \verbatim
00146 *>          RANK is INTEGER
00147 *>          The effective rank of A, i.e., the order of the submatrix
00148 *>          R11.  This is the same as the order of the submatrix T11
00149 *>          in the complete orthogonal factorization of A.
00150 *> \endverbatim
00151 *>
00152 *> \param[out] WORK
00153 *> \verbatim
00154 *>          WORK is REAL array, dimension
00155 *>                      (max( min(M,N)+3*N, 2*min(M,N)+NRHS )),
00156 *> \endverbatim
00157 *>
00158 *> \param[out] INFO
00159 *> \verbatim
00160 *>          INFO is INTEGER
00161 *>          = 0:  successful exit
00162 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00163 *> \endverbatim
00164 *
00165 *  Authors:
00166 *  ========
00167 *
00168 *> \author Univ. of Tennessee 
00169 *> \author Univ. of California Berkeley 
00170 *> \author Univ. of Colorado Denver 
00171 *> \author NAG Ltd. 
00172 *
00173 *> \date November 2011
00174 *
00175 *> \ingroup realGEsolve
00176 *
00177 *  =====================================================================
00178       SUBROUTINE SGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
00179      $                   WORK, INFO )
00180 *
00181 *  -- LAPACK driver routine (version 3.4.0) --
00182 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00183 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00184 *     November 2011
00185 *
00186 *     .. Scalar Arguments ..
00187       INTEGER            INFO, LDA, LDB, M, N, NRHS, RANK
00188       REAL               RCOND
00189 *     ..
00190 *     .. Array Arguments ..
00191       INTEGER            JPVT( * )
00192       REAL               A( LDA, * ), B( LDB, * ), WORK( * )
00193 *     ..
00194 *
00195 *  =====================================================================
00196 *
00197 *     .. Parameters ..
00198       INTEGER            IMAX, IMIN
00199       PARAMETER          ( IMAX = 1, IMIN = 2 )
00200       REAL               ZERO, ONE, DONE, NTDONE
00201       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, DONE = ZERO,
00202      $                   NTDONE = ONE )
00203 *     ..
00204 *     .. Local Scalars ..
00205       INTEGER            I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
00206       REAL               ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
00207      $                   SMAXPR, SMIN, SMINPR, SMLNUM, T1, T2
00208 *     ..
00209 *     .. External Functions ..
00210       REAL               SLAMCH, SLANGE
00211       EXTERNAL           SLAMCH, SLANGE
00212 *     ..
00213 *     .. External Subroutines ..
00214       EXTERNAL           SGEQPF, SLABAD, SLAIC1, SLASCL, SLASET, SLATZM,
00215      $                   SORM2R, STRSM, STZRQF, XERBLA
00216 *     ..
00217 *     .. Intrinsic Functions ..
00218       INTRINSIC          ABS, MAX, MIN
00219 *     ..
00220 *     .. Executable Statements ..
00221 *
00222       MN = MIN( M, N )
00223       ISMIN = MN + 1
00224       ISMAX = 2*MN + 1
00225 *
00226 *     Test the input arguments.
00227 *
00228       INFO = 0
00229       IF( M.LT.0 ) THEN
00230          INFO = -1
00231       ELSE IF( N.LT.0 ) THEN
00232          INFO = -2
00233       ELSE IF( NRHS.LT.0 ) THEN
00234          INFO = -3
00235       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00236          INFO = -5
00237       ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
00238          INFO = -7
00239       END IF
00240 *
00241       IF( INFO.NE.0 ) THEN
00242          CALL XERBLA( 'SGELSX', -INFO )
00243          RETURN
00244       END IF
00245 *
00246 *     Quick return if possible
00247 *
00248       IF( MIN( M, N, NRHS ).EQ.0 ) THEN
00249          RANK = 0
00250          RETURN
00251       END IF
00252 *
00253 *     Get machine parameters
00254 *
00255       SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' )
00256       BIGNUM = ONE / SMLNUM
00257       CALL SLABAD( SMLNUM, BIGNUM )
00258 *
00259 *     Scale A, B if max elements outside range [SMLNUM,BIGNUM]
00260 *
00261       ANRM = SLANGE( 'M', M, N, A, LDA, WORK )
00262       IASCL = 0
00263       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00264 *
00265 *        Scale matrix norm up to SMLNUM
00266 *
00267          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00268          IASCL = 1
00269       ELSE IF( ANRM.GT.BIGNUM ) THEN
00270 *
00271 *        Scale matrix norm down to BIGNUM
00272 *
00273          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00274          IASCL = 2
00275       ELSE IF( ANRM.EQ.ZERO ) THEN
00276 *
00277 *        Matrix all zero. Return zero solution.
00278 *
00279          CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
00280          RANK = 0
00281          GO TO 100
00282       END IF
00283 *
00284       BNRM = SLANGE( 'M', M, NRHS, B, LDB, WORK )
00285       IBSCL = 0
00286       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00287 *
00288 *        Scale matrix norm up to SMLNUM
00289 *
00290          CALL SLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00291          IBSCL = 1
00292       ELSE IF( BNRM.GT.BIGNUM ) THEN
00293 *
00294 *        Scale matrix norm down to BIGNUM
00295 *
00296          CALL SLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00297          IBSCL = 2
00298       END IF
00299 *
00300 *     Compute QR factorization with column pivoting of A:
00301 *        A * P = Q * R
00302 *
00303       CALL SGEQPF( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ), INFO )
00304 *
00305 *     workspace 3*N. Details of Householder rotations stored
00306 *     in WORK(1:MN).
00307 *
00308 *     Determine RANK using incremental condition estimation
00309 *
00310       WORK( ISMIN ) = ONE
00311       WORK( ISMAX ) = ONE
00312       SMAX = ABS( A( 1, 1 ) )
00313       SMIN = SMAX
00314       IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
00315          RANK = 0
00316          CALL SLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
00317          GO TO 100
00318       ELSE
00319          RANK = 1
00320       END IF
00321 *
00322    10 CONTINUE
00323       IF( RANK.LT.MN ) THEN
00324          I = RANK + 1
00325          CALL SLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
00326      $                A( I, I ), SMINPR, S1, C1 )
00327          CALL SLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
00328      $                A( I, I ), SMAXPR, S2, C2 )
00329 *
00330          IF( SMAXPR*RCOND.LE.SMINPR ) THEN
00331             DO 20 I = 1, RANK
00332                WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
00333                WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
00334    20       CONTINUE
00335             WORK( ISMIN+RANK ) = C1
00336             WORK( ISMAX+RANK ) = C2
00337             SMIN = SMINPR
00338             SMAX = SMAXPR
00339             RANK = RANK + 1
00340             GO TO 10
00341          END IF
00342       END IF
00343 *
00344 *     Logically partition R = [ R11 R12 ]
00345 *                             [  0  R22 ]
00346 *     where R11 = R(1:RANK,1:RANK)
00347 *
00348 *     [R11,R12] = [ T11, 0 ] * Y
00349 *
00350       IF( RANK.LT.N )
00351      $   CALL STZRQF( RANK, N, A, LDA, WORK( MN+1 ), INFO )
00352 *
00353 *     Details of Householder rotations stored in WORK(MN+1:2*MN)
00354 *
00355 *     B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
00356 *
00357       CALL SORM2R( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
00358      $             B, LDB, WORK( 2*MN+1 ), INFO )
00359 *
00360 *     workspace NRHS
00361 *
00362 *     B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
00363 *
00364       CALL STRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
00365      $            NRHS, ONE, A, LDA, B, LDB )
00366 *
00367       DO 40 I = RANK + 1, N
00368          DO 30 J = 1, NRHS
00369             B( I, J ) = ZERO
00370    30    CONTINUE
00371    40 CONTINUE
00372 *
00373 *     B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
00374 *
00375       IF( RANK.LT.N ) THEN
00376          DO 50 I = 1, RANK
00377             CALL SLATZM( 'Left', N-RANK+1, NRHS, A( I, RANK+1 ), LDA,
00378      $                   WORK( MN+I ), B( I, 1 ), B( RANK+1, 1 ), LDB,
00379      $                   WORK( 2*MN+1 ) )
00380    50    CONTINUE
00381       END IF
00382 *
00383 *     workspace NRHS
00384 *
00385 *     B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
00386 *
00387       DO 90 J = 1, NRHS
00388          DO 60 I = 1, N
00389             WORK( 2*MN+I ) = NTDONE
00390    60    CONTINUE
00391          DO 80 I = 1, N
00392             IF( WORK( 2*MN+I ).EQ.NTDONE ) THEN
00393                IF( JPVT( I ).NE.I ) THEN
00394                   K = I
00395                   T1 = B( K, J )
00396                   T2 = B( JPVT( K ), J )
00397    70             CONTINUE
00398                   B( JPVT( K ), J ) = T1
00399                   WORK( 2*MN+K ) = DONE
00400                   T1 = T2
00401                   K = JPVT( K )
00402                   T2 = B( JPVT( K ), J )
00403                   IF( JPVT( K ).NE.I )
00404      $               GO TO 70
00405                   B( I, J ) = T1
00406                   WORK( 2*MN+K ) = DONE
00407                END IF
00408             END IF
00409    80    CONTINUE
00410    90 CONTINUE
00411 *
00412 *     Undo scaling
00413 *
00414       IF( IASCL.EQ.1 ) THEN
00415          CALL SLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00416          CALL SLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
00417      $                INFO )
00418       ELSE IF( IASCL.EQ.2 ) THEN
00419          CALL SLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00420          CALL SLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
00421      $                INFO )
00422       END IF
00423       IF( IBSCL.EQ.1 ) THEN
00424          CALL SLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00425       ELSE IF( IBSCL.EQ.2 ) THEN
00426          CALL SLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00427       END IF
00428 *
00429   100 CONTINUE
00430 *
00431       RETURN
00432 *
00433 *     End of SGELSX
00434 *
00435       END
 All Files Functions