version 1.3, 2010/08/06 15:28:49
|
version 1.11, 2012/08/22 09:48:27
|
Line 1
|
Line 1
|
|
*> \brief \b DTGSYL |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DTGSYL + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsyl.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsyl.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsyl.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, |
|
* LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, |
|
* IWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* CHARACTER TRANS |
|
* INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, |
|
* $ LWORK, M, N |
|
* DOUBLE PRECISION DIF, SCALE |
|
* .. |
|
* .. Array Arguments .. |
|
* INTEGER IWORK( * ) |
|
* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), |
|
* $ D( LDD, * ), E( LDE, * ), F( LDF, * ), |
|
* $ WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DTGSYL solves the generalized Sylvester equation: |
|
*> |
|
*> A * R - L * B = scale * C (1) |
|
*> D * R - L * E = scale * F |
|
*> |
|
*> where R and L are unknown m-by-n matrices, (A, D), (B, E) and |
|
*> (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n, |
|
*> respectively, with real entries. (A, D) and (B, E) must be in |
|
*> generalized (real) Schur canonical form, i.e. A, B are upper quasi |
|
*> triangular and D, E are upper triangular. |
|
*> |
|
*> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output |
|
*> scaling factor chosen to avoid overflow. |
|
*> |
|
*> In matrix notation (1) is equivalent to solve Zx = scale b, where |
|
*> Z is defined as |
|
*> |
|
*> Z = [ kron(In, A) -kron(B**T, Im) ] (2) |
|
*> [ kron(In, D) -kron(E**T, Im) ]. |
|
*> |
|
*> Here Ik is the identity matrix of size k and X**T is the transpose of |
|
*> X. kron(X, Y) is the Kronecker product between the matrices X and Y. |
|
*> |
|
*> If TRANS = 'T', DTGSYL solves the transposed system Z**T*y = scale*b, |
|
*> which is equivalent to solve for R and L in |
|
*> |
|
*> A**T * R + D**T * L = scale * C (3) |
|
*> R * B**T + L * E**T = scale * -F |
|
*> |
|
*> This case (TRANS = 'T') is used to compute an one-norm-based estimate |
|
*> of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D) |
|
*> and (B,E), using DLACON. |
|
*> |
|
*> If IJOB >= 1, DTGSYL computes a Frobenius norm-based estimate |
|
*> of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the |
|
*> reciprocal of the smallest singular value of Z. See [1-2] for more |
|
*> information. |
|
*> |
|
*> This is a level 3 BLAS algorithm. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] TRANS |
|
*> \verbatim |
|
*> TRANS is CHARACTER*1 |
|
*> = 'N', solve the generalized Sylvester equation (1). |
|
*> = 'T', solve the 'transposed' system (3). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] IJOB |
|
*> \verbatim |
|
*> IJOB is INTEGER |
|
*> Specifies what kind of functionality to be performed. |
|
*> =0: solve (1) only. |
|
*> =1: The functionality of 0 and 3. |
|
*> =2: The functionality of 0 and 4. |
|
*> =3: Only an estimate of Dif[(A,D), (B,E)] is computed. |
|
*> (look ahead strategy IJOB = 1 is used). |
|
*> =4: Only an estimate of Dif[(A,D), (B,E)] is computed. |
|
*> ( DGECON on sub-systems is used ). |
|
*> Not referenced if TRANS = 'T'. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] M |
|
*> \verbatim |
|
*> M is INTEGER |
|
*> The order of the matrices A and D, and the row dimension of |
|
*> the matrices C, F, R and L. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] N |
|
*> \verbatim |
|
*> N is INTEGER |
|
*> The order of the matrices B and E, and the column dimension |
|
*> of the matrices C, F, R and L. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] A |
|
*> \verbatim |
|
*> A is DOUBLE PRECISION array, dimension (LDA, M) |
|
*> The upper quasi triangular matrix A. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDA |
|
*> \verbatim |
|
*> LDA is INTEGER |
|
*> The leading dimension of the array A. LDA >= max(1, M). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] B |
|
*> \verbatim |
|
*> B is DOUBLE PRECISION array, dimension (LDB, N) |
|
*> The upper quasi triangular matrix B. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDB |
|
*> \verbatim |
|
*> LDB is INTEGER |
|
*> The leading dimension of the array B. LDB >= max(1, N). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] C |
|
*> \verbatim |
|
*> C is DOUBLE PRECISION array, dimension (LDC, N) |
|
*> On entry, C contains the right-hand-side of the first matrix |
|
*> equation in (1) or (3). |
|
*> On exit, if IJOB = 0, 1 or 2, C has been overwritten by |
|
*> the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R, |
|
*> the solution achieved during the computation of the |
|
*> Dif-estimate. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDC |
|
*> \verbatim |
|
*> LDC is INTEGER |
|
*> The leading dimension of the array C. LDC >= max(1, M). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] D |
|
*> \verbatim |
|
*> D is DOUBLE PRECISION array, dimension (LDD, M) |
|
*> The upper triangular matrix D. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDD |
|
*> \verbatim |
|
*> LDD is INTEGER |
|
*> The leading dimension of the array D. LDD >= max(1, M). |
|
*> \endverbatim |
|
*> |
|
*> \param[in] E |
|
*> \verbatim |
|
*> E is DOUBLE PRECISION array, dimension (LDE, N) |
|
*> The upper triangular matrix E. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDE |
|
*> \verbatim |
|
*> LDE is INTEGER |
|
*> The leading dimension of the array E. LDE >= max(1, N). |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] F |
|
*> \verbatim |
|
*> F is DOUBLE PRECISION array, dimension (LDF, N) |
|
*> On entry, F contains the right-hand-side of the second matrix |
|
*> equation in (1) or (3). |
|
*> On exit, if IJOB = 0, 1 or 2, F has been overwritten by |
|
*> the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L, |
|
*> the solution achieved during the computation of the |
|
*> Dif-estimate. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDF |
|
*> \verbatim |
|
*> LDF is INTEGER |
|
*> The leading dimension of the array F. LDF >= max(1, M). |
|
*> \endverbatim |
|
*> |
|
*> \param[out] DIF |
|
*> \verbatim |
|
*> DIF is DOUBLE PRECISION |
|
*> On exit DIF is the reciprocal of a lower bound of the |
|
*> reciprocal of the Dif-function, i.e. DIF is an upper bound of |
|
*> Dif[(A,D), (B,E)] = sigma_min(Z), where Z as in (2). |
|
*> IF IJOB = 0 or TRANS = 'T', DIF is not touched. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] SCALE |
|
*> \verbatim |
|
*> SCALE is DOUBLE PRECISION |
|
*> On exit SCALE is the scaling factor in (1) or (3). |
|
*> If 0 < SCALE < 1, C and F hold the solutions R and L, resp., |
|
*> to a slightly perturbed system but the input matrices A, B, D |
|
*> and E have not been changed. If SCALE = 0, C and F hold the |
|
*> solutions R and L, respectively, to the homogeneous system |
|
*> with C = F = 0. Normally, SCALE = 1. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
|
*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LWORK |
|
*> \verbatim |
|
*> LWORK is INTEGER |
|
*> The dimension of the array WORK. LWORK > = 1. |
|
*> If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N). |
|
*> |
|
*> If LWORK = -1, then a workspace query is assumed; the routine |
|
*> only calculates the optimal size of the WORK array, returns |
|
*> this value as the first entry of the WORK array, and no error |
|
*> message related to LWORK is issued by XERBLA. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IWORK |
|
*> \verbatim |
|
*> IWORK is INTEGER array, dimension (M+N+6) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> =0: successful exit |
|
*> <0: If INFO = -i, the i-th argument had an illegal value. |
|
*> >0: (A, D) and (B, E) have common or close eigenvalues. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date November 2011 |
|
* |
|
*> \ingroup doubleSYcomputational |
|
* |
|
*> \par Contributors: |
|
* ================== |
|
*> |
|
*> Bo Kagstrom and Peter Poromaa, Department of Computing Science, |
|
*> Umea University, S-901 87 Umea, Sweden. |
|
* |
|
*> \par References: |
|
* ================ |
|
*> |
|
*> \verbatim |
|
*> |
|
*> [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software |
|
*> for Solving the Generalized Sylvester Equation and Estimating the |
|
*> Separation between Regular Matrix Pairs, Report UMINF - 93.23, |
|
*> Department of Computing Science, Umea University, S-901 87 Umea, |
|
*> Sweden, December 1993, Revised April 1994, Also as LAPACK Working |
|
*> Note 75. To appear in ACM Trans. on Math. Software, Vol 22, |
|
*> No 1, 1996. |
|
*> |
|
*> [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester |
|
*> Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal. |
|
*> Appl., 15(4):1045-1060, 1994 |
|
*> |
|
*> [3] B. Kagstrom and L. Westin, Generalized Schur Methods with |
|
*> Condition Estimators for Solving the Generalized Sylvester |
|
*> Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, |
|
*> July 1989, pp 745-751. |
|
*> \endverbatim |
|
*> |
|
* ===================================================================== |
SUBROUTINE DTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, |
SUBROUTINE DTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, |
$ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, |
$ LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, |
$ IWORK, INFO ) |
$ IWORK, INFO ) |
* |
* |
* -- LAPACK routine (version 3.2) -- |
* -- LAPACK computational routine (version 3.4.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2006 |
* November 2011 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER TRANS |
CHARACTER TRANS |
Line 20
|
Line 317
|
$ WORK( * ) |
$ WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DTGSYL solves the generalized Sylvester equation: |
|
* |
|
* A * R - L * B = scale * C (1) |
|
* D * R - L * E = scale * F |
|
* |
|
* where R and L are unknown m-by-n matrices, (A, D), (B, E) and |
|
* (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n, |
|
* respectively, with real entries. (A, D) and (B, E) must be in |
|
* generalized (real) Schur canonical form, i.e. A, B are upper quasi |
|
* triangular and D, E are upper triangular. |
|
* |
|
* The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output |
|
* scaling factor chosen to avoid overflow. |
|
* |
|
* In matrix notation (1) is equivalent to solve Zx = scale b, where |
|
* Z is defined as |
|
* |
|
* Z = [ kron(In, A) -kron(B', Im) ] (2) |
|
* [ kron(In, D) -kron(E', Im) ]. |
|
* |
|
* Here Ik is the identity matrix of size k and X' is the transpose of |
|
* X. kron(X, Y) is the Kronecker product between the matrices X and Y. |
|
* |
|
* If TRANS = 'T', DTGSYL solves the transposed system Z'*y = scale*b, |
|
* which is equivalent to solve for R and L in |
|
* |
|
* A' * R + D' * L = scale * C (3) |
|
* R * B' + L * E' = scale * (-F) |
|
* |
|
* This case (TRANS = 'T') is used to compute an one-norm-based estimate |
|
* of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D) |
|
* and (B,E), using DLACON. |
|
* |
|
* If IJOB >= 1, DTGSYL computes a Frobenius norm-based estimate |
|
* of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the |
|
* reciprocal of the smallest singular value of Z. See [1-2] for more |
|
* information. |
|
* |
|
* This is a level 3 BLAS algorithm. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* TRANS (input) CHARACTER*1 |
|
* = 'N', solve the generalized Sylvester equation (1). |
|
* = 'T', solve the 'transposed' system (3). |
|
* |
|
* IJOB (input) INTEGER |
|
* Specifies what kind of functionality to be performed. |
|
* =0: solve (1) only. |
|
* =1: The functionality of 0 and 3. |
|
* =2: The functionality of 0 and 4. |
|
* =3: Only an estimate of Dif[(A,D), (B,E)] is computed. |
|
* (look ahead strategy IJOB = 1 is used). |
|
* =4: Only an estimate of Dif[(A,D), (B,E)] is computed. |
|
* ( DGECON on sub-systems is used ). |
|
* Not referenced if TRANS = 'T'. |
|
* |
|
* M (input) INTEGER |
|
* The order of the matrices A and D, and the row dimension of |
|
* the matrices C, F, R and L. |
|
* |
|
* N (input) INTEGER |
|
* The order of the matrices B and E, and the column dimension |
|
* of the matrices C, F, R and L. |
|
* |
|
* A (input) DOUBLE PRECISION array, dimension (LDA, M) |
|
* The upper quasi triangular matrix A. |
|
* |
|
* LDA (input) INTEGER |
|
* The leading dimension of the array A. LDA >= max(1, M). |
|
* |
|
* B (input) DOUBLE PRECISION array, dimension (LDB, N) |
|
* The upper quasi triangular matrix B. |
|
* |
|
* LDB (input) INTEGER |
|
* The leading dimension of the array B. LDB >= max(1, N). |
|
* |
|
* C (input/output) DOUBLE PRECISION array, dimension (LDC, N) |
|
* On entry, C contains the right-hand-side of the first matrix |
|
* equation in (1) or (3). |
|
* On exit, if IJOB = 0, 1 or 2, C has been overwritten by |
|
* the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R, |
|
* the solution achieved during the computation of the |
|
* Dif-estimate. |
|
* |
|
* LDC (input) INTEGER |
|
* The leading dimension of the array C. LDC >= max(1, M). |
|
* |
|
* D (input) DOUBLE PRECISION array, dimension (LDD, M) |
|
* The upper triangular matrix D. |
|
* |
|
* LDD (input) INTEGER |
|
* The leading dimension of the array D. LDD >= max(1, M). |
|
* |
|
* E (input) DOUBLE PRECISION array, dimension (LDE, N) |
|
* The upper triangular matrix E. |
|
* |
|
* LDE (input) INTEGER |
|
* The leading dimension of the array E. LDE >= max(1, N). |
|
* |
|
* F (input/output) DOUBLE PRECISION array, dimension (LDF, N) |
|
* On entry, F contains the right-hand-side of the second matrix |
|
* equation in (1) or (3). |
|
* On exit, if IJOB = 0, 1 or 2, F has been overwritten by |
|
* the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L, |
|
* the solution achieved during the computation of the |
|
* Dif-estimate. |
|
* |
|
* LDF (input) INTEGER |
|
* The leading dimension of the array F. LDF >= max(1, M). |
|
* |
|
* DIF (output) DOUBLE PRECISION |
|
* On exit DIF is the reciprocal of a lower bound of the |
|
* reciprocal of the Dif-function, i.e. DIF is an upper bound of |
|
* Dif[(A,D), (B,E)] = sigma_min(Z), where Z as in (2). |
|
* IF IJOB = 0 or TRANS = 'T', DIF is not touched. |
|
* |
|
* SCALE (output) DOUBLE PRECISION |
|
* On exit SCALE is the scaling factor in (1) or (3). |
|
* If 0 < SCALE < 1, C and F hold the solutions R and L, resp., |
|
* to a slightly perturbed system but the input matrices A, B, D |
|
* and E have not been changed. If SCALE = 0, C and F hold the |
|
* solutions R and L, respectively, to the homogeneous system |
|
* with C = F = 0. Normally, SCALE = 1. |
|
* |
|
* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) |
|
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. |
|
* |
|
* LWORK (input) INTEGER |
|
* The dimension of the array WORK. LWORK > = 1. |
|
* If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N). |
|
* |
|
* If LWORK = -1, then a workspace query is assumed; the routine |
|
* only calculates the optimal size of the WORK array, returns |
|
* this value as the first entry of the WORK array, and no error |
|
* message related to LWORK is issued by XERBLA. |
|
* |
|
* IWORK (workspace) INTEGER array, dimension (M+N+6) |
|
* |
|
* INFO (output) INTEGER |
|
* =0: successful exit |
|
* <0: If INFO = -i, the i-th argument had an illegal value. |
|
* >0: (A, D) and (B, E) have common or close eigenvalues. |
|
* |
|
* Further Details |
|
* =============== |
|
* |
|
* Based on contributions by |
|
* Bo Kagstrom and Peter Poromaa, Department of Computing Science, |
|
* Umea University, S-901 87 Umea, Sweden. |
|
* |
|
* [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software |
|
* for Solving the Generalized Sylvester Equation and Estimating the |
|
* Separation between Regular Matrix Pairs, Report UMINF - 93.23, |
|
* Department of Computing Science, Umea University, S-901 87 Umea, |
|
* Sweden, December 1993, Revised April 1994, Also as LAPACK Working |
|
* Note 75. To appear in ACM Trans. on Math. Software, Vol 22, |
|
* No 1, 1996. |
|
* |
|
* [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester |
|
* Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal. |
|
* Appl., 15(4):1045-1060, 1994 |
|
* |
|
* [3] B. Kagstrom and L. Westin, Generalized Schur Methods with |
|
* Condition Estimators for Solving the Generalized Sylvester |
|
* Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7, |
|
* July 1989, pp 745-751. |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* Replaced various illegal calls to DCOPY by calls to DLASET. |
* Replaced various illegal calls to DCOPY by calls to DLASET. |
* Sven Hammarling, 1/5/02. |
* Sven Hammarling, 1/5/02. |
Line 485
|
Line 610
|
ELSE |
ELSE |
* |
* |
* Solve transposed (I, J)-subsystem |
* Solve transposed (I, J)-subsystem |
* A(I, I)' * R(I, J) + D(I, I)' * L(I, J) = C(I, J) |
* A(I, I)**T * R(I, J) + D(I, I)**T * L(I, J) = C(I, J) |
* R(I, J) * B(J, J)' + L(I, J) * E(J, J)' = -F(I, J) |
* R(I, J) * B(J, J)**T + L(I, J) * E(J, J)**T = -F(I, J) |
* for I = 1,2,..., P; J = Q, Q-1,..., 1 |
* for I = 1,2,..., P; J = Q, Q-1,..., 1 |
* |
* |
SCALE = ONE |
SCALE = ONE |