LAPACK 3.3.0

cgglse.f

Go to the documentation of this file.
00001       SUBROUTINE CGGLSE( 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       COMPLEX            A( LDA, * ), B( LDB, * ), C( * ), D( * ),
00014      $                   WORK( * ), X( * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  CGGLSE 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) COMPLEX 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) COMPLEX 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) COMPLEX 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) COMPLEX 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) COMPLEX array, dimension (N)
00079 *          On exit, X is the solution of the LSE problem.
00080 *
00081 *  WORK    (workspace/output) COMPLEX 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 *          CGEQRF, CGERQF, CUNMQR and CUNMRQ.
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       COMPLEX            CONE
00113       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+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           CAXPY, CCOPY, CGEMV, CGGRQF, CTRMV, CTRTRS,
00122      $                   CUNMQR, CUNMRQ, 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, 'CGEQRF', ' ', M, N, -1, -1 )
00158             NB2 = ILAENV( 1, 'CGERQF', ' ', M, N, -1, -1 )
00159             NB3 = ILAENV( 1, 'CUNMQR', ' ', M, N, P, -1 )
00160             NB4 = ILAENV( 1, 'CUNMRQ', ' ', 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( 'CGGLSE', -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 *     unitary.
00192 *
00193       CALL CGGRQF( 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 CUNMQR( 'Left', 'Conjugate Transpose', M, 1, MN, A, LDA,
00201      $             WORK( P+1 ), C, MAX( 1, M ), WORK( P+MN+1 ),
00202      $             LWORK-P-MN, INFO )
00203       LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
00204 *
00205 *     Solve T12*x2 = d for x2
00206 *
00207       IF( P.GT.0 ) THEN
00208          CALL CTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
00209      $                B( 1, N-P+1 ), LDB, D, P, INFO )
00210 *
00211          IF( INFO.GT.0 ) THEN
00212             INFO = 1
00213             RETURN
00214          END IF
00215 *
00216 *        Put the solution in X
00217 *
00218       CALL CCOPY( P, D, 1, X( N-P+1 ), 1 )
00219 *
00220 *        Update c1
00221 *
00222          CALL CGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), LDA,
00223      $               D, 1, CONE, C, 1 )
00224       END IF
00225 *
00226 *     Solve R11*x1 = c1 for x1
00227 *
00228       IF( N.GT.P ) THEN
00229          CALL CTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
00230      $                A, LDA, C, N-P, INFO )
00231 *
00232          IF( INFO.GT.0 ) THEN
00233             INFO = 2
00234             RETURN
00235          END IF
00236 *
00237 *        Put the solutions in X
00238 *
00239          CALL CCOPY( N-P, C, 1, X, 1 )
00240       END IF
00241 *
00242 *     Compute the residual vector:
00243 *
00244       IF( M.LT.N ) THEN
00245          NR = M + P - N
00246          IF( NR.GT.0 )
00247      $      CALL CGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, M+1 ),
00248      $                  LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 )
00249       ELSE
00250          NR = P
00251       END IF
00252       IF( NR.GT.0 ) THEN
00253          CALL CTRMV( 'Upper', 'No transpose', 'Non unit', NR,
00254      $               A( N-P+1, N-P+1 ), LDA, D, 1 )
00255          CALL CAXPY( NR, -CONE, D, 1, C( N-P+1 ), 1 )
00256       END IF
00257 *
00258 *     Backward transformation x = Q'*x
00259 *
00260       CALL CUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB,
00261      $             WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
00262       WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
00263 *
00264       RETURN
00265 *
00266 *     End of CGGLSE
00267 *
00268       END
 All Files Functions