Return to zgglse.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief <b> ZGGLSE 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 ZGGLSE + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgglse.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgglse.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgglse.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGGLSE( 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*16         A( LDA, * ), B( LDB, * ), C( * ), D( * ),
   29: *      $                   WORK( * ), X( * )
   30: *       ..
   31: *  
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> ZGGLSE 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*16 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*16 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*16 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*16 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*16 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*16 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: *>          ZGEQRF, CGERQF, ZUNMQR 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 complex16OTHERsolve
  178: *
  179: *  =====================================================================
  180:       SUBROUTINE ZGGLSE( 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*16         A( LDA, * ), B( LDB, * ), C( * ), D( * ),
  193:      $                   WORK( * ), X( * )
  194: *     ..
  195: *
  196: *  =====================================================================
  197: *
  198: *     .. Parameters ..
  199:       COMPLEX*16         CONE
  200:       PARAMETER          ( CONE = ( 1.0D+0, 0.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           XERBLA, ZAXPY, ZCOPY, ZGEMV, ZGGRQF, ZTRMV,
  209:      $                   ZTRTRS, ZUNMQR, ZUNMRQ
  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, 'ZGEQRF', ' ', M, N, -1, -1 )
  245:             NB2 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
  246:             NB3 = ILAENV( 1, 'ZUNMQR', ' ', M, N, P, -1 )
  247:             NB4 = ILAENV( 1, 'ZUNMRQ', ' ', 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( 'ZGGLSE', -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 ZGGRQF( 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 ZUNMQR( '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 ZTRTRS( '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 ZCOPY( P, D, 1, X( N-P+1 ), 1 )
  306: *
  307: *        Update c1
  308: *
  309:          CALL ZGEMV( '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 ZTRTRS( '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 ZCOPY( 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 ZGEMV( '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 ZTRMV( 'Upper', 'No transpose', 'Non unit', NR,
  341:      $               A( N-P+1, N-P+1 ), LDA, D, 1 )
  342:          CALL ZAXPY( NR, -CONE, D, 1, C( N-P+1 ), 1 )
  343:       END IF
  344: *
  345: *     Backward transformation x = Q**H*x
  346: *
  347:       CALL ZUNMRQ( '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 ZGGLSE
  354: *
  355:       END