◆ Zgglse()

Sub Zgglse	(	M As	Long,
		N As	Long,
		P As	Long,
		A() As	Complex,
		B() As	Complex,
		C() As	Complex,
		D() As	Complex,
		X() As	Complex,
		Info As	Long
	)

Linear equality-constrained least squares (LSE) problem of complex matrices

Purpose: This routine solves the linear equality-constrained least squares (LSE) problem:
minimize || c - Ax ||_2 subject to B*x = d

where A is an M x N matrix, B is a P x N matrix, c is a given M-vector, and d is a given P-vector. It is assumed that P <= N <= M + P, and
rank(B) = P and rank( (A) ) = N

( (B) )

These conditions ensure that the LSE problem has a unique solution, which is obtained using a generalized RQ factorization of the matrices (B, A) given by
B = (0 R)*Q, A = Z*T*Q

Parameters

[in]	M	Number of rows of the matrix A. (M >= 0)
[in]	N	Number of columns of the matrices A and B. (N >= 0) (If N = 0, returns without
[in]	P	Number of rows of the matrix B. (0 <= P <= N <= M + P)
[in,out]	A()	Array A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] M x N matrix A [out] The elements on and above the diagonal of the array contain the min(M, N) x N upper trapezoidal matrix T.
[in,out]	B()	Array B(LB1 - 1, LB2 - 1) (LB1 >= P, LB2 >= N) [in] P x N matrix B [out] The upper triangle of the subarray B(0 to P-1, N-P to N-1) contains the P x P upper triangular matrix R.
[in,out]	C()	Array C(LC - 1) (LC >= M) [in] The right hand side vector for the least squares part of the LSE problem. [out] The residual sum of squares for the solution is given by the sum of squares of elements N-P to M-1 of vector C.
[in,out]	D()	Array D(LD - 1) (LD >= P) [in] The right hand side vector for the constrained equation. [out] D() is destroyed.
[out]	X()	Array X(LX - 1) (LX >= N) The solution of the LSE problem.
[out]	Info	= 0: Successful exit. = -1: The argument M had an illegal value. (M < 0) = -2: The argument N had an illegal value. (N < 0) = -3: The argument P had an illegal value. (P < 0 or P > N or P < N - M) = -4: The argument A() is invalid. = -5: The argument B() is invalid. = -6: The argument C() is invalid. = -7: The argument D() is invalid. = -8: The argument X() is invalid. = 1: The upper triangular factor R associated with B in the generalized RQ factorization of the pair (B, A) is singular, so that rank(B) < P. The least squares solution could not be computed. = 2: The N-P x N-P part of the upper trapezoidal factor T associated with A in the generalized RQ factorization of the pair (B, A) is singular, so that rank((A^T B^T)^T) < N. The least squares solution could not be computed.

Reference: LAPACK

Example Program: Solve the linear equality-constrained least squares (LSE) problem, i.e. minimize || c - Ax ||_2 subject to B*x = d, where
( -0.82+0.83i 0.18-0.94i -0.18-0.12i )

A = ( -0.76-0.24i 0.57-0.16i -0.08-0.27i )

( 1.90+0.26i -0.98+0.54i 0.21+0.28i )

( 0.50-0.30i -0.31+0.37i 0.22+0.19i )

( 0.57-0.91i -0.28-0.45i 0.25+0.91i )

B = ( 0.83+0.46i 0.63-0.19i -0.69+0.09i )

( 0.24-1.33i -0.56-0.67i 0.90+1.25i )

( 1.7126-0.6648i )

c = ( 0.8697+0.7604i )

( -2.1048-1.6171i )

( -0.9297+0.1252i )

( -1.5111+0.3107i )

d = ( -0.0941-1.2737i )

( -1.5579+1.0462i )

Sub Ex_Zgglse()

Const M = 4, N = 3, P = 3

Dim A(M - 1, N - 1) As Complex, B(P - 1, N - 1) As Complex

Dim C(M - 1) As Complex, D(P - 1) As Complex, X(N - 1) As Complex

Dim S As Double, Info As Long

A(0, 0) = Cmplx(-0.82, 0.83): A(0, 1) = Cmplx(0.18, -0.94): A(0, 2) = Cmplx(-0.18, -0.12)

A(1, 0) = Cmplx(-0.76, -0.24): A(1, 1) = Cmplx(0.57, -0.16): A(1, 2) = Cmplx(-0.08, -0.27)

A(2, 0) = Cmplx(1.9, 0.26): A(2, 1) = Cmplx(-0.98, 0.54): A(2, 2) = Cmplx(0.21, 0.28)

A(3, 0) = Cmplx(0.5, -0.3): A(3, 1) = Cmplx(-0.31, 0.37): A(3, 2) = Cmplx(0.22, 0.19)

B(0, 0) = Cmplx(0.57, -0.91): B(0, 1) = Cmplx(-0.28, -0.45): B(0, 2) = Cmplx(0.25, 0.91)

B(1, 0) = Cmplx(0.83, 0.46): B(1, 1) = Cmplx(0.63, -0.19): B(1, 2) = Cmplx(-0.69, 0.09)

B(2, 0) = Cmplx(0.24, -1.33): B(2, 1) = Cmplx(-0.56, -0.67): B(2, 2) = Cmplx(0.9, 1.25)

C(0) = Cmplx(1.7126, -0.6648): C(1) = Cmplx(0.8697, 0.7604)

C(2) = Cmplx(-2.1048, -1.6171): C(3) = Cmplx(-0.9297, 0.1252)

D(0) = Cmplx(-1.5111, 0.3107): D(1) = Cmplx(-0.0941, -1.2737)

D(2) = Cmplx(-1.5579, 1.0462)

Call Zgglse(M, N, P, A(), B(), C(), D(), X(), Info)

S = Dznrm2(M - N + P, C(N - P))

Debug.Print "X ="

Debug.Print Creal(X(0)), Cimag(X(0)), Creal(X(1)), Cimag(X(1))

Debug.Print Creal(X(2)), Cimag(X(2))

Debug.Print "SumSq =", S, "Info =", Info

End Sub

Example Results: X =

-0.82 -0.940000000000001 0.739999999999999 0.2

0.48 0.209999999999999

SumSq = 5.2372425074627E-15 Info = 0