LAPACK 3.3.0

dgglse.f

Go to the documentation of this file.
00001       SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
00002      $                   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, P
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( * ), D( * ),
00014      $                   WORK( * ), X( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  DGGLSE solves the linear equality-constrained least squares (LSE)
00021 *  problem:
00022 *
00023 *          minimize || c - A*x ||_2   subject to   B*x = d
00024 *
00025 *  where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
00026 *  M-vector, and d is a given P-vector. It is assumed that
00027 *  P <= N <= M+P, and
00028 *
00029 *           rank(B) = P and  rank( (A) ) = N.
00030 *                                ( (B) )
00031 *
00032 *  These conditions ensure that the LSE problem has a unique solution,
00033 *  which is obtained using a generalized RQ factorization of the
00034 *  matrices (B, A) given by
00035 *
00036 *     B = (0 R)*Q,   A = Z*T*Q.
00037 *
00038 *  Arguments
00039 *  =========
00040 *
00041 *  M       (input) INTEGER
00042 *          The number of rows of the matrix A.  M >= 0.
00043 *
00044 *  N       (input) INTEGER
00045 *          The number of columns of the matrices A and B. N >= 0.
00046 *
00047 *  P       (input) INTEGER
00048 *          The number of rows of the matrix B. 0 <= P <= N <= M+P.
00049 *
00050 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00051 *          On entry, the M-by-N matrix A.
00052 *          On exit, the elements on and above the diagonal of the array
00053 *          contain the min(M,N)-by-N upper trapezoidal matrix T.
00054 *
00055 *  LDA     (input) INTEGER
00056 *          The leading dimension of the array A. LDA >= max(1,M).
00057 *
00058 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
00059 *          On entry, the P-by-N matrix B.
00060 *          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
00061 *          contains the P-by-P upper triangular matrix R.
00062 *
00063 *  LDB     (input) INTEGER
00064 *          The leading dimension of the array B. LDB >= max(1,P).
00065 *
00066 *  C       (input/output) DOUBLE PRECISION array, dimension (M)
00067 *          On entry, C contains the right hand side vector for the
00068 *          least squares part of the LSE problem.
00069 *          On exit, the residual sum of squares for the solution
00070 *          is given by the sum of squares of elements N-P+1 to M of
00071 *          vector C.
00072 *
00073 *  D       (input/output) DOUBLE PRECISION array, dimension (P)
00074 *          On entry, D contains the right hand side vector for the
00075 *          constrained equation.
00076 *          On exit, D is destroyed.
00077 *
00078 *  X       (output) DOUBLE PRECISION array, dimension (N)
00079 *          On exit, X is the solution of the LSE problem.
00080 *
00081 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00082 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00083 *
00084 *  LWORK   (input) INTEGER
00085 *          The dimension of the array WORK. LWORK >= max(1,M+N+P).
00086 *          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
00087 *          where NB is an upper bound for the optimal blocksizes for
00088 *          DGEQRF, SGERQF, DORMQR and SORMRQ.
00089 *
00090 *          If LWORK = -1, then a workspace query is assumed; the routine
00091 *          only calculates the optimal size of the WORK array, returns
00092 *          this value as the first entry of the WORK array, and no error
00093 *          message related to LWORK is issued by XERBLA.
00094 *
00095 *  INFO    (output) INTEGER
00096 *          = 0:  successful exit.
00097 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00098 *          = 1:  the upper triangular factor R associated with B in the
00099 *                generalized RQ factorization of the pair (B, A) is
00100 *                singular, so that rank(B) < P; the least squares
00101 *                solution could not be computed.
00102 *          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
00103 *                T associated with A in the generalized RQ factorization
00104 *                of the pair (B, A) is singular, so that
00105 *                rank( (A) ) < N; the least squares solution could not
00106 *                    ( (B) )
00107 *                be computed.
00108 *
00109 *  =====================================================================
00110 *
00111 *     .. Parameters ..
00112       DOUBLE PRECISION   ONE
00113       PARAMETER          ( ONE = 1.0D+0 )
00114 *     ..
00115 *     .. Local Scalars ..
00116       LOGICAL            LQUERY
00117       INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
00118      $                   NB4, NR
00119 *     ..
00120 *     .. External Subroutines ..
00121       EXTERNAL           DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ,
00122      $                   DTRMV, DTRTRS, XERBLA
00123 *     ..
00124 *     .. External Functions ..
00125       INTEGER            ILAENV
00126       EXTERNAL           ILAENV
00127 *     ..
00128 *     .. Intrinsic Functions ..
00129       INTRINSIC          INT, MAX, MIN
00130 *     ..
00131 *     .. Executable Statements ..
00132 *
00133 *     Test the input parameters
00134 *
00135       INFO = 0
00136       MN = MIN( M, N )
00137       LQUERY = ( LWORK.EQ.-1 )
00138       IF( M.LT.0 ) THEN
00139          INFO = -1
00140       ELSE IF( N.LT.0 ) THEN
00141          INFO = -2
00142       ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
00143          INFO = -3
00144       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00145          INFO = -5
00146       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
00147          INFO = -7
00148       END IF
00149 *
00150 *     Calculate workspace
00151 *
00152       IF( INFO.EQ.0) THEN
00153          IF( N.EQ.0 ) THEN
00154             LWKMIN = 1
00155             LWKOPT = 1
00156          ELSE
00157             NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
00158             NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
00159             NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 )
00160             NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, P, -1 )
00161             NB = MAX( NB1, NB2, NB3, NB4 )
00162             LWKMIN = M + N + P
00163             LWKOPT = P + MN + MAX( M, N )*NB
00164          END IF
00165          WORK( 1 ) = LWKOPT
00166 *
00167          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00168             INFO = -12
00169          END IF
00170       END IF
00171 *
00172       IF( INFO.NE.0 ) THEN
00173          CALL XERBLA( 'DGGLSE', -INFO )
00174          RETURN
00175       ELSE IF( LQUERY ) THEN
00176          RETURN
00177       END IF
00178 *
00179 *     Quick return if possible
00180 *
00181       IF( N.EQ.0 )
00182      $   RETURN
00183 *
00184 *     Compute the GRQ factorization of matrices B and A:
00185 *
00186 *            B*Q' = (  0  T12 ) P   Z'*A*Q' = ( R11 R12 ) N-P
00187 *                     N-P  P                  (  0  R22 ) M+P-N
00188 *                                               N-P  P
00189 *
00190 *     where T12 and R11 are upper triangular, and Q and Z are
00191 *     orthogonal.
00192 *
00193       CALL DGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
00194      $             WORK( P+MN+1 ), LWORK-P-MN, INFO )
00195       LOPT = WORK( P+MN+1 )
00196 *
00197 *     Update c = Z'*c = ( c1 ) N-P
00198 *                       ( c2 ) M+P-N
00199 *
00200       CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ),
00201      $             C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO )
00202       LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
00203 *
00204 *     Solve T12*x2 = d for x2
00205 *
00206       IF( P.GT.0 ) THEN
00207          CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
00208      $                B( 1, N-P+1 ), LDB, D, P, INFO )
00209 *
00210          IF( INFO.GT.0 ) THEN
00211             INFO = 1
00212             RETURN
00213          END IF
00214 *
00215 *        Put the solution in X
00216 *
00217          CALL DCOPY( P, D, 1, X( N-P+1 ), 1 )
00218 *
00219 *        Update c1
00220 *
00221          CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
00222      $               D, 1, ONE, C, 1 )
00223       END IF
00224 *
00225 *     Solve R11*x1 = c1 for x1
00226 *
00227       IF( N.GT.P ) THEN
00228          CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
00229      $                A, LDA, C, N-P, INFO )
00230 *
00231          IF( INFO.GT.0 ) THEN
00232             INFO = 2
00233             RETURN
00234          END IF
00235 *
00236 *        Put the solutions in X
00237 *
00238          CALL DCOPY( N-P, C, 1, X, 1 )
00239       END IF
00240 *
00241 *     Compute the residual vector:
00242 *
00243       IF( M.LT.N ) THEN
00244          NR = M + P - N
00245          IF( NR.GT.0 )
00246      $      CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ),
00247      $                  LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
00248       ELSE
00249          NR = P
00250       END IF
00251       IF( NR.GT.0 ) THEN
00252          CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR,
00253      $               A( N-P+1, N-P+1 ), LDA, D, 1 )
00254          CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 )
00255       END IF
00256 *
00257 *     Backward transformation x = Q'*x
00258 *
00259       CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X,
00260      $             N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
00261       WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
00262 *
00263       RETURN
00264 *
00265 *     End of DGGLSE
00266 *
00267       END
 All Files Functions