Diff of /src/SB03OU.f [000000] .. [6a0e2a]  Maximize  Restore

Switch to unified view

a b/src/SB03OU.f
1
      SUBROUTINE SB03OU( DISCR, LTRANS, N, M, A, LDA, B, LDB, TAU, U,
2
     $                   LDU, SCALE, DWORK, LDWORK, INFO )
3
C
4
C     SLICOT RELEASE 5.0.
5
C
6
C     Copyright (c) 2002-2009 NICONET e.V.
7
C
8
C     This program is free software: you can redistribute it and/or
9
C     modify it under the terms of the GNU General Public License as
10
C     published by the Free Software Foundation, either version 2 of
11
C     the License, or (at your option) any later version.
12
C
13
C     This program is distributed in the hope that it will be useful,
14
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
15
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16
C     GNU General Public License for more details.
17
C
18
C     You should have received a copy of the GNU General Public License
19
C     along with this program.  If not, see
20
C     <http://www.gnu.org/licenses/>.
21
C
22
C     PURPOSE
23
C
24
C     To solve for X = op(U)'*op(U) either the stable non-negative
25
C     definite continuous-time Lyapunov equation
26
C                                   2
27
C        op(A)'*X + X*op(A) = -scale *op(B)'*op(B)                   (1)
28
C
29
C     or the convergent non-negative definite discrete-time Lyapunov
30
C     equation
31
C                                   2
32
C        op(A)'*X*op(A) - X = -scale *op(B)'*op(B)                   (2)
33
C
34
C     where op(K) = K or K' (i.e., the transpose of the matrix K), A is
35
C     an N-by-N matrix in real Schur form, op(B) is an M-by-N matrix,
36
C     U is an upper triangular matrix containing the Cholesky factor of
37
C     the solution matrix X, X = op(U)'*op(U), and scale is an output
38
C     scale factor, set less than or equal to 1 to avoid overflow in X.
39
C     If matrix B has full rank then the solution matrix X will be
40
C     positive-definite and hence the Cholesky factor U will be
41
C     nonsingular, but if B is rank deficient then X may only be
42
C     positive semi-definite and U will be singular.
43
C
44
C     In the case of equation (1) the matrix A must be stable (that
45
C     is, all the eigenvalues of A must have negative real parts),
46
C     and for equation (2) the matrix A must be convergent (that is,
47
C     all the eigenvalues of A must lie inside the unit circle).
48
C
49
C     ARGUMENTS
50
C
51
C     Mode Parameters
52
C
53
C     DISCR   LOGICAL
54
C             Specifies the type of Lyapunov equation to be solved as
55
C             follows:
56
C             = .TRUE. :  Equation (2), discrete-time case;
57
C             = .FALSE.:  Equation (1), continuous-time case.
58
C
59
C     LTRANS  LOGICAL
60
C             Specifies the form of op(K) to be used, as follows:
61
C             = .FALSE.:  op(K) = K    (No transpose);
62
C             = .TRUE. :  op(K) = K**T (Transpose).
63
C
64
C     Input/Output Parameters
65
C
66
C     N       (input) INTEGER
67
C             The order of the matrix A and the number of columns in
68
C             matrix op(B).  N >= 0.
69
C
70
C     M       (input) INTEGER
71
C             The number of rows in matrix op(B).  M >= 0.
72
C
73
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
74
C             The leading N-by-N upper Hessenberg part of this array
75
C             must contain a real Schur form matrix S. The elements
76
C             below the upper Hessenberg part of the array A are not
77
C             referenced. The 2-by-2 blocks must only correspond to
78
C             complex conjugate pairs of eigenvalues (not to real
79
C             eigenvalues).
80
C
81
C     LDA     INTEGER
82
C             The leading dimension of array A.  LDA >= MAX(1,N).
83
C
84
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
85
C             if LTRANS = .FALSE., and dimension (LDB,M), if
86
C             LTRANS = .TRUE..
87
C             On entry, if LTRANS = .FALSE., the leading M-by-N part of
88
C             this array must contain the coefficient matrix B of the
89
C             equation.
90
C             On entry, if LTRANS = .TRUE., the leading N-by-M part of
91
C             this array must contain the coefficient matrix B of the
92
C             equation.
93
C             On exit, if LTRANS = .FALSE., the leading
94
C             MIN(M,N)-by-MIN(M,N) upper triangular part of this array
95
C             contains the upper triangular matrix R (as defined in
96
C             METHOD), and the M-by-MIN(M,N) strictly lower triangular
97
C             part together with the elements of the array TAU are
98
C             overwritten by details of the matrix P (also defined in
99
C             METHOD). When M < N, columns (M+1),...,N of the array B
100
C             are overwritten by the matrix Z (see METHOD).
101
C             On exit, if LTRANS = .TRUE., the leading
102
C             MIN(M,N)-by-MIN(M,N) upper triangular part of
103
C             B(1:N,M-N+1), if M >= N, or of B(N-M+1:N,1:M), if M < N,
104
C             contains the upper triangular matrix R (as defined in
105
C             METHOD), and the remaining elements (below the diagonal
106
C             of R) together with the elements of the array TAU are
107
C             overwritten by details of the matrix P (also defined in
108
C             METHOD). When M < N, rows 1,...,(N-M) of the array B
109
C             are overwritten by the matrix Z (see METHOD).
110
C
111
C     LDB     INTEGER
112
C             The leading dimension of array B.
113
C             LDB >= MAX(1,M), if LTRANS = .FALSE.,
114
C             LDB >= MAX(1,N), if LTRANS = .TRUE..
115
C
116
C     TAU     (output) DOUBLE PRECISION array of dimension (MIN(N,M))
117
C             This array contains the scalar factors of the elementary
118
C             reflectors defining the matrix P.
119
C
120
C     U       (output) DOUBLE PRECISION array of dimension (LDU,N)
121
C             The leading N-by-N upper triangular part of this array
122
C             contains the Cholesky factor of the solution matrix X of
123
C             the problem, X = op(U)'*op(U).
124
C             The array U may be identified with B in the calling
125
C             statement, if B is properly dimensioned, and the
126
C             intermediate results returned in B are not needed.
127
C
128
C     LDU     INTEGER
129
C             The leading dimension of array U.  LDU >= MAX(1,N).
130
C
131
C     SCALE   (output) DOUBLE PRECISION
132
C             The scale factor, scale, set less than or equal to 1 to
133
C             prevent the solution overflowing.
134
C
135
C     Workspace
136
C
137
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
138
C             On exit, if INFO = 0, or INFO = 1, DWORK(1) returns the
139
C             optimal value of LDWORK.
140
C
141
C     LDWORK  INTEGER
142
C             The length of the array DWORK. LDWORK >= MAX(1,4*N).
143
C             For optimum performance LDWORK should sometimes be larger.
144
C
145
C     Error Indicator
146
C
147
C     INFO    INTEGER
148
C             = 0:  successful exit;
149
C             < 0:  if INFO = -i, the i-th argument had an illegal
150
C                   value;
151
C             = 1:  if the Lyapunov equation is (nearly) singular
152
C                   (warning indicator);
153
C                   if DISCR = .FALSE., this means that while the matrix
154
C                   A has computed eigenvalues with negative real parts,
155
C                   it is only just stable in the sense that small
156
C                   perturbations in A can make one or more of the
157
C                   eigenvalues have a non-negative real part;
158
C                   if DISCR = .TRUE., this means that while the matrix
159
C                   A has computed eigenvalues inside the unit circle,
160
C                   it is nevertheless only just convergent, in the
161
C                   sense that small perturbations in A can make one or
162
C                   more of the eigenvalues lie outside the unit circle;
163
C                   perturbed values were used to solve the equation
164
C                   (but the matrix A is unchanged);
165
C             = 2:  if matrix A is not stable (that is, one or more of
166
C                   the eigenvalues of A has a non-negative real part),
167
C                   if DISCR = .FALSE., or not convergent (that is, one
168
C                   or more of the eigenvalues of A lies outside the
169
C                   unit circle), if DISCR = .TRUE.;
170
C             = 3:  if matrix A has two or more consecutive non-zero
171
C                   elements on the first sub-diagonal, so that there is
172
C                   a block larger than 2-by-2 on the diagonal;
173
C             = 4:  if matrix A has a 2-by-2 diagonal block with real
174
C                   eigenvalues instead of a complex conjugate pair.
175
C
176
C     METHOD
177
C
178
C     The method used by the routine is based on the Bartels and
179
C     Stewart method [1], except that it finds the upper triangular
180
C     matrix U directly without first finding X and without the need
181
C     to form the normal matrix op(B)'*op(B) [2].
182
C
183
C     If LTRANS = .FALSE., the matrix B is factored as
184
C
185
C        B = P ( R ),  M >= N,   B = P ( R  Z ),  M < N,
186
C              ( 0 )
187
C
188
C     (QR factorization), where P is an M-by-M orthogonal matrix and
189
C     R is a square upper triangular matrix.
190
C
191
C     If LTRANS = .TRUE., the matrix B is factored as
192
C
193
C        B = ( 0  R ) P,  M >= N,  B = ( Z ) P,  M < N,
194
C                                      ( R )
195
C
196
C     (RQ factorization), where P is an M-by-M orthogonal matrix and
197
C     R is a square upper triangular matrix.
198
C
199
C     These factorizations are used to solve the continuous-time
200
C     Lyapunov equation in the canonical form
201
C                                                        2
202
C       op(A)'*op(U)'*op(U) + op(U)'*op(U)*op(A) = -scale *op(F)'*op(F),
203
C
204
C     or the discrete-time Lyapunov equation in the canonical form
205
C                                                        2
206
C       op(A)'*op(U)'*op(U)*op(A) - op(U)'*op(U) = -scale *op(F)'*op(F),
207
C
208
C     where U and F are N-by-N upper triangular matrices, and
209
C
210
C        F = R,                                  if M >= N, or
211
C
212
C        F = ( R ),    if LTRANS = .FALSE.,  or
213
C            ( 0 )
214
C
215
C        F = ( 0  Z ), if LTRANS = .TRUE.,       if M < N.
216
C            ( 0  R )
217
C
218
C     The canonical equation is solved for U.
219
C
220
C     REFERENCES
221
C
222
C     [1] Bartels, R.H. and Stewart, G.W.
223
C         Solution of the matrix equation  A'X + XB = C.
224
C         Comm. A.C.M., 15, pp. 820-826, 1972.
225
C
226
C     [2] Hammarling, S.J.
227
C         Numerical solution of the stable, non-negative definite
228
C         Lyapunov equation.
229
C         IMA J. Num. Anal., 2, pp. 303-325, 1982.
230
C
231
C     NUMERICAL ASPECTS
232
C                               3
233
C     The algorithm requires 0(N ) operations and is backward stable.
234
C
235
C     FURTHER COMMENTS
236
C
237
C     The Lyapunov equation may be very ill-conditioned. In particular,
238
C     if A is only just stable (or convergent) then the Lyapunov
239
C     equation will be ill-conditioned. "Large" elements in U relative
240
C     to those of A and B, or a "small" value for scale, are symptoms
241
C     of ill-conditioning. A condition estimate can be computed using
242
C     SLICOT Library routine SB03MD.
243
C
244
C     CONTRIBUTOR
245
C
246
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, May 1997.
247
C     Supersedes Release 2.0 routine SB03CZ by Sven Hammarling,
248
C     NAG Ltd, United Kingdom.
249
C     Partly based on routine PLYAPS by A. Varga, University of Bochum,
250
C     May 1992.
251
C
252
C     REVISIONS
253
C
254
C     Dec. 1997, April 1998, May 1999.
255
C
256
C     KEYWORDS
257
C
258
C     Lyapunov equation, orthogonal transformation, real Schur form.
259
C
260
C     ******************************************************************
261
C
262
C     .. Parameters ..
263
      DOUBLE PRECISION  ZERO, ONE
