◆ Zgels()

Sub Zgels	(	Trans As	String,
		M As	Long,
		N As	Long,
		A() As	Complex,
		B() As	Complex,
		Info As	Long,
		Optional Nrhs As	Long = `1`
	)

Solution to overdetermined or underdetermined linear equations Ax = b for complex matrices (full rank)

Purpose

This routine solves overdetermined or underdetermined complex linear systems involving an M x N matrix A, or its conjugate-transpose, using a QR or LQ factorization of A. It is assumed that A has full rank.

The following options are provided:

If Trans = "N" and M >= N: find the least squares solution of an overdetermined system, i.e., solve the least squares problem
minimize || B - A*X ||.
If Trans = "N" and M < N: find the minimum norm solution of an underdetermined system A * X = B.
If Trans = "T" and M >= N: find the minimum norm solution of an underdetermined system A^H * X = B.
If Trans = "T" and M < N: find the least squares solution of an overdetermined system, i.e., solve the least squares problem
minimize || B - A^H*X ||.

Several right hand side vectors b and solution vectors x can be handled in a single call; they are stored as the columns of the M x Nrhs right hand side matrix B and the N x Nrhs solution matrix X.

Parameters

[in]	Trans	= "N": The linear system involves A. = "T": The linear system involves A^H.
[in]	M	Number of rows of the matrix A. (M >= 0) (If M = 0, returns zero vectors in B())
[in]	N	Number of columns of the matrix A. (N >= 0) (If N = 0, returns zero vectors in B())
[in,out]	A()	Array A(LA1 - 1, LA2 - 1) (LA1 >= M, LA2 >= N) [in] M x N matrix A. [out] M >= N: A() is overwritten by details of its QR factorization as returned by Zgeqrf. M < N: A() is overwritten by details of its LQ factorization as returned by Zgelqf.
[in,out]	B()	Array B(LB1 - 1, LB2 - 1) (LB1 >= max(M, N), LB2 >= Nrhs) (2D array) or B(LB - 1) (LB >= max(M, N), Nrhs = 1) (1D array) [in] Matrix B of right hand side vectors, stored columnwise; B is M x Nrhs if Trans = "N", or N x Nrhs if Trans = "T". [out] If Info = 0, B() is overwritten by the solution vectors, stored columnwise: Trans = "N" and M >= N: Rows 0 to N-1 of B() contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements N to M-1 in that column. Trans = "N" and M < N: Rows 0 to N-1 of B() contain the minimum norm solution vectors. Trans = "T" and M >= N: Rows 0 to M-1 of B() contain the minimum norm solution vectors. Trans = "T" and M < N: Rows 0 to M-1 of B() contain the least squares solution vectors. The residual sum of squares for the solution in each column is given by the sum of squares of elements M to N-1 in that column.
[out]	Info	= 0: Successful exit. = -1: The argument Trans had an illegal value. (Trans <> "T" nor "N") = -2: The argument M had an illegal value. (M < 0) = -3: The argument N had an illegal value. (N < 0) = -4: The argument A() is invalid. = -5: The argument B() is invalid. = -7: The argument Nrhs had an illegal value. (Nrhs < 0) = i > 0: The i-th diagonal element of the triangular factor of A is zero, so that A does not have full rank. The least squares solution could not be computed.
[in]	Nrhs	(Optional) Number of right hand sides, i.e., number of columns of the matrices B and X. (Nrhs >= 0) (If Nrhs = 0, returns without computation) (default = 1)

Reference: LAPACK

Example Program: Compute the least squares solution of the overdetermined linear equations Ax = b and its variance, 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 )

( 1.7126-0.6648i )

B = ( 0.8697+0.7604i )

( -2.1048-1.6171i )

( -0.9297+0.1252i )

Sub Ex_Zgels()

Const M = 4, N = 3

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

Dim 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) = Cmplx(1.7126, -0.6648): B(1) = Cmplx(0.8697, 0.7604)

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

Call Zgels("N", M, N, A(), B(), Info)

If Info <> 0 Then

Debug.Print "Error in Zgels: Info =", Info

Exit Sub

End If

Debug.Print "X ="

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

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

Call Zgecov(0, N, A(), Ci(), Info)

Debug.Print "Var ="

Debug.Print Creal(Ci(0)), Creal(Ci(1)), Creal(Ci(2))

Debug.Print "Info =", Info

End Sub

Example Results: X =

-0.819999999999999 -0.939999999999999 0.740000000000001 0.200000000000002

0.480000000000002 0.21

Var =

5.46169501938982 15.1880504061464 21.4241290120714

Info = 0