LAPACK: SRC/sgelsy.f Source File

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