264
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
265
C     .. Scalar Arguments ..
266
      LOGICAL           DISCR, LTRANS
267
      INTEGER           INFO, LDA, LDB, LDU, LDWORK, M, N
268
      DOUBLE PRECISION  SCALE
269
C     .. Array Arguments ..
270
      DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), TAU(*), U(LDU,*)
271
C     .. Local Scalars ..
272
      INTEGER           I, J, K, L, MN, WRKOPT
273
C     .. External Subroutines ..
274
      EXTERNAL          DCOPY, DGEQRF, DGERQF, DLACPY, DLASET, SB03OT,
275
     $                  XERBLA
276
C     .. Intrinsic Functions ..
277
      INTRINSIC         MAX, MIN
278
C     .. Executable Statements ..
279
C
280
      INFO = 0
281
C
282
C     Test the input scalar arguments.
283
C
284
      IF( N.LT.0 ) THEN
285
         INFO = -3
286
      ELSE IF( M.LT.0 ) THEN
287
         INFO = -4
288
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
289
         INFO = -6
290
      ELSE IF( ( LDB.LT.MAX( 1, M ) .AND. .NOT.LTRANS ) .OR.
291
     $         ( LDB.LT.MAX( 1, N ) .AND.      LTRANS ) ) THEN
292
         INFO = -8
