SRC/cgglse.f

   1 *> \brief <b> CGGLSE 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 CGGLSE + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgglse.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgglse.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgglse.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CGGLSE( 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 *       COMPLEX            A( LDA, * ), B( LDB, * ), C( * ), D( * ),
  29 *      $                   WORK( * ), X( * )
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *  =============
  35 *>
  36 *> \verbatim
  37 *>
  38 *> CGGLSE 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 COMPLEX 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 *>          CGEQRF, CGERQF, CUNMQR and CUNMRQ.
 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 complexOTHERsolve
 178 *
 179 *  =====================================================================
 180       SUBROUTINE CGGLSE( 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       COMPLEX            A( LDA, * ), B( LDB, * ), C( * ), D( * ),
 193      $                   WORK( * ), X( * )
 194 *     ..
 195 *
 196 *  =====================================================================
 197 *
 198 *     .. Parameters ..
 199       COMPLEX            CONE
 200       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+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           CAXPY, CCOPY, CGEMV, CGGRQF, CTRMV, CTRTRS,
 209      $                   CUNMQR, CUNMRQ, 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, 'CGEQRF', ' ', M, N, -1, -1 )
 245             NB2 = ILAENV( 1, 'CGERQF', ' ', M, N, -1, -1 )
 246             NB3 = ILAENV( 1, 'CUNMQR', ' ', M, N, P, -1 )
 247             NB4 = ILAENV( 1, 'CUNMRQ', ' ', 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( 'CGGLSE', -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**H = (  0  T12 ) P   Z**H*A*Q**H = ( 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 *     unitary.
 279 *
 280       CALL CGGRQF( 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**H *c = ( c1 ) N-P
 285 *                       ( c2 ) M+P-N
 286 *
 287       CALL CUNMQR( 'Left', 'Conjugate Transpose', M, 1, MN, A, LDA,
 288      $             WORK( P+1 ), C, MAX( 1, M ), WORK( P+MN+1 ),
 289      $             LWORK-P-MN, INFO )
 290       LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
 291 *
 292 *     Solve T12*x2 = d for x2
 293 *
 294       IF( P.GT.0 ) THEN
 295          CALL CTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
 296      $                B( 1, N-P+1 ), LDB, D, P, INFO )
 297 *
 298          IF( INFO.GT.0 ) THEN
 299             INFO = 1
 300             RETURN
 301          END IF
 302 *
 303 *        Put the solution in X
 304 *
 305       CALL CCOPY( P, D, 1, X( N-P+1 ), 1 )
 306 *
 307 *        Update c1
 308 *
 309          CALL CGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), LDA,
 310      $               D, 1, CONE, C, 1 )
 311       END IF
 312 *
 313 *     Solve R11*x1 = c1 for x1
 314 *
 315       IF( N.GT.P ) THEN
 316          CALL CTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
 317      $                A, LDA, C, N-P, INFO )
 318 *
 319          IF( INFO.GT.0 ) THEN
 320             INFO = 2
 321             RETURN
 322          END IF
 323 *
 324 *        Put the solutions in X
 325 *
 326          CALL CCOPY( N-P, C, 1, X, 1 )
 327       END IF
 328 *
 329 *     Compute the residual vector:
 330 *
 331       IF( M.LT.N ) THEN
 332          NR = M + P - N
 333          IF( NR.GT.0 )
 334      $      CALL CGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, M+1 ),
 335      $                  LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 )
 336       ELSE
 337          NR = P
 338       END IF
 339       IF( NR.GT.0 ) THEN
 340          CALL CTRMV( 'Upper', 'No transpose', 'Non unit', NR,
 341      $               A( N-P+1, N-P+1 ), LDA, D, 1 )
 342          CALL CAXPY( NR, -CONE, D, 1, C( N-P+1 ), 1 )
 343       END IF
 344 *
 345 *     Backward transformation x = Q**H*x
 346 *
 347       CALL CUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB,
 348      $             WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
 349       WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
 350 *
 351       RETURN
 352 *
 353 *     End of CGGLSE
 354 *
 355       END