Linear Equality Constrained Least Squares Problem

The linear equality constrained least squares (LSE) problem is

$$\min_x \Vert Ax - b \Vert \quad \mbox{subject to} \quad Bx = d$$

where A is an m-by-n matrix, B is a p-by-n matrix, b is an m vector, and d is a p vector, with $p \leq n \leq p+m$.

The LSE problem is solved by the driver routine xGGLSE (see section 4.6). Let $\widehat {x}$ be the value of x computed by xGGLSE. The approximate error bound

$$\begin{eqnarray*} % latex2html id marker 10929\frac{ \Vert x - \widehat {x} \... ... ...$$

can be computed by the following code fragment

      EPSMCH = SLAMCH( 'E' )
*     Get the 2-norm of the right hand side C
      CNORM = SNRM2( M, C, 1 )
      print*,'CNORM = ',CNORM
*     Solve the least squares problem with equality constraints
      CALL SGGLSE( M, N, P, A, LDA, B, LDA, C, D, Xc, WORK, LWORK, IWORK, INFO )
*     Get the Frobenius norm of A and B
      ANORM = SLANTR( 'F', 'U', 'N', N, N, A, LDA, WORK )
      BNORM = SLANTR( 'F', 'U', 'N', P, P, B( 1, N-P+1 ), LDA, WORK )
      MN = MIN( M, N )
      IF( N.EQ.P ) THEN
         APPSNM = ZERO
         RNORM = SLANTR( '1', 'U', 'N', N, N, B, LDB, WORK(P+MN+N+1) )
         CALL STRCON( '1', 'U', 'N', N, B, LDB, RCOND, WORK( P+MN+N+1 ),
     $                IWORK, INFO )
         BAPSNM = ONE/ (RCOND * RNORM )
      ELSE
*        Estimate norm of (AP)^+
         RNORM = SLANTR( '1', 'U', 'N', N-P, N-P, A, LDA, WORK(P+MN+1) )
         CALL STRCON( '1', 'U', 'N', N-P, A, LDA, RCOND, WORK( P+MN+1 ),
     $                IWORK, INFO )
         APPSNM = ONE/ (RCOND * RNORM )
*        Estimate norm of B^+_A
         KASE = 0
         CALL SLACON( P, WORK( P+MN+1 ), WORK( P+MN+N+1 ), IWORK, EST, KASE )
   30    CONTINUE
            CALL STRSV( 'Upper', 'No trans', 'Non unit', P, B( 1, N-P+1 ),
     $                  LDB, WORK( P+MN+N+1 ), 1 )
            CALL SGEMV( 'No trans', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
     $                  WORK( P+MN+N+1 ), 1, ZERO, WORK( P+MN+P+1 ), 1 )
            CALL STRSV( 'Upper', 'No transpose', 'Non unit', N-P, A, LDA,
     $                  WORK( P+MN+P+1 ), 1 )
            DO I = 1, P
               WORK( P+MN+I ) = WORK( P+MN+N+I )
            END DO
            CALL SLACON( N, WORK( P+MN+N+1 ), WORK( P+MN+1 ), IWORK, EST, KASE )
*
            IF( KASE.EQ.0 ) GOTO 40
            DO I = 1, P
               WORK( P+MN+N+I ) = WORK( MN+N+I )
            END DO
            CALL STRSV( 'Upper', 'Trans', 'Non unit', N-P, A, LDA,
     $                  WORK( P+MN+1 ), 1 )
            CALL SGEMV( 'Trans', N-P, P, -ONE, A( 1, N-P+1 ), LDA,
     $                  WORK( P+MN+1 ), 1, ONE, WORK( P+MN+N+1 ), 1 )
            CALL STRSV( 'Upper', 'Trans', 'Non unit', P, B( 1, N-P+1 ),
     $                  LDB, WORK( P+MN+N+1 ), 1 )
            CALL SLACON( P, WORK( P+MN+1 ), WORK( P+MN+N+1 ), IWORK, EST, KASE )
*
            IF( KASE.EQ.0 ) GOTO 40
         GOTO 30
   40    CONTINUE
         BAPSNM = EST
*
      END IF
*     Estimate norm of A*B^+_A
      IF( P+M.EQ.N ) THEN
         EST = ZERO
      ELSE
         R22RS = MIN( P, M-N+P )
         KASE = 0
         CALL SLACON( P, WORK( P+MN+P+1 ), WORK( P+MN+1 ), IWORK, EST, KASE )
   50    CONTINUE
            CALL STRSV( 'Upper', 'No trans', 'Non unit', P, B( 1, N-P+1 ),
     $                  LDB, WORK( P+MN+1 ), 1 )
            DO I = 1, R22RS
               WORK( P+MN+P+I ) = WORK( P+MN+I )
            END DO
            CALL STRMV( 'Upper', 'No trans', 'Non unit', R22RS,
     $                  A( N-P+1, N-P+1 ), LDA, WORK( P+MN+P+1 ), 1 )
            IF( M.LT.N ) THEN
               CALL SGEMV( 'No trans', R22RS, N-M, ONE, A( N-P+1, M+1 ), LDA,
     $                     WORK( P+MN+R22RS+1 ), 1, ONE, WORK( P+MN+P+1 ), 1 )
            END IF
            CALL SLACON( R22RS, WORK( P+MN+1 ), WORK( P+MN+P+1 ), IWORK, EST,
     $                   KASE  )
*
            IF( KASE.EQ.0 ) GOTO 60
            DO I = 1, R22RS
               WORK( P+MN+I ) = WORK( P+MN+P+I )
            END DO
            CALL STRMV( 'Upper', 'Trans', 'Non Unit', R22RS,
     $                  A( N-P+1, N-P+1 ), LDA, WORK( P+MN+1 ), 1 )
            IF( M.LT.N ) THEN
               CALL SGEMV( 'Trans', R22RS, N-M, ONE, A( N-P+1, M+1 ), LDA,
     $                     WORK( P+MN+P+1 ), 1, ZERO, WORK( P+MN+R22RS+1 ), 1 )
            END IF
            CALL STRSV( 'Upper', 'Trans', 'Non unit', P, B( 1, N-P+1 ), LDB,
     $                  WORK( P+MN+1 ), 1 )
            CALL SLACON( P, WORK( P+MN+P+1 ), WORK( P+MN+1 ), IWORK, EST, KASE )
*
            IF( KASE.EQ.0 ) GOTO 60
            GOTO 50
   60    CONTINUE
      END IF
      ABAPSN = EST
*     Get the 2-norm of Xc
      XNORM = SNRM2( N, Xc, 1 )
      IF( APPSNM.EQ.0.0E0 ) THEN
*        B is square and nonsingular
         ERRBD = EPSMCH*BNORM*BAPSNM
      ELSE
*        Get the 2-norm of the residual A*Xc - C
         RNORM = SNRM2( M-N+P, C( N-P+1 ), 1 )
*        Get the 2-norm of Xc
         XNORM = SNRM2( N, Xc, 1 )
*        Get the condition numbers
         CNDBA = BNORM*BAPSNM
         CNDAB = ANORM*APPSNM
*        Get the approximate error bound
         ERRBD = EPSMCH*( (1.0E0 + CNORM/(ANORM*XNORM))*CNDAB +
     $               RNORM/(ANORM*XNORM)*(1.0E0 + BNORM*ABAPSN/ANORM)*
     $               (CNDAB*CNDAB) + 2.0E0*CNDBA )
      END IF

For example, if ${\tt SLAMCH('E')} = 2^{-24} = 5.961 \cdot 10^{-8}$,

$$A = \left( \begin{tabular}{rrrr} 1 & 1 & 1 & 1 \\ 1 & 3 & ... ...left( \begin{array}{r} 1 \\ 3 \\ -1 \end{array} \right),$$

then (to 7 decimal places),

$$\widehat {x} = \left( \begin{array}{r} 0.5000000 \\ -0.5000001 \\ 1.4999999 \\ 0.4999998 \end{array} \right).$$

The computed error bound ${\tt ERRBD} = 5.7\cdot10^{-7}$, where ${\tt CNDAB} = 2.09$, ${\tt CNDBA} = 3.12$. The true error $= 1.2\cdot10^{-7}$. The exact solution is x = [0.5, -0.5, 1.5, 0.5]^T.