SRC/dgglse.f

   1 *> \brief <b> DGGLSE solves overdetermined or underdetermined systems for OTHER matrices</b>
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGGLSE + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgglse.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgglse.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgglse.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
  22 *                          INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( * ), D( * ),
  29 *      $                   WORK( * ), X( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> DGGLSE solves the linear equality-constrained least squares (LSE)
  39 *> problem:
  40 *>
  41 *>         minimize || c - A*x ||_2   subject to   B*x = d
  42 *>
  43 *> where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
  44 *> M-vector, and d is a given P-vector. It is assumed that
  45 *> P <= N <= M+P, and
  46 *>
  47 *>          rank(B) = P and  rank( (A) ) = N.
  48 *>                               ( (B) )
  49 *>
  50 *> These conditions ensure that the LSE problem has a unique solution,
  51 *> which is obtained using a generalized RQ factorization of the
  52 *> matrices (B, A) given by
  53 *>
  54 *>    B = (0 R)*Q,   A = Z*T*Q.
  55 *> \endverbatim
  56 *
  57 *  Arguments:
  58 *  ==========
  59 *
  60 *> \param[in] M
  61 *> \verbatim
  62 *>          M is INTEGER
  63 *>          The number of rows of the matrix A.  M >= 0.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] N
  67 *> \verbatim
  68 *>          N is INTEGER
  69 *>          The number of columns of the matrices A and B. N >= 0.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] P
  73 *> \verbatim
  74 *>          P is INTEGER
  75 *>          The number of rows of the matrix B. 0 <= P <= N <= M+P.
  76 *> \endverbatim
  77 *>
  78 *> \param[in,out] A
  79 *> \verbatim
  80 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
  81 *>          On entry, the M-by-N matrix A.
  82 *>          On exit, the elements on and above the diagonal of the array
  83 *>          contain the min(M,N)-by-N upper trapezoidal matrix T.
  84 *> \endverbatim
  85 *>
  86 *> \param[in] LDA
  87 *> \verbatim
  88 *>          LDA is INTEGER
  89 *>          The leading dimension of the array A. LDA >= max(1,M).
  90 *> \endverbatim
  91 *>
  92 *> \param[in,out] B
  93 *> \verbatim
  94 *>          B is DOUBLE PRECISION array, dimension (LDB,N)
  95 *>          On entry, the P-by-N matrix B.
  96 *>          On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
  97 *>          contains the P-by-P upper triangular matrix R.
  98 *> \endverbatim
  99 *>
 100 *> \param[in] LDB
 101 *> \verbatim
 102 *>          LDB is INTEGER
 103 *>          The leading dimension of the array B. LDB >= max(1,P).
 104 *> \endverbatim
 105 *>
 106 *> \param[in,out] C
 107 *> \verbatim
 108 *>          C is DOUBLE PRECISION array, dimension (M)
 109 *>          On entry, C contains the right hand side vector for the
 110 *>          least squares part of the LSE problem.
 111 *>          On exit, the residual sum of squares for the solution
 112 *>          is given by the sum of squares of elements N-P+1 to M of
 113 *>          vector C.
 114 *> \endverbatim
 115 *>
 116 *> \param[in,out] D
 117 *> \verbatim
 118 *>          D is DOUBLE PRECISION array, dimension (P)
 119 *>          On entry, D contains the right hand side vector for the
 120 *>          constrained equation.
 121 *>          On exit, D is destroyed.
 122 *> \endverbatim
 123 *>
 124 *> \param[out] X
 125 *> \verbatim
 126 *>          X is DOUBLE PRECISION array, dimension (N)
 127 *>          On exit, X is the solution of the LSE problem.
 128 *> \endverbatim
 129 *>
 130 *> \param[out] WORK
 131 *> \verbatim
 132 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 133 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] LWORK
 137 *> \verbatim
 138 *>          LWORK is INTEGER
 139 *>          The dimension of the array WORK. LWORK >= max(1,M+N+P).
 140 *>          For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
 141 *>          where NB is an upper bound for the optimal blocksizes for
 142 *>          DGEQRF, SGERQF, DORMQR and SORMRQ.
 143 *>
 144 *>          If LWORK = -1, then a workspace query is assumed; the routine
 145 *>          only calculates the optimal size of the WORK array, returns
 146 *>          this value as the first entry of the WORK array, and no error
 147 *>          message related to LWORK is issued by XERBLA.
 148 *> \endverbatim
 149 *>
 150 *> \param[out] INFO
 151 *> \verbatim
 152 *>          INFO is INTEGER
 153 *>          = 0:  successful exit.
 154 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 155 *>          = 1:  the upper triangular factor R associated with B in the
 156 *>                generalized RQ factorization of the pair (B, A) is
 157 *>                singular, so that rank(B) < P; the least squares
 158 *>                solution could not be computed.
 159 *>          = 2:  the (N-P) by (N-P) part of the upper trapezoidal factor
 160 *>                T associated with A in the generalized RQ factorization
 161 *>                of the pair (B, A) is singular, so that
 162 *>                rank( (A) ) < N; the least squares solution could not
 163 *>                    ( (B) )
 164 *>                be computed.
 165 *> \endverbatim
 166 *
 167 *  Authors:
 168 *  ========
 169 *
 170 *> \author Univ. of Tennessee
 171 *> \author Univ. of California Berkeley
 172 *> \author Univ. of Colorado Denver
 173 *> \author NAG Ltd.
 174 *
 175 *> \date November 2011
 176 *
 177 *> \ingroup doubleOTHERsolve
 178 *
 179 *  =====================================================================
 180       SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
 181      $                   INFO )
 182 *
 183 *  -- LAPACK driver routine (version 3.4.0) --
 184 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 185 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 186 *     November 2011
 187 *
 188 *     .. Scalar Arguments ..
 189       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
 190 *     ..
 191 *     .. Array Arguments ..
 192       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( * ), D( * ),
 193      $                   WORK( * ), X( * )
 194 *     ..
 195 *
 196 *  =====================================================================
 197 *
 198 *     .. Parameters ..
 199       DOUBLE PRECISION   ONE
 200       PARAMETER          ( ONE = 1.0D+0 )
 201 *     ..
 202 *     .. Local Scalars ..
 203       LOGICAL            LQUERY
 204       INTEGER            LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
 205      $                   NB4, NR
 206 *     ..
 207 *     .. External Subroutines ..
 208       EXTERNAL           DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ,
 209      $                   DTRMV, DTRTRS, XERBLA
 210 *     ..
 211 *     .. External Functions ..
 212       INTEGER            ILAENV
 213       EXTERNAL           ILAENV
 214 *     ..
 215 *     .. Intrinsic Functions ..
 216       INTRINSIC          INT, MAX, MIN
 217 *     ..
 218 *     .. Executable Statements ..
 219 *
 220 *     Test the input parameters
 221 *
 222       INFO = 0
 223       MN = MIN( M, N )
 224       LQUERY = ( LWORK.EQ.-1 )
 225       IF( M.LT.0 ) THEN
 226          INFO = -1
 227       ELSE IF( N.LT.0 ) THEN
 228          INFO = -2
 229       ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
 230          INFO = -3
 231       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 232          INFO = -5
 233       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
 234          INFO = -7
 235       END IF
 236 *
 237 *     Calculate workspace
 238 *
 239       IF( INFO.EQ.0) THEN
 240          IF( N.EQ.0 ) THEN
 241             LWKMIN = 1
 242             LWKOPT = 1
 243          ELSE
 244             NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
 245             NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
 246             NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, P, -1 )
 247             NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, P, -1 )
 248             NB = MAX( NB1, NB2, NB3, NB4 )
 249             LWKMIN = M + N + P
 250             LWKOPT = P + MN + MAX( M, N )*NB
 251          END IF
 252          WORK( 1 ) = LWKOPT
 253 *
 254          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
 255             INFO = -12
 256          END IF
 257       END IF
 258 *
 259       IF( INFO.NE.0 ) THEN
 260          CALL XERBLA( 'DGGLSE', -INFO )
 261          RETURN
 262       ELSE IF( LQUERY ) THEN
 263          RETURN
 264       END IF
 265 *
 266 *     Quick return if possible
 267 *
 268       IF( N.EQ.0 )
 269      $   RETURN
 270 *
 271 *     Compute the GRQ factorization of matrices B and A:
 272 *
 273 *            B*Q**T = (  0  T12 ) P   Z**T*A*Q**T = ( R11 R12 ) N-P
 274 *                        N-P  P                     (  0  R22 ) M+P-N
 275 *                                                      N-P  P
 276 *
 277 *     where T12 and R11 are upper triangular, and Q and Z are
 278 *     orthogonal.
 279 *
 280       CALL DGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
 281      $             WORK( P+MN+1 ), LWORK-P-MN, INFO )
 282       LOPT = WORK( P+MN+1 )
 283 *
 284 *     Update c = Z**T *c = ( c1 ) N-P
 285 *                          ( c2 ) M+P-N
 286 *
 287       CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ),
 288      $             C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO )
 289       LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
 290 *
 291 *     Solve T12*x2 = d for x2
 292 *
 293       IF( P.GT.0 ) THEN
 294          CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
 295      $                B( 1, N-P+1 ), LDB, D, P, INFO )
 296 *
 297          IF( INFO.GT.0 ) THEN
 298             INFO = 1
 299             RETURN
 300          END IF
 301 *
 302 *        Put the solution in X
 303 *
 304          CALL DCOPY( P, D, 1, X( N-P+1 ), 1 )
 305 *
 306 *        Update c1
 307 *
 308          CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
 309      $               D, 1, ONE, C, 1 )
 310       END IF
 311 *
 312 *     Solve R11*x1 = c1 for x1
 313 *
 314       IF( N.GT.P ) THEN
 315          CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
 316      $                A, LDA, C, N-P, INFO )
 317 *
 318          IF( INFO.GT.0 ) THEN
 319             INFO = 2
 320             RETURN
 321          END IF
 322 *
 323 *        Put the solutions in X
 324 *
 325          CALL DCOPY( N-P, C, 1, X, 1 )
 326       END IF
 327 *
 328 *     Compute the residual vector:
 329 *
 330       IF( M.LT.N ) THEN
 331          NR = M + P - N
 332          IF( NR.GT.0 )
 333      $      CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ),
 334      $                  LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
 335       ELSE
 336          NR = P
 337       END IF
 338       IF( NR.GT.0 ) THEN
 339          CALL DTRMV( 'Upper', 'No transpose', 'Non unit', NR,
 340      $               A( N-P+1, N-P+1 ), LDA, D, 1 )
 341          CALL DAXPY( NR, -ONE, D, 1, C( N-P+1 ), 1 )
 342       END IF
 343 *
 344 *     Backward transformation x = Q**T*x
 345 *
 346       CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X,
 347      $             N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
 348       WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
 349 *
 350       RETURN
 351 *
 352 *     End of DGGLSE
 353 *
 354       END