293
      ELSE IF( LDU.LT.MAX( 1, N ) ) THEN
294
         INFO = -11
295
      ELSE IF( LDWORK.LT.MAX( 1, 4*N ) ) THEN
296
         INFO = -14
297
      END IF
298
C
299
      IF ( INFO.NE.0 ) THEN
300
C
301
C        Error return.
302
C
303
         CALL XERBLA( 'SB03OU', -INFO )
304
         RETURN
305
      END IF
306
C
307
C     Quick return if possible.
308
C
309
      MN = MIN( N, M )
310
      IF ( MN.EQ.0 ) THEN
311
         SCALE = ONE
312
         DWORK(1) = ONE
313
         RETURN
314
      END IF
315
C
316
C     (Note: Comments in the code beginning "Workspace:" describe the
317
C     minimal amount of real workspace needed at that point in the
318
C     code, as well as the preferred amount for good performance.
319
C     NB refers to the optimal block size for the immediately
320
C     following subroutine, as returned by ILAENV.)
321
C
322
      IF ( LTRANS ) THEN
323
C
324
C        Case op(K) = K'.
325
C
326
C        Perform the RQ factorization of B.
327
C        Workspace: need   N;
328
C                   prefer N*NB.
329
C
330
         CALL DGERQF( N, M, B, LDB, TAU, DWORK, LDWORK, INFO )
