Return to zggglm.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
    2:      $                   INFO )
    3: *
    4: *  -- LAPACK driver routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            INFO, LDA, LDB, LWORK, M, N, P
   11: *     ..
   12: *     .. Array Arguments ..
   13:       COMPLEX*16         A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
   14:      $                   X( * ), Y( * )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  ZGGGLM solves a general Gauss-Markov linear model (GLM) problem:
   21: *
   22: *          minimize || y ||_2   subject to   d = A*x + B*y
   23: *              x
   24: *
   25: *  where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
   26: *  given N-vector. It is assumed that M <= N <= M+P, and
   27: *
   28: *             rank(A) = M    and    rank( A B ) = N.
   29: *
   30: *  Under these assumptions, the constrained equation is always
   31: *  consistent, and there is a unique solution x and a minimal 2-norm
   32: *  solution y, which is obtained using a generalized QR factorization
   33: *  of the matrices (A, B) given by
   34: *
   35: *     A = Q*(R),   B = Q*T*Z.
   36: *           (0)
   37: *
   38: *  In particular, if matrix B is square nonsingular, then the problem
   39: *  GLM is equivalent to the following weighted linear least squares
   40: *  problem
   41: *
   42: *               minimize || inv(B)*(d-A*x) ||_2
   43: *                   x
   44: *
   45: *  where inv(B) denotes the inverse of B.
   46: *
   47: *  Arguments
   48: *  =========
   49: *
   50: *  N       (input) INTEGER
   51: *          The number of rows of the matrices A and B.  N >= 0.
   52: *
   53: *  M       (input) INTEGER
   54: *          The number of columns of the matrix A.  0 <= M <= N.
   55: *
   56: *  P       (input) INTEGER
   57: *          The number of columns of the matrix B.  P >= N-M.
   58: *
   59: *  A       (input/output) COMPLEX*16 array, dimension (LDA,M)
   60: *          On entry, the N-by-M matrix A.
   61: *          On exit, the upper triangular part of the array A contains
   62: *          the M-by-M upper triangular matrix R.
   63: *
   64: *  LDA     (input) INTEGER
   65: *          The leading dimension of the array A. LDA >= max(1,N).
   66: *
   67: *  B       (input/output) COMPLEX*16 array, dimension (LDB,P)
   68: *          On entry, the N-by-P matrix B.
   69: *          On exit, if N <= P, the upper triangle of the subarray
   70: *          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
   71: *          if N > P, the elements on and above the (N-P)th subdiagonal
   72: *          contain the N-by-P upper trapezoidal matrix T.
   73: *
   74: *  LDB     (input) INTEGER
   75: *          The leading dimension of the array B. LDB >= max(1,N).
   76: *
   77: *  D       (input/output) COMPLEX*16 array, dimension (N)
   78: *          On entry, D is the left hand side of the GLM equation.
   79: *          On exit, D is destroyed.
   80: *
   81: *  X       (output) COMPLEX*16 array, dimension (M)
   82: *  Y       (output) COMPLEX*16 array, dimension (P)
   83: *          On exit, X and Y are the solutions of the GLM problem.
   84: *
   85: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
   86: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   87: *
   88: *  LWORK   (input) INTEGER
   89: *          The dimension of the array WORK. LWORK >= max(1,N+M+P).
   90: *          For optimum performance, LWORK >= M+min(N,P)+max(N,P)*NB,
   91: *          where NB is an upper bound for the optimal blocksizes for
   92: *          ZGEQRF, ZGERQF, ZUNMQR and ZUNMRQ.
   93: *
   94: *          If LWORK = -1, then a workspace query is assumed; the routine
   95: *          only calculates the optimal size of the WORK array, returns
   96: *          this value as the first entry of the WORK array, and no error
   97: *          message related to LWORK is issued by XERBLA.
   98: *
   99: *  INFO    (output) INTEGER
  100: *          = 0:  successful exit.
  101: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  102: *          = 1:  the upper triangular factor R associated with A in the
  103: *                generalized QR factorization of the pair (A, B) is
  104: *                singular, so that rank(A) < M; the least squares
  105: *                solution could not be computed.
  106: *          = 2:  the bottom (N-M) by (N-M) part of the upper trapezoidal
  107: *                factor T associated with B in the generalized QR
  108: *                factorization of the pair (A, B) is singular, so that
  109: *                rank( A B ) < N; the least squares solution could not
  110: *                be computed.
  111: *
  112: *  ===================================================================
  113: *
  114: *     .. Parameters ..
  115:       COMPLEX*16         CZERO, CONE
  116:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
  117:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
  118: *     ..
  119: *     .. Local Scalars ..
  120:       LOGICAL            LQUERY
  121:       INTEGER            I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
  122:      $                   NB4, NP
  123: *     ..
  124: *     .. External Subroutines ..
  125:       EXTERNAL           XERBLA, ZCOPY, ZGEMV, ZGGQRF, ZTRTRS, ZUNMQR,
  126:      $                   ZUNMRQ
  127: *     ..
  128: *     .. External Functions ..
  129:       INTEGER            ILAENV
  130:       EXTERNAL           ILAENV 
  131: *     ..
  132: *     .. Intrinsic Functions ..
  133:       INTRINSIC          INT, MAX, MIN
  134: *     ..
  135: *     .. Executable Statements ..
  136: *
  137: *     Test the input parameters
  138: *
  139:       INFO = 0
  140:       NP = MIN( N, P )
  141:       LQUERY = ( LWORK.EQ.-1 )
  142:       IF( N.LT.0 ) THEN
  143:          INFO = -1
  144:       ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
  145:          INFO = -2
  146:       ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN
  147:          INFO = -3
  148:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  149:          INFO = -5
  150:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  151:          INFO = -7
  152:       END IF
  153: *
  154: *     Calculate workspace
  155: *
  156:       IF( INFO.EQ.0) THEN
  157:          IF( N.EQ.0 ) THEN
  158:             LWKMIN = 1
  159:             LWKOPT = 1
  160:          ELSE
  161:             NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 )
  162:             NB2 = ILAENV( 1, 'ZGERQF', ' ', N, M, -1, -1 )
  163:             NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 )
  164:             NB4 = ILAENV( 1, 'ZUNMRQ', ' ', N, M, P, -1 )
  165:             NB = MAX( NB1, NB2, NB3, NB4 )
  166:             LWKMIN = M + N + P
  167:             LWKOPT = M + NP + MAX( N, P )*NB
  168:          END IF
  169:          WORK( 1 ) = LWKOPT
  170: *
  171:          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
  172:             INFO = -12
  173:          END IF
  174:       END IF
  175: *
  176:       IF( INFO.NE.0 ) THEN
  177:          CALL XERBLA( 'ZGGGLM', -INFO )
  178:          RETURN
  179:       ELSE IF( LQUERY ) THEN
  180:          RETURN
  181:       END IF
  182: *
  183: *     Quick return if possible
  184: *
  185:       IF( N.EQ.0 )
  186:      $   RETURN
  187: *
  188: *     Compute the GQR factorization of matrices A and B:
  189: *
  190: *            Q'*A = ( R11 ) M,    Q'*B*Z' = ( T11   T12 ) M
  191: *                   (  0  ) N-M             (  0    T22 ) N-M
  192: *                      M                     M+P-N  N-M
  193: *
  194: *     where R11 and T22 are upper triangular, and Q and Z are
  195: *     unitary.
  196: *
  197:       CALL ZGGQRF( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ),
  198:      $             WORK( M+NP+1 ), LWORK-M-NP, INFO )
  199:       LOPT = WORK( M+NP+1 )
  200: *
  201: *     Update left-hand-side vector d = Q'*d = ( d1 ) M
  202: *                                             ( d2 ) N-M
  203: *
  204:       CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK,
  205:      $             D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
  206:       LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) )
  207: *
  208: *     Solve T22*y2 = d2 for y2
  209: *
  210:       IF( N.GT.M ) THEN
  211:          CALL ZTRTRS( 'Upper', 'No transpose', 'Non unit', N-M, 1,
  212:      $                B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO )
  213: *
  214:          IF( INFO.GT.0 ) THEN
  215:             INFO = 1
  216:             RETURN
  217:          END IF
  218: *
  219:          CALL ZCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 )
  220:       END IF
  221: *
  222: *     Set y1 = 0
  223: *
  224:       DO 10 I = 1, M + P - N
  225:          Y( I ) = CZERO
  226:    10 CONTINUE
  227: *
  228: *     Update d1 = d1 - T12*y2
  229: *
  230:       CALL ZGEMV( 'No transpose', M, N-M, -CONE, B( 1, M+P-N+1 ), LDB,
  231:      $            Y( M+P-N+1 ), 1, CONE, D, 1 )
  232: *
  233: *     Solve triangular system: R11*x = d1
  234: *
  235:       IF( M.GT.0 ) THEN
  236:          CALL ZTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA,
  237:      $                D, M, INFO )
  238: *
  239:          IF( INFO.GT.0 ) THEN
  240:             INFO = 2
  241:             RETURN
  242:          END IF
  243: *
  244: *        Copy D to X
  245: *
  246:          CALL ZCOPY( M, D, 1, X, 1 )
  247:       END IF
  248: *
  249: *     Backward transformation y = Z'*y
  250: *
  251:       CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP,
  252:      $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
  253:      $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
  254:       WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
  255: *
  256:       RETURN
  257: *
  258: *     End of ZGGGLM
  259: *
  260:       END