Annotation of rpl/lapack/lapack/zgglse.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK,
! 2: $ INFO )
! 3: *
! 4: * -- LAPACK driver routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * February 2007
! 8: *
! 9: * .. Scalar Arguments ..
! 10: INTEGER INFO, LDA, LDB, LWORK, M, N, P
! 11: * ..
! 12: * .. Array Arguments ..
! 13: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( * ), D( * ),
! 14: $ WORK( * ), X( * )
! 15: * ..
! 16: *
! 17: * Purpose
! 18: * =======
! 19: *
! 20: * ZGGLSE solves the linear equality-constrained least squares (LSE)
! 21: * problem:
! 22: *
! 23: * minimize || c - A*x ||_2 subject to B*x = d
! 24: *
! 25: * where A is an M-by-N matrix, B is a P-by-N matrix, c is a given
! 26: * M-vector, and d is a given P-vector. It is assumed that
! 27: * P <= N <= M+P, and
! 28: *
! 29: * rank(B) = P and rank( ( A ) ) = N.
! 30: * ( ( B ) )
! 31: *
! 32: * These conditions ensure that the LSE problem has a unique solution,
! 33: * which is obtained using a generalized RQ factorization of the
! 34: * matrices (B, A) given by
! 35: *
! 36: * B = (0 R)*Q, A = Z*T*Q.
! 37: *
! 38: * Arguments
! 39: * =========
! 40: *
! 41: * M (input) INTEGER
! 42: * The number of rows of the matrix A. M >= 0.
! 43: *
! 44: * N (input) INTEGER
! 45: * The number of columns of the matrices A and B. N >= 0.
! 46: *
! 47: * P (input) INTEGER
! 48: * The number of rows of the matrix B. 0 <= P <= N <= M+P.
! 49: *
! 50: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
! 51: * On entry, the M-by-N matrix A.
! 52: * On exit, the elements on and above the diagonal of the array
! 53: * contain the min(M,N)-by-N upper trapezoidal matrix T.
! 54: *
! 55: * LDA (input) INTEGER
! 56: * The leading dimension of the array A. LDA >= max(1,M).
! 57: *
! 58: * B (input/output) COMPLEX*16 array, dimension (LDB,N)
! 59: * On entry, the P-by-N matrix B.
! 60: * On exit, the upper triangle of the subarray B(1:P,N-P+1:N)
! 61: * contains the P-by-P upper triangular matrix R.
! 62: *
! 63: * LDB (input) INTEGER
! 64: * The leading dimension of the array B. LDB >= max(1,P).
! 65: *
! 66: * C (input/output) COMPLEX*16 array, dimension (M)
! 67: * On entry, C contains the right hand side vector for the
! 68: * least squares part of the LSE problem.
! 69: * On exit, the residual sum of squares for the solution
! 70: * is given by the sum of squares of elements N-P+1 to M of
! 71: * vector C.
! 72: *
! 73: * D (input/output) COMPLEX*16 array, dimension (P)
! 74: * On entry, D contains the right hand side vector for the
! 75: * constrained equation.
! 76: * On exit, D is destroyed.
! 77: *
! 78: * X (output) COMPLEX*16 array, dimension (N)
! 79: * On exit, X is the solution of the LSE problem.
! 80: *
! 81: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
! 82: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 83: *
! 84: * LWORK (input) INTEGER
! 85: * The dimension of the array WORK. LWORK >= max(1,M+N+P).
! 86: * For optimum performance LWORK >= P+min(M,N)+max(M,N)*NB,
! 87: * where NB is an upper bound for the optimal blocksizes for
! 88: * ZGEQRF, CGERQF, ZUNMQR and CUNMRQ.
! 89: *
! 90: * If LWORK = -1, then a workspace query is assumed; the routine
! 91: * only calculates the optimal size of the WORK array, returns
! 92: * this value as the first entry of the WORK array, and no error
! 93: * message related to LWORK is issued by XERBLA.
! 94: *
! 95: * INFO (output) INTEGER
! 96: * = 0: successful exit.
! 97: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 98: * = 1: the upper triangular factor R associated with B in the
! 99: * generalized RQ factorization of the pair (B, A) is
! 100: * singular, so that rank(B) < P; the least squares
! 101: * solution could not be computed.
! 102: * = 2: the (N-P) by (N-P) part of the upper trapezoidal factor
! 103: * T associated with A in the generalized RQ factorization
! 104: * of the pair (B, A) is singular, so that
! 105: * rank( (A) ) < N; the least squares solution could not
! 106: * ( (B) )
! 107: * be computed.
! 108: *
! 109: * =====================================================================
! 110: *
! 111: * .. Parameters ..
! 112: COMPLEX*16 CONE
! 113: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
! 114: * ..
! 115: * .. Local Scalars ..
! 116: LOGICAL LQUERY
! 117: INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3,
! 118: $ NB4, NR
! 119: * ..
! 120: * .. External Subroutines ..
! 121: EXTERNAL XERBLA, ZAXPY, ZCOPY, ZGEMV, ZGGRQF, ZTRMV,
! 122: $ ZTRTRS, ZUNMQR, ZUNMRQ
! 123: * ..
! 124: * .. External Functions ..
! 125: INTEGER ILAENV
! 126: EXTERNAL ILAENV
! 127: * ..
! 128: * .. Intrinsic Functions ..
! 129: INTRINSIC INT, MAX, MIN
! 130: * ..
! 131: * .. Executable Statements ..
! 132: *
! 133: * Test the input parameters
! 134: *
! 135: INFO = 0
! 136: MN = MIN( M, N )
! 137: LQUERY = ( LWORK.EQ.-1 )
! 138: IF( M.LT.0 ) THEN
! 139: INFO = -1
! 140: ELSE IF( N.LT.0 ) THEN
! 141: INFO = -2
! 142: ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
! 143: INFO = -3
! 144: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 145: INFO = -5
! 146: ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
! 147: INFO = -7
! 148: END IF
! 149: *
! 150: * Calculate workspace
! 151: *
! 152: IF( INFO.EQ.0) THEN
! 153: IF( N.EQ.0 ) THEN
! 154: LWKMIN = 1
! 155: LWKOPT = 1
! 156: ELSE
! 157: NB1 = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
! 158: NB2 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 )
! 159: NB3 = ILAENV( 1, 'ZUNMQR', ' ', M, N, P, -1 )
! 160: NB4 = ILAENV( 1, 'ZUNMRQ', ' ', M, N, P, -1 )
! 161: NB = MAX( NB1, NB2, NB3, NB4 )
! 162: LWKMIN = M + N + P
! 163: LWKOPT = P + MN + MAX( M, N )*NB
! 164: END IF
! 165: WORK( 1 ) = LWKOPT
! 166: *
! 167: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
! 168: INFO = -12
! 169: END IF
! 170: END IF
! 171: *
! 172: IF( INFO.NE.0 ) THEN
! 173: CALL XERBLA( 'ZGGLSE', -INFO )
! 174: RETURN
! 175: ELSE IF( LQUERY ) THEN
! 176: RETURN
! 177: END IF
! 178: *
! 179: * Quick return if possible
! 180: *
! 181: IF( N.EQ.0 )
! 182: $ RETURN
! 183: *
! 184: * Compute the GRQ factorization of matrices B and A:
! 185: *
! 186: * B*Q' = ( 0 T12 ) P Z'*A*Q' = ( R11 R12 ) N-P
! 187: * N-P P ( 0 R22 ) M+P-N
! 188: * N-P P
! 189: *
! 190: * where T12 and R11 are upper triangular, and Q and Z are
! 191: * unitary.
! 192: *
! 193: CALL ZGGRQF( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
! 194: $ WORK( P+MN+1 ), LWORK-P-MN, INFO )
! 195: LOPT = WORK( P+MN+1 )
! 196: *
! 197: * Update c = Z'*c = ( c1 ) N-P
! 198: * ( c2 ) M+P-N
! 199: *
! 200: CALL ZUNMQR( 'Left', 'Conjugate Transpose', M, 1, MN, A, LDA,
! 201: $ WORK( P+1 ), C, MAX( 1, M ), WORK( P+MN+1 ),
! 202: $ LWORK-P-MN, INFO )
! 203: LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
! 204: *
! 205: * Solve T12*x2 = d for x2
! 206: *
! 207: IF( P.GT.0 ) THEN
! 208: CALL ZTRTRS( 'Upper', 'No transpose', 'Non-unit', P, 1,
! 209: $ B( 1, N-P+1 ), LDB, D, P, INFO )
! 210: *
! 211: IF( INFO.GT.0 ) THEN
! 212: INFO = 1
! 213: RETURN
! 214: END IF
! 215: *
! 216: * Put the solution in X
! 217: *
! 218: CALL ZCOPY( P, D, 1, X( N-P+1 ), 1 )
! 219: *
! 220: * Update c1
! 221: *
! 222: CALL ZGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), LDA,
! 223: $ D, 1, CONE, C, 1 )
! 224: END IF
! 225: *
! 226: * Solve R11*x1 = c1 for x1
! 227: *
! 228: IF( N.GT.P ) THEN
! 229: CALL ZTRTRS( 'Upper', 'No transpose', 'Non-unit', N-P, 1,
! 230: $ A, LDA, C, N-P, INFO )
! 231: *
! 232: IF( INFO.GT.0 ) THEN
! 233: INFO = 2
! 234: RETURN
! 235: END IF
! 236: *
! 237: * Put the solutions in X
! 238: *
! 239: CALL ZCOPY( N-P, C, 1, X, 1 )
! 240: END IF
! 241: *
! 242: * Compute the residual vector:
! 243: *
! 244: IF( M.LT.N ) THEN
! 245: NR = M + P - N
! 246: IF( NR.GT.0 )
! 247: $ CALL ZGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, M+1 ),
! 248: $ LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 )
! 249: ELSE
! 250: NR = P
! 251: END IF
! 252: IF( NR.GT.0 ) THEN
! 253: CALL ZTRMV( 'Upper', 'No transpose', 'Non unit', NR,
! 254: $ A( N-P+1, N-P+1 ), LDA, D, 1 )
! 255: CALL ZAXPY( NR, -CONE, D, 1, C( N-P+1 ), 1 )
! 256: END IF
! 257: *
! 258: * Backward transformation x = Q'*x
! 259: *
! 260: CALL ZUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB,
! 261: $ WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
! 262: WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) )
! 263: *
! 264: RETURN
! 265: *
! 266: * End of ZGGLSE
! 267: *
! 268: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgglse.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack