← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/slicot/sb03ou.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

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