Dggglm Dgglse

Dggglm() Sub Dggglm ( N As  Long, M As  Long, P As  Long, A() As  Double, B() As  Double, D() As  Double, X() As  Double, Y() As  Double, Info As  Long  ) General Gauss-Markov linear model (GLM) problem Purpose This routine solves a general Gauss-Markov linear model (GLM) problem: minimize || y ||_2 subject to d = A*x + B*y x where A is an N x M matrix, B is an N x P matrix, and d is a given N vector. It is assumed that M <= N <= M + P, and rank(A) = M and rank(A B) = N Under these assumptions, the constrained equation is always consistent, and there is a unique solution x and a minimal 2-norm solution y, which is obtained using a generalized QR factorization of the matrices (A, B) given by A = Q*(R), B = Q*T*Z (0) In particular, if matrix B is square nonsingular, then the problem GLM is equivalent to the following weighted linear least squares problem minimize || inv(B)*(d - A*x) ||_2 x Parameters [in] N Number of rows of the matrices A and B. (N >= 0) (If N = 0, returns without computation) [in] M Number of columns of the matrix A. (0 <= M <= N) [in] P Number of columns of the matrix B. (N - M <= P) [in,out] A() Array A(LA1 - 1, LA2 - 1) (LA1 >= N, LA2 >= M) [in] N x M matrix A. [out] The upper triangular part of the array A() contains the M x M upper triangular matrix R. [in,out] B() Array B(LB1 - 1, LB2 - 1) (LB1 >= N, LB2 >= P) [in] N x P matrix B. [out] If N <= P, the upper triangle of the subarray B(0 to N-1, P-N to P-1) contains the N x N upper triangular matrix T. If N > P, the elements on and above the (N-P)-th subdiagonal contain the N x P upper trapezoidal matrix T. [in,out] D() Array D(LD - 1) (LD >= N) [in] D() is the left hand side of the GLM equation. [out] D() is destroyed. [out] X() Array X(LX - 1) (LX >= M) [out] Y() Array Y(LY - 1) (LY >= P) X() and Y() are the solutions of the GLM problem. [out] Info = 0: Successful exit. = -1: The argument N had an illegal value. (N < 0) = -2: The argument M had an illegal value. (M < 0 or M > N) = -3: The argument P had an illegal value. (P < 0 or P < N - M) = -4: The argument A() is invalid. = -5: The argument B() is invalid. = -6: The argument D() is invalid. = -7: The argument X() is invalid. = -8: The argument Y() is invalid. = 1: The upper triangular factor R associated with A in the generalized QR factorization of the pair (A, B) is singular, so that rank(A) < M. The least squares solution could not be computed. = 2: The bottom N-M x N-M part of the upper trapezoidal factor T associated with B in the generalized QR factorization of the pair (A, B) is singular, so that rank(A B) < N. The least squares solution could not be computed. Reference LAPACK Example Program Solve a general Gauss-Markov linear model (GLM) problem, i.e. find x which minimizes || y ||_2 subject to d = A*x + B*y, where ( -1.06 0.48 -0.04 ) A = ( -1.19 0.73 -0.24 ) ( 1.97 -0.89 0.56 ) ( 0.68 -0.53 0.08 ) ( 1 0 0 0 ) B = ( 0 1 0 0 ) ( 0 0 1 0 ) ( 0 0 0 1 ) ( 0.3884 ) d = ( 0.1120 ) ( -0.3644 ) ( -0.0002 ) Sub EX_Dggglm() Const M = 3, N = 4, P = 4 Dim A(N - 1, M - 1) As Double, B(N - 1, P - 1) As Double Dim C(M - 1) As Double, D(P - 1) As Double Dim X(M - 1) As Double, Y(P - 1) As Double Dim S As Double, Info As Long A(0, 0) = -1.06: A(0, 1) = 0.48: A(0, 2) = -0.04 A(1, 0) = -1.19: A(1, 1) = 0.73: A(1, 2) = -0.24 A(2, 0) = 1.97: A(2, 1) = -0.89: A(2, 2) = 0.56 A(3, 0) = 0.68: A(3, 1) = -0.53: A(3, 2) = 0.08 B(0, 0) = 1: B(0, 1) = 0: B(0, 2) = 0: B(0, 3) = 0 B(1, 0) = 0: B(1, 1) = 1: B(1, 2) = 0: B(1, 3) = 0 B(2, 0) = 0: B(2, 1) = 0: B(2, 2) = 1: B(2, 3) = 0 B(3, 0) = 0: B(3, 1) = 0: B(3, 2) = 0: B(3, 3) = 1 D(0) = 0.3884: D(1) = 0.112: D(2) = -0.3644: D(3) = -0.0002 Call Dggglm(N, M, P, A(), B(), D(), X(), Y(), Info) S = Dnrm2(P, Y(0)) Debug.Print "X =", X(0), X(1), X(2) Debug.Print "SumSq =", S, "Info =", Info End Sub Example Results X = -0.820000000000001 -0.940000000000001 0.740000000000001 SumSq = 1.45716771982052E-16 Info = 0