331
C
332
C        The triangular matrix F is constructed in the array U so that
333
C        U can share the same memory as B.
334
C
335
         IF ( M.GE.N ) THEN
336
            CALL DLACPY( 'Upper', MN, N, B(1,M-N+1), LDB, U, LDU )
337
         ELSE
338
C
339
            DO 10 I = M, 1, -1
340
               CALL DCOPY( N-M+I, B(1,I), 1, U(1,N-M+I), 1 )
341
   10       CONTINUE
342
C
343
            CALL DLASET( 'Full', N, N-M, ZERO, ZERO, U, LDU )
344
         END IF
345
      ELSE
346
C
347
C        Case op(K) = K.
348
C
349
C        Perform the QR factorization of B.
350
C        Workspace: need   N;
351
C                   prefer N*NB.
352
C
353
         CALL DGEQRF( M, N, B, LDB, TAU, DWORK, LDWORK, INFO )
354
         CALL DLACPY( 'Upper', MN, N, B, LDB, U, LDU )
355
         IF ( M.LT.N )
356
     $      CALL DLASET( 'Upper', N-M, N-M, ZERO, ZERO, U(M+1,M+1),
357
     $                   LDU )
358
      END IF
359
      WRKOPT = DWORK(1)
360
C
361
C     Solve the canonical Lyapunov equation
362
C                                                      2
363
C     op(A)'*op(U)'*op(U) + op(U)'*op(U)*op(A) = -scale *op(F)'*op(F),
364
C
365
C     or
366
C                                                      2
367
C     op(A)'*op(U)'*op(U)*op(A) - op(U)'*op(U) = -scale *op(F)'*op(F)
368
C
369
C     for U.
370
C
371
      CALL SB03OT( DISCR, LTRANS, N, A, LDA, U, LDU, SCALE, DWORK,
372
     $             INFO )
373
      IF ( INFO.NE.0 .AND. INFO.NE.1 )
374
     $   RETURN
375
C
376
C     Make the diagonal elements of U non-negative.
377
C
378
      IF ( LTRANS ) THEN
379
C
380
         DO 30 J = 1, N
381
            IF ( U(J,J).LT.ZERO ) THEN
382
C
383
               DO 20 I = 1, J
384
                  U(I,J) = -U(I,J)
385
   20          CONTINUE
386
C
387
            END IF
388
   30    CONTINUE
389
C
390
      ELSE
391
         K = 1
392
C
393
         DO 50 J = 1, N
394
            DWORK(K) = U(J,J)
395
            L = 1
396
C
397
            DO 40 I = 1, J
398
               IF ( DWORK(L).LT.ZERO ) U(I,J) = -U(I,J)
399
               L = L + 1
400
   40       CONTINUE
401
C
402
            K = K + 1
403
   50    CONTINUE
404
C
405
      END IF
406
C
407
      DWORK(1) = MAX( WRKOPT, 4*N )
408
      RETURN
409
C *** Last line of SB03OU ***
410
      END

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks