Purpose
To compute the following matrices -1 G = B*R *B', - -1 A = A - B*R *L', - -1 Q = Q - L*R *L', where A, B, Q, R, L, and G are N-by-N, N-by-M, N-by-N, M-by-M, N-by-M, and N-by-N matrices, respectively, with Q, R and G symmetric matrices. When R is well-conditioned with respect to inversion, standard algorithms for solving linear-quadratic optimization problems will then also solve optimization problems with coupling weighting matrix L. Moreover, a gain in efficiency is possible using matrix G in the deflating subspace algorithms (see SLICOT Library routine SB02OD).Specification
SUBROUTINE SB02MT( JOBG, JOBL, FACT, UPLO, N, M, A, LDA, B, LDB, $ Q, LDQ, R, LDR, L, LDL, IPIV, OUFACT, G, LDG, $ IWORK, DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER FACT, JOBG, JOBL, UPLO INTEGER INFO, LDA, LDB, LDG, LDL, LDQ, LDR, LDWORK, M, $ N, OUFACT C .. Array Arguments .. INTEGER IPIV(*), IWORK(*) DOUBLE PRECISION A(LDA,*), B(LDB,*), DWORK(*), G(LDG,*), $ L(LDL,*), Q(LDQ,*), R(LDR,*)Arguments
Mode Parameters
JOBG CHARACTER*1 Specifies whether or not the matrix G is to be computed, as follows: = 'G': Compute G; = 'N': Do not compute G. JOBL CHARACTER*1 Specifies whether or not the matrix L is zero, as follows: = 'Z': L is zero; = 'N': L is nonzero. FACT CHARACTER*1 Specifies how the matrix R is given (factored or not), as follows: = 'N': Array R contains the matrix R; = 'C': Array R contains the Cholesky factor of R; = 'U': Array R contains the symmetric indefinite UdU' or LdL' factorization of R. UPLO CHARACTER*1 Specifies which triangle of the matrices R and Q (if JOBL = 'N') is stored, as follows: = 'U': Upper triangle is stored; = 'L': Lower triangle is stored.Input/Output Parameters
N (input) INTEGER The order of the matrices A, Q, and G, and the number of rows of the matrices B and L. N >= 0. M (input) INTEGER The order of the matrix R, and the number of columns of the matrices B and L. M >= 0. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, if JOBL = 'N', the leading N-by-N part of this array must contain the matrix A. On exit, if JOBL = 'N', and INFO = 0, the leading N-by-N - -1 part of this array contains the matrix A = A - B*R L'. If JOBL = 'Z', this array is not referenced. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,N) if JOBL = 'N'; LDA >= 1 if JOBL = 'Z'. B (input/output) DOUBLE PRECISION array, dimension (LDB,M) On entry, the leading N-by-M part of this array must contain the matrix B. On exit, if OUFACT = 1, and INFO = 0, the leading N-by-M -1 part of this array contains the matrix B*chol(R) . On exit, B is unchanged if OUFACT = 2 (hence also when FACT = 'U'). LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N). Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) On entry, if JOBL = 'N', the leading N-by-N upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array must contain the upper triangular part or lower triangular part, respectively, of the symmetric matrix Q. The stricly lower triangular part (if UPLO = 'U') or stricly upper triangular part (if UPLO = 'L') is not referenced. On exit, if JOBL = 'N' and INFO = 0, the leading N-by-N upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array contains the upper triangular part or lower triangular part, respectively, of - -1 the symmetric matrix Q = Q - L*R *L'. If JOBL = 'Z', this array is not referenced. LDQ INTEGER The leading dimension of array Q. LDQ >= MAX(1,N) if JOBL = 'N'; LDQ >= 1 if JOBL = 'Z'. R (input/output) DOUBLE PRECISION array, dimension (LDR,M) On entry, if FACT = 'N', the leading M-by-M upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array must contain the upper triangular part or lower triangular part, respectively, of the symmetric input weighting matrix R. On entry, if FACT = 'C', the leading M-by-M upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array must contain the Cholesky factor of the positive definite input weighting matrix R (as produced by LAPACK routine DPOTRF). On entry, if FACT = 'U', the leading M-by-M upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array must contain the factors of the UdU' or LdL' factorization, respectively, of the symmetric indefinite input weighting matrix R (as produced by LAPACK routine DSYTRF). If FACT = 'N', the stricly lower triangular part (if UPLO = 'U') or stricly upper triangular part (if UPLO = 'L') of this array is used as workspace. On exit, if OUFACT = 1, and INFO = 0 (or INFO = M+1), the leading M-by-M upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array contains the Cholesky factor of the given input weighting matrix. On exit, if OUFACT = 2, and INFO = 0 (or INFO = M+1), the leading M-by-M upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array contains the factors of the UdU' or LdL' factorization, respectively, of the given input weighting matrix. On exit R is unchanged if FACT = 'C' or 'U'. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,M). L (input/output) DOUBLE PRECISION array, dimension (LDL,M) On entry, if JOBL = 'N', the leading N-by-M part of this array must contain the matrix L. On exit, if JOBL = 'N', OUFACT = 1, and INFO = 0, the leading N-by-M part of this array contains the matrix -1 L*chol(R) . On exit, L is unchanged if OUFACT = 2 (hence also when FACT = 'U'). L is not referenced if JOBL = 'Z'. LDL INTEGER The leading dimension of array L. LDL >= MAX(1,N) if JOBL = 'N'; LDL >= 1 if JOBL = 'Z'. IPIV (input/output) INTEGER array, dimension (M) On entry, if FACT = 'U', this array must contain details of the interchanges performed and the block structure of the d factor in the UdU' or LdL' factorization of matrix R (as produced by LAPACK routine DSYTRF). On exit, if OUFACT = 2, this array contains details of the interchanges performed and the block structure of the d factor in the UdU' or LdL' factorization of matrix R, as produced by LAPACK routine DSYTRF. This array is not referenced if FACT = 'C'. OUFACT (output) INTEGER Information about the factorization finally used. OUFACT = 1: Cholesky factorization of R has been used; OUFACT = 2: UdU' (if UPLO = 'U') or LdL' (if UPLO = 'L') factorization of R has been used. G (output) DOUBLE PRECISION array, dimension (LDG,N) If JOBG = 'G', and INFO = 0, the leading N-by-N upper triangular part (if UPLO = 'U') or lower triangular part (if UPLO = 'L') of this array contains the upper triangular part (if UPLO = 'U') or lower triangular part -1 (if UPLO = 'L'), respectively, of the matrix G = B*R B'. If JOBG = 'N', this array is not referenced. LDG INTEGER The leading dimension of array G. LDG >= MAX(1,N) if JOBG = 'G', LDG >= 1 if JOBG = 'N'.Workspace
IWORK INTEGER array, dimension (M) DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal value of LDWORK; if FACT = 'N', DWORK(2) contains the reciprocal condition number of the given matrix R. LDWORK INTEGER The length of the array DWORK. LDWORK >= 1 if FACT = 'C'; LDWORK >= MAX(2,3*M,N*M) if FACT = 'N'; LDWORK >= MAX(1,N*M) if FACT = 'U'. For optimum performance LDWORK should be larger than 3*M, if FACT = 'N'. The N*M workspace is not needed for FACT = 'N', if matrix R is positive definite.Error Indicator
INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; = i: if the i-th element (1 <= i <= M) of the d factor is exactly zero; the UdU' (or LdL') factorization has been completed, but the block diagonal matrix d is exactly singular; = M+1: if the matrix R is numerically singular.Method
- - The matrices G, and/or A and Q are evaluated using the given or computed symmetric factorization of R.Numerical Aspects
The routine should not be used when R is ill-conditioned.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None
Click here to get a compressed (gzip) tar file containing the source code of the routine, the example program, data, documentation, and related files.
Return to index