LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgglse.f
Go to the documentation of this file.
00001 *> \brief <b> SGGLSE solves overdetermined or underdetermined systems for OTHER 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 SGGLSE + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgglse.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgglse.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgglse.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
00022 *                          INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       REAL               A( LDA, * ), B( LDB, * ), C( * ), D( * ),
00029 *      $                   WORK( * ), X( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> SGGLSE solves the linear equality-constrained least squares (LSE)
00039 *> problem:
00040 *>
00041 *>         minimize || c - A*x ||_2   subject to   B*x = d
00042 *>
00043 *> where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
00044 *> M-vector, and d is a given P-vector. It is assumed that
00045 *> P <= N <= M+P, and
00046 *>
00047 *>          rank(B) = P and  rank( (A) ) = N.
00048 *>                               ( (B) )
00049 *>
00050 *> These conditions ensure that the LSE problem has a unique solution,
00051 *> which is obtained using a generalized RQ factorization of the
00052 *> matrices (B, A) given by
00053 *>
00054 *>    B = (0 R)*Q,   A = Z*T*Q.
00055 *> \endverbatim
00056 *
00057 *  Arguments:
00058 *  ==========
00059 *
00060 *> \param[in] M
00061 *> \verbatim
00062 *>          M is INTEGER
00063 *>          The number of rows of the matrix A.  M >= 0.
00064 *> \endverbatim
00065 *>
00066 *> \param[in] N
00067 *> \verbatim
00068 *>          N is INTEGER
00069 *>          The number of columns of the matrices A and B. N >= 0.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] P
00073 *> \verbatim
00074 *>          P is INTEGER
00075 *>          The number of rows of the matrix B. 0 <= P <= N <= M+P.
00076 *> \endverbatim
00077 *>
00078 *> \param[in,out] A
00079 *> \verbatim
00080 *>          A is REAL array, dimension (LDA,N)
00081 *>          On entry, the M-by-N matrix A.
00082 *>          On exit, the elements on and above the diagonal of the array
00083 *>          contain the min(M,N)-by-N upper trapezoidal matrix T.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] LDA
00087 *> \verbatim
00088 *>          LDA is INTEGER
00089 *>          The leading dimension of the array A. LDA >= max(1,M).
00090 *> \endverbatim
00091 *>
00092 *> \param[in,out] B
00093 *> \verbatim
00094 *>          B is REAL array, dimension (LDB,N)
00095 *>          On entry, the P-by-N matrix B.
00096 *>          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
00097 *>          contains the P-by-P upper triangular matrix R.
00098 *> \endverbatim
00099 *>
00100 *> \param[in] LDB
00101 *> \verbatim
00102 *>          LDB is INTEGER
00103 *>          The leading dimension of the array B. LDB >= max(1,P).
00104 *> \endverbatim
00105 *>
00106 *> \param[in,out] C
00107 *> \verbatim
00108 *>          C is REAL array, dimension (M)
00109 *>          On entry, C contains the right hand side vector for the
00110 *>          least squares part of the LSE problem.
00111 *>          On exit, the residual sum of squares for the solution
00112 *>          is given by the sum of squares of elements N-P+1 to M of
00113 *>          vector C.
00114 *> \endverbatim
00115 *>
00116 *> \param[in,out] D
00117 *> \verbatim
00118 *>          D is REAL array, dimension (P)
00119 *>          On entry, D contains the right hand side vector for the
00120 *>          constrained equation.
00121 *>          On exit, D is destroyed.
00122 *> \endverbatim
00123 *>
00124 *> \param[out] X
00125 *> \verbatim
00126 *>          X is REAL array, dimension (N)
00127 *>          On exit, X is the solution of the LSE problem.
00128 *> \endverbatim
00129 *>
00130 *> \param[out] WORK
00131 *> \verbatim
00132 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00133 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00134 *> \endverbatim
00135 *>
00136 *> \param[in] LWORK
00137 *> \verbatim
00138 *>          LWORK is INTEGER
00139 *>          The dimension of the array WORK. LWORK >= max(1,M+N+P).
00140 *>          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
00141 *>          where NB is an upper bound for the optimal blocksizes for
00142 *>          SGEQRF, SGERQF, SORMQR and SORMRQ.
00143 *>
00144 *>          If LWORK = -1, then a workspace query is assumed; the routine
00145 *>          only calculates the optimal size of the WORK array, returns
00146 *>          this value as the first entry of the WORK array, and no error
00147 *>          message related to LWORK is issued by XERBLA.
00148 *> \endverbatim
00149 *>
00150 *> \param[out] INFO
00151 *> \verbatim
00152 *>          INFO is INTEGER
00153 *>          = 0:  successful exit.
00154 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00155 *>          = 1:  the upper triangular factor R associated with B in the
00156 *>                generalized RQ factorization of the pair (B, A) is
00157 *>                singular, so that rank(B) < P; the least squares
00158 *>                solution could not be computed.
00159 *>          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
00160 *>                T associated with A in the generalized RQ factorization
00161 *>                of the pair (B, A) is singular, so that
00162 *>                rank( (A) ) < N; the least squares solution could not
00163 *>                    ( (B) )
00164 *>                be computed.
00165 *> \endverbatim
00166 *
00167 *  Authors:
00168 *  ========
00169 *
00170 *> \author Univ. of Tennessee 
00171 *> \author Univ. of California Berkeley 
00172 *> \author Univ. of Colorado Denver 
00173 *> \author NAG Ltd. 
00174 *
00175 *> \date November 2011
00176 *
00177 *> \ingroup realOTHERsolve
00178 *
00179 *  =====================================================================
00180       SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
00181      $                   INFO )
00182 *
00183 *  -- LAPACK driver routine (version 3.4.0) --
00184 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00185 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00186 *     November 2011
00187 *
00188 *     .. Scalar Arguments ..
00189       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
00190 *     ..
00191 *     .. Array Arguments ..
00192       REAL               A( LDA, * ), B( LDB, * ), C( * ), D( * ),
00193      $                   WORK( * ), X( * )
00194 *     ..
00195 *
00196 *  =====================================================================
00197 *
00198 *     .. Parameters ..
00199       REAL               ONE
00200       PARAMETER          ( ONE = 1.0E+0 )
00201 *     ..
00202 *     .. Local Scalars ..
00203       LOGICAL            LQUERY
00204       INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
00205      $                   NB4, NR
00206 *     ..
00207 *     .. External Subroutines ..
00208       EXTERNAL           SAXPY, SCOPY, SGEMV, SGGRQF, SORMQR, SORMRQ,
00209      $                   STRMV, STRTRS, XERBLA
00210 *     ..
00211 *     .. External Functions ..
00212       INTEGER            ILAENV
00213       EXTERNAL           ILAENV 
00214 *     ..
00215 *     .. Intrinsic Functions ..
00216       INTRINSIC          INT, MAX, MIN
00217 *     ..
00218 *     .. Executable Statements ..
00219 *
00220 *     Test the input parameters
00221 *
00222       INFO = 0
00223       MN = MIN( M, N )
00224       LQUERY = ( LWORK.EQ.-1 )
00225       IF( M.LT.0 ) THEN
00226          INFO = -1
00227       ELSE IF( N.LT.0 ) THEN
00228          INFO = -2
00229       ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
00230          INFO = -3
00231       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00232          INFO = -5
00233       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
00234          INFO = -7
00235       END IF
00236 *
00237 *     Calculate workspace
00238 *
00239       IF( INFO.EQ.0) THEN
00240          IF( N.EQ.0 ) THEN
00241             LWKMIN = 1
00242             LWKOPT = 1
00243          ELSE
00244             NB1 = ILAENV( 1, 'SGEQRF', ' ', M, N, -1, -1 )
00245             NB2 = ILAENV( 1, 'SGERQF', ' ', M, N, -1, -1 )
00246             NB3 = ILAENV( 1, 'SORMQR', ' ', M, N, P, -1 )
00247             NB4 = ILAENV( 1, 'SORMRQ', ' ', M, N, P, -1 )
00248             NB = MAX( NB1, NB2, NB3, NB4 )
00249             LWKMIN = M + N + P
00250             LWKOPT = P + MN + MAX( M, N )*NB
00251          END IF
00252          WORK( 1 ) = LWKOPT
00253 *
00254          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00255             INFO = -12
00256          END IF
00257       END IF
00258 *
00259       IF( INFO.NE.0 ) THEN
00260          CALL XERBLA( 'SGGLSE', -INFO )
00261          RETURN
00262       ELSE IF( LQUERY ) THEN
00263          RETURN
00264       END IF
00265 *
00266 *     Quick return if possible
00267 *
00268       IF( N.EQ.0 )
00269      $   RETURN
00270 *
00271 *     Compute the GRQ factorization of matrices B and A:
00272 *
00273 *            B*Q**T = (  0  T12 ) P   Z**T*A*Q**T = ( R11 R12 ) N-P
00274 *                        N-P  P                     (  0  R22 ) M+P-N
00275 *                                                      N-P  P
00276 *
00277 *     where T12 and R11 are upper triangular, and Q and Z are
00278 *     orthogonal.
00279 *
00280       CALL SGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
00281      $             WORK( P+MN+1 ), LWORK-P-MN, INFO )
00282       LOPT = WORK( P+MN+1 )
00283 *
00284 *     Update c = Z**T *c = ( c1 ) N-P
00285 *                          ( c2 ) M+P-N
00286 *
00287       CALL SORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ),
00288      $             C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO )
00289       LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
00290 *
00291 *     Solve T12*x2 = d for x2
00292 *
00293       IF( P.GT.0 ) THEN
00294          CALL STRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
00295      $                B( 1, N-P+1 ), LDB, D, P, INFO )
00296 *
00297          IF( INFO.GT.0 ) THEN
00298             INFO = 1
00299             RETURN
00300          END IF
00301 *
00302 *        Put the solution in X
00303 *
00304          CALL SCOPY( P, D, 1, X( N-P+1 ), 1 )
00305 *
00306 *        Update c1
00307 *
00308          CALL SGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
00309      $               D, 1, ONE, C, 1 )
00310       END IF
00311 *
00312 *     Solve R11*x1 = c1 for x1
00313 *
00314       IF( N.GT.P ) THEN
00315          CALL STRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
00316      $                A, LDA, C, N-P, INFO )
00317 *
00318          IF( INFO.GT.0 ) THEN
00319             INFO = 2
00320             RETURN
00321          END IF
00322 *
00323 *        Put the solutions in X
00324 *
00325          CALL SCOPY( N-P, C, 1, X, 1 )
00326       END IF
00327 *
00328 *     Compute the residual vector:
00329 *
00330       IF( M.LT.N ) THEN
00331          NR = M + P - N
00332          IF( NR.GT.0 )
00333      $      CALL SGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ),
00334      $                  LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
00335       ELSE
00336          NR = P
00337       END IF
00338       IF( NR.GT.0 ) THEN
00339          CALL STRMV( 'Upper', 'No transpose', 'Non unit', NR,
00340      $               A( N-P+1, N-P+1 ), LDA, D, 1 )
00341          CALL SAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 )
00342       END IF
00343 *
00344 *     Backward transformation x = Q**T*x
00345 *
00346       CALL SORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X,
00347      $             N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
00348       WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
00349 *
00350       RETURN
00351 *
00352 *     End of SGGLSE
00353 *
00354       END
 All Files Functions