version 1.6, 2010/08/13 21:04:15
|
version 1.15, 2017/06/17 10:54:30
|
Line 1
|
Line 1
|
|
*> \brief \b ZTGSYL |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download ZTGSYL + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsyl.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsyl.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsyl.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE ZTGSYL( 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( * ) |
|
* COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ), |
|
* $ D( LDD, * ), E( LDE, * ), F( LDF, * ), |
|
* $ WORK( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> ZTGSYL 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 complex entries. A, B, D and E are upper |
|
*> triangular (i.e., (A,D) and (B,E) in generalized Schur form). |
|
*> |
|
*> 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**H, Im) ] (2) |
|
*> [ kron(In, D) -kron(E**H, Im) ], |
|
*> |
|
*> Here Ix is the identity matrix of size x and X**H is the conjugate |
|
*> transpose of X. Kron(X, Y) is the Kronecker product between the |
|
*> matrices X and Y. |
|
*> |
|
*> If TRANS = 'C', y in the conjugate transposed system Z**H *y = scale*b |
|
*> is solved for, which is equivalent to solve for R and L in |
|
*> |
|
*> A**H * R + D**H * L = scale * C (3) |
|
*> R * B**H + L * E**H = scale * -F |
|
*> |
|
*> This case (TRANS = 'C') 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 ZLACON. |
|
*> |
|
*> If IJOB >= 1, ZTGSYL 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. |
|
*> |
|
*> This is a level-3 BLAS algorithm. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] TRANS |
|
*> \verbatim |
|
*> TRANS is CHARACTER*1 |
|
*> = 'N': solve the generalized sylvester equation (1). |
|
*> = 'C': solve the "conjugate 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 is used). |
|
*> =4: Only an estimate of Dif[(A,D), (B,E)] is computed. |
|
*> (ZGECON on sub-systems is used). |
|
*> Not referenced if TRANS = 'C'. |
|
*> \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 COMPLEX*16 array, dimension (LDA, M) |
|
*> The upper 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 COMPLEX*16 array, dimension (LDB, N) |
|
*> The upper 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 = 'C', DIF is not referenced. |
|
*> \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, R and L will |
|
*> hold the solutions to the homogenious system with C = F = 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is COMPLEX*16 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+2) |
|
*> \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 very close |
|
*> eigenvalues. |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date December 2016 |
|
* |
|
*> \ingroup complex16SYcomputational |
|
* |
|
*> \par Contributors: |
|
* ================== |
|
*> |
|
*> Bo Kagstrom and Peter Poromaa, Department of Computing Science, |
|
*> Umea University, S-901 87 Umea, Sweden. |
|
* |
|
*> \par References: |
|
* ================ |
|
*> |
|
*> [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. |
|
*> \n |
|
*> [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. |
|
*> \n |
|
*> [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. |
|
*> |
|
* ===================================================================== |
SUBROUTINE ZTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, |
SUBROUTINE ZTGSYL( 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.7.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..-- |
* January 2007 |
* December 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER TRANS |
CHARACTER TRANS |
Line 20
|
Line 313
|
$ WORK( * ) |
$ WORK( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* ZTGSYL 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 complex entries. A, B, D and E are upper |
|
* triangular (i.e., (A,D) and (B,E) in generalized Schur form). |
|
* |
|
* 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 Ix is the identity matrix of size x and X' is the conjugate |
|
* transpose of X. Kron(X, Y) is the Kronecker product between the |
|
* matrices X and Y. |
|
* |
|
* If TRANS = 'C', y in the conjugate transposed system Z'*y = scale*b |
|
* is solved for, 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 = 'C') 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 ZLACON. |
|
* |
|
* If IJOB >= 1, ZTGSYL 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. |
|
* |
|
* This is a level-3 BLAS algorithm. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* TRANS (input) CHARACTER*1 |
|
* = 'N': solve the generalized sylvester equation (1). |
|
* = 'C': solve the "conjugate 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 is used). |
|
* =4: Only an estimate of Dif[(A,D), (B,E)] is computed. |
|
* (ZGECON on sub-systems is used). |
|
* Not referenced if TRANS = 'C'. |
|
* |
|
* 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) COMPLEX*16 array, dimension (LDA, M) |
|
* The upper triangular matrix A. |
|
* |
|
* LDA (input) INTEGER |
|
* The leading dimension of the array A. LDA >= max(1, M). |
|
* |
|
* B (input) COMPLEX*16 array, dimension (LDB, N) |
|
* The upper triangular matrix B. |
|
* |
|
* LDB (input) INTEGER |
|
* The leading dimension of the array B. LDB >= max(1, N). |
|
* |
|
* C (input/output) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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) COMPLEX*16 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 = 'C', DIF is not referenced. |
|
* |
|
* 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, R and L will |
|
* hold the solutions to the homogenious system with C = F = 0. |
|
* |
|
* WORK (workspace/output) COMPLEX*16 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+2) |
|
* |
|
* 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 very 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 CCOPY by calls to CLASET. |
* Replaced various illegal calls to CCOPY by calls to CLASET. |
* Sven Hammarling, 1/5/02. |
* Sven Hammarling, 1/5/02. |
Line 492
|
Line 614
|
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)**H * R(I, J) + D(I, I)**H * 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) + L(I, J) * E(J, J) = -F(I, J) |
* for I = 1,2,..., P; J = Q, Q-1,..., 1 |
* for I = 1,2,..., P; J = Q, Q-1,..., 1 |
* |
* |