← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/slicot/sb10dd.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/sb10dd.f
 
1
      SUBROUTINE SB10DD( N, M, NP, NCON, NMEAS, GAMMA, A, LDA, B, LDB, 
 
2
     $                   C, LDC, D, LDD, AK, LDAK, BK, LDBK, CK, LDCK,
 
3
     $                   DK, LDDK, X, LDX, Z, LDZ, RCOND, TOL, IWORK,
 
4
     $                   DWORK, LDWORK, BWORK, INFO )
 
5
C
 
6
C     RELEASE 4.0, WGS COPYRIGHT 1999.
 
7
C
 
8
C     PURPOSE
 
9
C
 
10
C     To compute the matrices of an H-infinity (sub)optimal n-state
 
11
C     controller
 
12
C
 
13
C                           | AK | BK |
 
14
C                       K = |----|----|,
 
15
C                           | CK | DK |
 
16
C
 
17
C     for the discrete-time system
 
18
C
 
19
C                   | A  | B1  B2  |   | A | B |
 
20
C               P = |----|---------| = |---|---|
 
21
C                   | C1 | D11 D12 |   | C | D | 
 
22
C                   | C2 | D21 D22 |
 
23
C
 
24
C     and for a given value of gamma, where B2 has as column size the
 
25
C     number of control inputs (NCON) and C2 has as row size the number
 
26
C     of measurements (NMEAS) being provided to the controller.
 
27
C
 
28
C     It is assumed that 
 
29
C
 
30
C     (A1) (A,B2) is stabilizable and (C2,A) is detectable,
 
31
C
 
32
C     (A2) D12 is full column rank and D21 is full row rank,
 
33
C
 
34
C               j*Theta 
 
35
C     (A3) | A-e       *I  B2  | has full column rank for all 
 
36
C          |    C1         D12 |
 
37
C
 
38
C          0 <= Theta < 2*Pi ,
 
39
C
 
40
C               j*Theta
 
41
C     (A4) | A-e       *I  B1  |  has full row rank for all
 
42
C          |    C2         D21 |  
 
43
C          
 
44
C          0 <= Theta < 2*Pi .
 
45
C
 
46
C     ARGUMENTS
 
47
C
 
48
C     Input/Output Parameters
 
49
C
 
50
C     N       (input) INTEGER
 
51
C             The order of the system.  N >= 0.
 
52
C
 
53
C     M       (input) INTEGER
 
54
C             The column size of the matrix B.  M >= 0.
 
55
C
 
56
C     NP      (input) INTEGER
 
57
C             The row size of the matrix C.  NP >= 0.
 
58
C
 
59
C     NCON    (input) INTEGER
 
60
C             The number of control inputs (M2).  M >= NCON >= 0,
 
61
C             NP-NMEAS >= NCON.
 
62
C
 
63
C     NMEAS   (input) INTEGER
 
64
C             The number of measurements (NP2).  NP >= NMEAS >= 0,
 
65
C             M-NCON >= NMEAS.
 
66
C
 
67
C     GAMMA   (input) DOUBLE PRECISION
 
68
C             The value of gamma. It is assumed that gamma is
 
69
C             sufficiently large so that the controller is admissible.
 
70
C             GAMMA > 0.
 
71
C
 
72
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
 
73
C             The leading N-by-N part of this array must contain the
 
74
C             system state matrix A.
 
75
C
 
76
C     LDA     INTEGER
 
77
C             The leading dimension of the array A.  LDA >= max(1,N).
 
78
C
 
79
C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
 
80
C             The leading N-by-M part of this array must contain the
 
81
C             system input matrix B.
 
82
C
 
83
C     LDB     INTEGER
 
84
C             The leading dimension of the array B.  LDB >= max(1,N).
 
85
C
 
86
C     C       (input) DOUBLE PRECISION array, dimension (LDC,N)
 
87
C             The leading NP-by-N part of this array must contain the
 
88
C             system output matrix C.
 
89
C
 
90
C     LDC     INTEGER
 
91
C             The leading dimension of the array C.  LDC >= max(1,NP).
 
92
C
 
93
C     D       (input) DOUBLE PRECISION array, dimension (LDD,M)
 
94
C             The leading NP-by-M part of this array must contain the
 
95
C             system input/output matrix D.
 
96
C           
 
97
C     LDD     INTEGER
 
98
C             The leading dimension of the array D.  LDD >= max(1,NP).
 
99
C
 
100
C     AK      (output) DOUBLE PRECISION array, dimension (LDAK,N)
 
101
C             The leading N-by-N part of this array contains the
 
102
C             controller state matrix AK.
 
103
C
 
104
C     LDAK    INTEGER
 
105
C             The leading dimension of the array AK.  LDAK >= max(1,N).
 
106
C
 
107
C     BK      (output) DOUBLE PRECISION array, dimension (LDBK,NMEAS)
 
108
C             The leading N-by-NMEAS part of this array contains the
 
109
C             controller input matrix BK.
 
110
C
 
111
C     LDBK    INTEGER
 
112
C             The leading dimension of the array BK.  LDBK >= max(1,N).
 
113
C
 
114
C     CK      (output) DOUBLE PRECISION array, dimension (LDCK,N)
 
115
C             The leading NCON-by-N part of this array contains the
 
116
C             controller output matrix CK.
 
117
C
 
118
C     LDCK    INTEGER
 
119
C             The leading dimension of the array CK.
 
120
C             LDCK >= max(1,NCON).
 
121
C
 
122
C     DK      (output) DOUBLE PRECISION array, dimension (LDDK,NMEAS)
 
123
C             The leading NCON-by-NMEAS part of this array contains the
 
124
C             controller input/output matrix DK.
 
125
C
 
126
C     LDDK    INTEGER
 
127
C             The leading dimension of the array DK.
 
128
C             LDDK >= max(1,NCON).
 
129
C
 
130
C     X       (output) DOUBLE PRECISION array, dimension (LDX,N)
 
131
C             The leading N-by-N part of this array contains the matrix
 
132
C             X, solution of the X-Riccati equation.
 
133
C
 
134
C     LDX     INTEGER
 
135
C             The leading dimension of the array X.  LDX >= max(1,N).
 
136
C
 
137
C     Z       (output) DOUBLE PRECISION array, dimension (LDZ,N)
 
138
C             The leading N-by-N part of this array contains the matrix
 
139
C             Z, solution of the Z-Riccati equation.
 
140
C
 
141
C     LDZ     INTEGER
 
142
C             The leading dimension of the array Z.  LDZ >= max(1,N).
 
143
C  
 
144
C     RCOND   (output) DOUBLE PRECISION array, dimension (8)
 
145
C             RCOND contains estimates of the reciprocal condition
 
146
C             numbers of the matrices which are to be inverted and 
 
147
C             estimates of the reciprocal condition numbers of the
 
148
C             Riccati equations which have to be solved during the
 
149
C             computation of the controller. (See the description of
 
150
C             the algorithm in [2].)
 
151
C             RCOND(1) contains the reciprocal condition number of the 
 
152
C                      matrix R3;
 
153
C             RCOND(2) contains the reciprocal condition number of the 
 
154
C                      matrix R1 - R2'*inv(R3)*R2;
 
155
C             RCOND(3) contains the reciprocal condition number of the 
 
156
C                      matrix V21;
 
157
C             RCOND(4) contains the reciprocal condition number of the 
 
158
C                      matrix St3;
 
159
C             RCOND(5) contains the reciprocal condition number of the 
 
160
C                      matrix V12;
 
161
C             RCOND(6) contains the reciprocal condition number of the 
 
162
C                      matrix Im2 + DKHAT*D22
 
163
C             RCOND(7) contains the reciprocal condition number of the 
 
164
C                      X-Riccati equation;
 
165
C             RCOND(8) contains the reciprocal condition number of the 
 
166
C                      Z-Riccati equation.
 
167
C
 
168
C     Tolerances
 
169
C
 
170
C     TOL     DOUBLE PRECISION
 
171
C             Tolerance used in neglecting the small singular values
 
172
C             in rank determination. If TOL <= 0, then a default value
 
173
C             equal to 1000*EPS is used, where EPS is the relative
 
174
C             machine precision.
 
175
C
 
176
C     Workspace
 
177
C
 
178
C     IWORK   INTEGER array, dimension max(2*max(M2,N),M,M2+NP2,N*N)
 
179
C
 
180
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
 
181
C             On exit, if INFO = 0, DWORK(1) contains the optimal
 
182
C             LDWORK.
 
183
C
 
184
C     LDWORK  INTEGER
 
185
C             The dimension of the array DWORK.
 
186
C             LDWORK >= max(LW1,LW2,LW3,LW4), where
 
187
C             LW1 = (N+NP1+1)*(N+M2) + max(3*(N+M2)+N+NP1,5*(N+M2));
 
188
C             LW2 = (N+NP2)*(N+M1+1) + max(3*(N+NP2)+N+M1,5*(N+NP2));
 
189
C             LW3 = 13*N*N + 2*M*M + N*(8*M+NP2) + M1*(M2+NP2) + 6*N + 
 
190
C                   max(14*N+23,16*N,2*N+M,3*M);
 
191
C             LW4 = 13*N*N + M*M + (8*N+M+M2+2*NP2)*(M2+NP2) + 6*N +
 
192
C                   N*(M+NP2) + max(14*N+23,16*N,2*N+M2+NP2,3*(M2+NP2));
 
193
C             For good performance, LDWORK must generally be larger.
 
194
C             Denoting Q = max(M1,M2,NP1,NP2), an upper bound is
 
195
C             max((N+Q)*(N+Q+6),13*N*N + M*M + 2*Q*Q + N*(M+Q) +
 
196
C                 max(M*(M+7*N),2*Q*(8*N+M+2*Q)) + 6*N + 
 
197
C                 max(14*N+23,16*N,2*N+max(M,2*Q),3*max(M,2*Q)).
 
198
C
 
199
C     BWORK   LOGICAL array, dimension (2*N)
 
200
C
 
201
C     Error Indicator
 
202
C
 
203
C     INFO    INTEGER
 
204
C             = 0:  successful exit;
 
205
C             < 0:  if INFO = -i, the i-th argument had an illegal
 
206
C                   value;
 
207
C                                      j*Theta
 
208
C             = 1:  if the matrix | A-e       *I  B2  | had not full
 
209
C                                 |      C1       D12 |
 
210
C                   column rank;
 
211
C                                      j*Theta
 
212
C             = 2:  if the matrix | A-e       *I  B1  | had not full
 
213
C                                 |      C2       D21 |
 
214
C                   row rank;
 
215
C             = 3:  if the matrix D12 had not full column rank;
 
216
C             = 4:  if the matrix D21 had not full row rank;
 
217
C             = 5:  if the controller is not admissible (too small value
 
218
C                   of gamma);
 
219
C             = 6:  if the X-Riccati equation was not solved
 
220
C                   successfully (the controller is not admissible or
 
221
C                   there are numerical difficulties);
 
222
C             = 7:  if the Z-Riccati equation was not solved
 
223
C                   successfully (the controller is not admissible or
 
224
C                   there are numerical difficulties);
 
225
C             = 8:  if the matrix Im2 + DKHAT*D22 is singular.
 
226
C             = 9:  if the singular value decomposition (SVD) algorithm
 
227
C                   did not converge (when computing the SVD of one of
 
228
C                   the matrices |A   B2 |, |A   B1 |, D12 or D21).
 
229
C                                |C1  D12|  |C2  D21|
 
230
C
 
231
C     METHOD
 
232
C
 
233
C     The routine implements the method presented in [1].
 
234
C
 
235
C     REFERENCES 
 
236
C
 
237
C     [1] Green, M. and Limebeer, D.J.N.
 
238
C         Linear Robust Control.
 
239
C         Prentice-Hall, Englewood Cliffs, NJ, 1995.
 
240
C
 
241
C     [2] Petkov, P.Hr., Gu, D.W., and Konstantinov, M.M.
 
242
C         Fortran 77 routines for Hinf and H2 design of linear 
 
243
C         discrete-time control systems.
 
244
C         Report 99-8, Department of Engineering, Leicester University,
 
245
C         April 1999.
 
246
C
 
247
C     NUMERICAL ASPECTS
 
248
C
 
249
C     With approaching the minimum value of gamma some of the matrices
 
250
C     which are to be inverted tend to become ill-conditioned and 
 
251
C     the X- or Z-Riccati equation may also become ill-conditioned
 
252
C     which may deteriorate the accuracy of the result. (The
 
253
C     corresponding reciprocal condition numbers are given in 
 
254
C     the output array RCOND.) 
 
255
C
 
256
C     CONTRIBUTORS
 
257
C
 
258
C     P.Hr. Petkov, D.W. Gu and M.M. Konstantinov, April 1999.
 
259
C
 
260
C     REVISIONS
 
261
C
 
262
C     V. Sima, Research Institute for Informatics, Bucharest, Sep. 1999.
 
263
C     V. Sima, Research Institute for Informatics, Bucharest, Feb. 2000.
 
264
C
 
265
C     KEYWORDS
 
266
C
 
267
C     Algebraic Riccati equation, discrete-time H-infinity optimal
 
268
C     control, robust control.
 
269
C
 
270
C     ******************************************************************
 
271
C
 
272
C     .. Parameters ..
 
273
      DOUBLE PRECISION   ZERO, ONE, THOUSN
 
274
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
 
275
     $                     THOUSN = 1.0D+3 )
 
276
C     ..
 
277
C     .. Scalar Arguments ..
 
278
      INTEGER            INFO, LDA, LDAK, LDB, LDBK, LDC, LDCK, LDD,
 
279
     $                   LDDK, LDWORK, LDX, LDZ, M, N, NCON, NMEAS, NP
 
280
      DOUBLE PRECISION   GAMMA, TOL
 
281
C     ..
 
282
C     .. Array Arguments ..
 
283
      INTEGER            IWORK( * )
 
284
      DOUBLE PRECISION   A( LDA, * ), AK( LDAK, * ), B( LDB, * ),
 
285
     $                   BK( LDBK, * ), C( LDC, * ), CK( LDCK, * ),
 
286
     $                   D( LDD, * ), DK( LDDK, * ), DWORK( * ),
 
287
     $                   RCOND( * ), X( LDX, * ), Z( LDZ, * )
 
288
      LOGICAL            BWORK( * )
 
289
C     ..
 
290
C     .. Local Scalars ..
 
291
      INTEGER            INFO2, IR2, IR3, IS2, IS3, IWB, IWC, IWD, IWG,
 
292
     $                   IWH, IWI, IWL, IWQ, IWR, IWRK, IWS, IWT, IWU,
 
293
     $                   IWV, IWW, J, LWAMAX, M1, M2, MINWRK, NP1, NP2
 
294
      DOUBLE PRECISION   ANORM, FERR, RCOND2, SEPD, TOLL
 
295
C
 
296
C     .. External Functions
 
297
      DOUBLE PRECISION   DLAMCH, DLANGE, DLANSY
 
298
      EXTERNAL           DLAMCH, DLANGE, DLANSY
 
299
C     ..
 
300
C     .. External Subroutines ..
 
301
      EXTERNAL           DGECON, DGEMM, DGESVD, DGETRF, DGETRS, DLACPY,
 
302
     $                   DLASET, DPOCON, DPOTRF, DSCAL, DSWAP, DSYRK,
 
303
     $                   DSYTRF, DSYTRS, DTRCON, DTRSM, MA02AD, MB01RU,
 
304
     $                   MB01RX, SB02OD, SB02SD, XERBLA
 
305
C     ..
 
306
C     .. Intrinsic Functions ..
 
307
      INTRINSIC          DBLE, INT, MAX
 
308
C     ..
 
309
C     .. Executable Statements ..      
 
310
C
 
311
C     Decode and Test input parameters.
 
312
C
 
313
      M1 = M - NCON
 
314
      M2 = NCON
 
315
      NP1 = NP - NMEAS
 
316
      NP2 = NMEAS
 
317
C
 
318
      INFO = 0
 
319
      IF( N.LT.0 ) THEN
 
320
         INFO = -1
 
321
      ELSE IF( M.LT.0 ) THEN
 
322
         INFO = -2
 
323
      ELSE IF( NP.LT.0 ) THEN
 
324
         INFO = -3
 
325
      ELSE IF( NCON.LT.0 .OR. M1.LT.0 .OR. M2.GT.NP1 ) THEN
 
326
         INFO = -4
 
327
      ELSE IF( NMEAS.LT.0 .OR. NP1.LT.0 .OR. NP2.GT.M1 ) THEN
 
328
         INFO = -5
 
329
      ELSE IF( GAMMA.LE.ZERO ) THEN
 
330
         INFO = -6
 
331
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 
332
         INFO = -8
 
333
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
 
334
         INFO = -10  
 
335
      ELSE IF( LDC.LT.MAX( 1, NP ) ) THEN
 
336
         INFO = -12
 
337
      ELSE IF( LDD.LT.MAX( 1, NP ) ) THEN
 
338
         INFO = -14
 
339
      ELSE IF( LDAK.LT.MAX( 1, N ) ) THEN
 
340
         INFO = -16
 
341
      ELSE IF( LDBK.LT.MAX( 1, N ) ) THEN
 
342
         INFO = -18
 
343
      ELSE IF( LDCK.LT.MAX( 1, M2 ) ) THEN
 
344
         INFO = -20
 
345
      ELSE IF( LDDK.LT.MAX( 1, M2 ) ) THEN
 
346
         INFO = -22
 
347
      ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
 
348
         INFO = -24
 
349
      ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
 
350
         INFO = -26 
 
351
      ELSE
 
352
C
 
353
C        Compute workspace.
 
354
C
 
355
         IWB = ( N + NP1 + 1 )*( N + M2 ) + 
 
356
     $         MAX( 3*( N + M2 ) + N + NP1, 5*( N + M2 ) )
 
357
         IWC = ( N + NP2 )*( N + M1 + 1 ) +  
 
358
     $         MAX( 3*( N + NP2 ) + N + M1, 5*( N + NP2 ) )
 
359
         IWD = 13*N*N + 2*M*M + N*( 8*M + NP2 ) + M1*( M2 + NP2 ) + 
 
360
     $         6*N + MAX( 14*N + 23, 16*N, 2*N + M, 3*M )
 
361
         IWG = 13*N*N + M*M + ( 8*N + M + M2 + 2*NP2 )*( M2 + NP2 ) +
 
362
     $         6*N + N*( M + NP2 ) + 
 
363
     $         MAX( 14*N + 23, 16*N, 2*N + M2 + NP2, 3*( M2 + NP2 ) )
 
364
         MINWRK = MAX( IWB, IWC, IWD, IWG )
 
365
         IF( LDWORK.LT.MINWRK )
 
366
     $      INFO = -31
 
367
      END IF       
 
368
      IF( INFO.NE.0 ) THEN
 
369
         CALL XERBLA( 'SB10DD', -INFO )
 
370
         RETURN
 
371
      END IF
 
372
C
 
373
C     Quick return if possible.
 
374
C
 
375
      IF( N.EQ.0 .OR. M.EQ.0 .OR. NP.EQ.0 .OR. M1.EQ.0 .OR. M2.EQ.0
 
376
     $    .OR. NP1.EQ.0 .OR. NP2.EQ.0 ) THEN
 
377
          RCOND( 1 ) = ONE
 
378
          RCOND( 2 ) = ONE
 
379
          RCOND( 3 ) = ONE
 
380
          RCOND( 4 ) = ONE
 
381
          RCOND( 5 ) = ONE
 
382
          RCOND( 6 ) = ONE
 
383
          RCOND( 7 ) = ONE
 
384
          RCOND( 8 ) = ONE
 
385
          DWORK( 1 ) = ONE
 
386
          RETURN
 
387
      END IF
 
388
C
 
389
      TOLL = TOL
 
390
      IF( TOLL.LE.ZERO ) THEN
 
391
C
 
392
C        Set the default value of the tolerance in rank determination.
 
393
C
 
394
         TOLL = THOUSN*DLAMCH( 'Epsilon' )
 
395
      END IF
 
396
C
 
397
C     Workspace usage.
 
398
C
 
399
      IWS  = (N+NP1)*(N+M2) + 1
 
400
      IWRK = IWS + (N+M2)
 
401
C
 
402
C                      jTheta
 
403
C     Determine if |A-e      I  B2 | has full column rank at
 
404
C                  |     C1     D12|
 
405
C     Theta = Pi/2 .
 
406
C     Workspace: need (N+NP1+1)*(N+M2) + MAX(3*(N+M2)+N+NP1,5*(N+M2));
 
407
C                prefer larger.
 
408
C
 
409
      CALL DLACPY( 'Full', N, N, A, LDA, DWORK, N+NP1 )
 
410
      CALL DLACPY( 'Full', NP1, N, C, LDC, DWORK( N+1 ), N+NP1 )
 
411
      CALL DLACPY( 'Full', N, M2, B( 1, M1+1 ), LDB,
 
412
     $             DWORK( (N+NP1)*N+1 ), N+NP1 )
 
413
      CALL DLACPY( 'Full', NP1, M2, D( 1, M1+1 ), LDD,
 
414
     $             DWORK( (N+NP1)*N+N+1 ), N+NP1 )
 
415
      CALL DGESVD( 'N', 'N', N+NP1, N+M2, DWORK, N+NP1, DWORK( IWS ),
 
416
     $             DWORK, N+NP1, DWORK, N+M2, DWORK( IWRK ),
 
417
     $             LDWORK-IWRK+1, INFO2 )
 
418
      IF( INFO2.GT.0 ) THEN
 
419
         INFO = 9
 
420
         RETURN
 
421
      END IF
 
422
      IF( DWORK( IWS+N+M2 ) / DWORK( IWS ).LE.TOLL ) THEN
 
423
         INFO = 1
 
424
         RETURN
 
425
      END IF
 
426
      LWAMAX = INT( DWORK( IWRK ) ) + IWRK - 1
 
427
C
 
428
C     Workspace usage.
 
429
C
 
430
      IWS  = (N+NP2)*(N+M1) + 1
 
431
      IWRK = IWS + (N+NP2)
 
432
C
 
433
C                      jTheta
 
434
C     Determine if |A-e      I   B1 | has full row rank at
 
435
C                  |     C2      D21|
 
436
C     Theta = Pi/2 .
 
437
C     Workspace: need   (N+NP2)*(N+M1+1) + 
 
438
C                       MAX(3*(N+NP2)+N+M1,5*(N+NP2));
 
439
C                prefer larger.
 
440
C
 
441
      CALL DLACPY( 'Full', N, N, A, LDA, DWORK, N+NP2 )
 
442
      CALL DLACPY( 'Full', NP2, N, C( NP1+1, 1), LDC, DWORK( N+1 ),
 
443
     $             N+NP2 )     
 
444
      CALL DLACPY( 'Full', N, M1, B, LDB, DWORK( (N+NP2)*N+1 ),
 
445
     $             N+NP2 )      
 
446
      CALL DLACPY( 'Full', NP2, M1, D( NP1+1, 1 ), LDD,
 
447
     $             DWORK( (N+NP2)*N+N+1 ), N+NP2 )
 
448
      CALL DGESVD( 'N', 'N', N+NP2, N+M1, DWORK, N+NP2, DWORK( IWS ),
 
449
     $             DWORK, N+NP2, DWORK, N+M1, DWORK( IWRK ),
 
450
     $             LDWORK-IWRK+1, INFO2 )
 
451
      IF( INFO2.GT.0 ) THEN
 
452
         INFO = 9
 
453
         RETURN
 
454
      END IF
 
455
      IF( DWORK( IWS+N+NP2 ) / DWORK( IWS ).LE.TOLL ) THEN
 
456
         INFO = 2
 
457
         RETURN
 
458
      END IF
 
459
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
460
C
 
461
C     Workspace usage.
 
462
C
 
463
      IWS  = NP1*M2 + 1
 
464
      IWRK = IWS + M2
 
465
C
 
466
C     Determine if D12 has full column rank.
 
467
C     Workspace: need   (NP1+1)*M2 + MAX(3*M2+NP1,5*M2);
 
468
C                prefer larger.
 
469
C
 
470
      CALL DLACPY( 'Full', NP1, M2, D( 1, M1+1 ), LDD, DWORK, NP1 )
 
471
      CALL DGESVD( 'N', 'N', NP1, M2, DWORK, NP1, DWORK( IWS ), DWORK,
 
472
     $             NP1, DWORK, M2, DWORK( IWRK ), LDWORK-IWRK+1, INFO2 )
 
473
      IF( INFO2.GT.0 ) THEN
 
474
         INFO = 9
 
475
         RETURN
 
476
      END IF
 
477
      IF( DWORK( IWS+M2 ) / DWORK( IWS ).LE.TOLL ) THEN
 
478
         INFO = 3
 
479
         RETURN
 
480
      END IF
 
481
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
482
C
 
483
C     Workspace usage.
 
484
C
 
485
      IWS  = NP2*M1 + 1
 
486
      IWRK = IWS + NP2
 
487
C
 
488
C     Determine if D21 has full row rank.
 
489
C     Workspace: need   NP2*(M1+1) + MAX(3*NP2+M1,5*NP2);
 
490
C                prefer larger.
 
491
C
 
492
      CALL DLACPY( 'Full', NP2, M1, D( NP1+1, 1 ), LDD, DWORK, NP2 )
 
493
      CALL DGESVD( 'N', 'N', NP2, M1, DWORK, NP2, DWORK( IWS ), DWORK,
 
494
     $             NP2, DWORK, M1, DWORK( IWRK ), LDWORK-IWRK+1, INFO2 )
 
495
      IF( INFO2.GT.0 ) THEN
 
496
         INFO = 9
 
497
         RETURN
 
498
      END IF
 
499
      IF( DWORK( IWS+NP2 ) / DWORK( IWS ).LE.TOLL ) THEN
 
500
         INFO = 4
 
501
         RETURN
 
502
      END IF
 
503
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
504
C
 
505
C     Workspace usage.
 
506
C
 
507
      IWV  = 1
 
508
      IWB  = IWV + M*M
 
509
      IWC  = IWB + N*M1
 
510
      IWD  = IWC + ( M2 + NP2 )*N
 
511
      IWQ  = IWD + ( M2 + NP2 )*M1 
 
512
      IWL  = IWQ + N*N
 
513
      IWR  = IWL + N*M
 
514
      IWI  = IWR + 2*N
 
515
      IWH  = IWI + 2*N
 
516
      IWS  = IWH + 2*N
 
517
      IWT  = IWS + ( 2*N + M )*( 2*N + M )
 
518
      IWU  = IWT + ( 2*N + M )*2*N 
 
519
      IWRK = IWU + 4*N*N
 
520
      IR2  = IWV + M1
 
521
      IR3  = IR2 + M*M1
 
522
C
 
523
C     Compute R0 = |D11'||D11 D12| -|gamma^2*Im1 0| .
 
524
C                  |D12'|           |   0        0|
 
525
C
 
526
      CALL DSYRK( 'Lower', 'Transpose', M, NP1, ONE, D, LDD, ZERO,
 
527
     $            DWORK, M )
 
528
      DO 10 J = 1, M*M1, M + 1
 
529
         DWORK( J ) = DWORK( J ) - GAMMA*GAMMA
 
530
   10 CONTINUE
 
531
C
 
532
C     Compute C1'*C1 .
 
533
C
 
534
      CALL DSYRK( 'Lower', 'Transpose', N, NP1, ONE, C, LDC, ZERO, 
 
535
     $            DWORK( IWQ ), N )
 
536
C
 
537
C     Compute C1'*|D11 D12| .
 
538
C
 
539
      CALL DGEMM( 'Transpose', 'NoTranspose', N, M, NP1, ONE, C, LDC, 
 
540
     $            D, LDD, ZERO, DWORK( IWL ), N )
 
541
C
 
542
C     Solution of the X-Riccati equation.
 
543
C     Workspace: need   13*N*N + 2*M*M + N*(8*M+NP2) + M1*(M2+NP2) + 
 
544
C                       6*N + max(14*N+23,16*N,2*N+M,3*M);
 
545
C                prefer larger.
 
546
C
 
547
      CALL SB02OD( 'D', 'B', 'N', 'L', 'N', 'S', N, M, NP, A, LDA, B,
 
548
     $             LDB, DWORK( IWQ ), N, DWORK, M, DWORK( IWL ), N,
 
549
     $             RCOND2, X, LDX, DWORK( IWR ), DWORK( IWI ),
 
550
     $             DWORK( IWH ), DWORK( IWS ), 2*N+M, DWORK( IWT ),
 
551
     $             2*N+M, DWORK( IWU ), 2*N, TOLL, IWORK,
 
552
     $             DWORK( IWRK ), LDWORK-IWRK+1, BWORK, INFO2 )
 
553
      IF( INFO2.GT.0 ) THEN
 
554
         INFO = 6
 
555
         RETURN
 
556
      END IF
 
557
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
558
C
 
559
C     Condition estimation.
 
560
C     Workspace: need   4*N*N + 2*M*M + N*(3*M+NP2) + M1*(M2+NP2) +
 
561
C                       max(5*N,max(3,2*N*N)+N*N);
 
562
C                prefer larger.
 
563
C
 
564
      IWS  = IWR
 
565
      IWH  = IWS + M*M
 
566
      IWT  = IWH + N*M
 
567
      IWU  = IWT + N*N 
 
568
      IWG  = IWU + N*N 
 
569
      IWRK = IWG + N*N
 
570
      CALL DLACPY( 'Lower', M, M, DWORK, M, DWORK( IWS ), M )
 
571
      CALL DSYTRF( 'Lower', M, DWORK( IWS ), M, IWORK, DWORK( IWRK ),
 
572
     $             LDWORK-IWRK+1, INFO2 )
 
573
      IF( INFO2.GT.0 ) THEN
 
574
         INFO = 5
 
575
         RETURN
 
576
      END IF
 
577
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
578
C
 
579
      CALL MA02AD( 'Full', N, M, B, LDB, DWORK( IWH ), M )
 
580
      CALL DSYTRS( 'Lower', M, N, DWORK( IWS ), M, IWORK, DWORK( IWH ), 
 
581
     $             M, INFO2 )
 
582
      CALL MB01RX( 'Left', 'Lower', 'NoTranspose', N, M, ZERO, ONE,
 
583
     $             DWORK( IWG ), N, B, LDB, DWORK( IWH ), M, INFO2 )
 
584
      CALL SB02SD( 'C', 'N', 'N', 'L', 'O', N, A, LDA, DWORK( IWT ), N,
 
585
     $             DWORK( IWU ), N, DWORK( IWG ), N, DWORK( IWQ ), N, X,
 
586
     $             LDX, SEPD, RCOND( 7 ), FERR, IWORK, DWORK( IWRK ),
 
587
     $             LDWORK-IWRK+1, INFO2 )
 
588
      IF( INFO2.GT.0 ) RCOND( 7 ) = ZERO
 
589
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
590
C
 
591
C     Workspace usage.
 
592
C
 
593
      IWRK = IWR
 
594
C
 
595
C     Compute the lower triangle of |R1  R2'| = R0 + B'*X*B .
 
596
C                                   |R2  R3 |
 
597
C
 
598
      CALL MB01RU( 'Lower', 'Transpose', M, N, ONE, ONE, DWORK, M,
 
599
     $             B, LDB, X, LDX, DWORK( IWRK ), M*N, INFO2 )
 
600
C
 
601
C     Compute the Cholesky factorization of R3, R3 = V12'*V12 .
 
602
C     Note that V12' is stored.
 
603
C
 
604
      ANORM = DLANSY( '1', 'Lower', M2, DWORK( IR3 ), M, DWORK( IWRK ) )
 
605
      CALL DPOTRF( 'Lower', M2, DWORK( IR3 ), M, INFO2 )
 
606
      IF( INFO2.GT.0 ) THEN
 
607
         INFO = 5
 
608
         RETURN
 
609
      END IF
 
610
      CALL DPOCON( 'Lower', M2, DWORK( IR3 ), M, ANORM, RCOND( 1 ),
 
611
     $             DWORK( IWRK ), IWORK, INFO2 )
 
612
C
 
613
C     Return if the matrix is singular to working precision.
 
614
C
 
615
      IF( RCOND( 1 ).LT.TOLL ) THEN
 
616
         INFO = 5
 
617
         RETURN
 
618
      END IF     
 
619
C
 
620
      CALL DTRCON( '1', 'Lower', 'NonUnit', M2, DWORK( IR3 ), M,
 
621
     $             RCOND( 5 ), DWORK( IWRK ), IWORK, INFO2 )
 
622
C
 
623
C     Return if the matrix is singular to working precision.
 
624
C
 
625
      IF( RCOND( 5 ).LT.TOLL ) THEN
 
626
         INFO = 5
 
627
         RETURN
 
628
      END IF     
 
629
C
 
630
C     Compute R2 <- inv(V12')*R2 .
 
631
C
 
632
      CALL DTRSM( 'Left', 'Lower', 'NoTranspose', 'NonUnit', M2, M1,
 
633
     $            ONE, DWORK( IR3 ), M, DWORK( IR2 ), M )
 
634
C                   
 
635
C     Compute -Nabla = R2'*inv(R3)*R2 - R1 .
 
636
C
 
637
      CALL DSYRK( 'Lower', 'Transpose', M1, M2, ONE, DWORK( IR2 ), M,
 
638
     $            -ONE, DWORK, M )
 
639
C
 
640
C     Compute the Cholesky factorization of -Nabla, -Nabla = V21t'*V21t.
 
641
C     Note that V21t' is stored.
 
642
C
 
643
      ANORM = DLANSY( '1', 'Lower', M1, DWORK, M, DWORK( IWRK ) )
 
644
      CALL DPOTRF( 'Lower', M1, DWORK, M, INFO2 )
 
645
      IF( INFO2.GT.0 ) THEN
 
646
         INFO = 5
 
647
         RETURN
 
648
      END IF
 
649
      CALL DPOCON( 'Lower', M1, DWORK, M, ANORM, RCOND( 2 ),
 
650
     $             DWORK( IWRK ), IWORK, INFO2 )
 
651
C
 
652
C     Return if the matrix is singular to working precision.
 
653
C
 
654
      IF( RCOND( 2 ).LT.TOLL ) THEN
 
655
         INFO = 5
 
656
         RETURN
 
657
      END IF     
 
658
C
 
659
      CALL DTRCON( '1', 'Lower', 'NonUnit', M1, DWORK, M, RCOND( 3 ),
 
660
     $             DWORK( IWRK ), IWORK, INFO2 )
 
661
C
 
662
C     Return if the matrix is singular to working precision.
 
663
C
 
664
      IF( RCOND( 3 ).LT.TOLL ) THEN
 
665
         INFO = 5
 
666
         RETURN
 
667
      END IF     
 
668
C
 
669
C     Compute X*A .
 
670
C
 
671
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, N, ONE, X, LDX,
 
672
     $            A, LDA, ZERO, DWORK( IWQ ), N )
 
673
C
 
674
C     Compute |L1| = |D11'|*C1 + B'*X*A .
 
675
C             |L2| = |D12'| 
 
676
C 
 
677
      CALL MA02AD( 'Full', N, M, DWORK( IWL ), N, DWORK( IWRK ), M )
 
678
      CALL DLACPY( 'Full', M, N, DWORK( IWRK ), M, DWORK( IWL ), M )
 
679
      CALL DGEMM( 'Transpose', 'NoTranspose', M, N, N, ONE, B, LDB, 
 
680
     $            DWORK( IWQ ), N, ONE, DWORK( IWL ), M )
 
681
C
 
682
C     Compute L2 <- inv(V12')*L2 .
 
683
C
 
684
      CALL DTRSM( 'Left', 'Lower', 'NoTranspose', 'NonUnit', M2, N, ONE,
 
685
     $            DWORK( IR3 ), M, DWORK( IWL+M1 ), M )
 
686
C
 
687
C     Compute L_Nabla = L1 - R2'*inv(R3)*L2 .
 
688
C
 
689
      CALL DGEMM( 'Transpose', 'NoTranspose', M1, N, M2, -ONE,
 
690
     $            DWORK( IR2 ), M, DWORK( IWL+M1 ), M, ONE, 
 
691
     $            DWORK( IWL ), M )
 
692
C
 
693
C     Compute L_Nabla <- inv(V21t')*L_Nabla .
 
694
C
 
695
      CALL DTRSM( 'Left', 'Lower', 'NoTranspose', 'NonUnit', M1, N, ONE,
 
696
     $            DWORK, M, DWORK( IWL ), M )
 
697
C
 
698
C     Compute Bt1 = B1*inv(V21t) .
 
699
C
 
700
      CALL DLACPY( 'Full', N, M1, B, LDB, DWORK( IWB ), N )
 
701
      CALL DTRSM( 'Right', 'Lower', 'Transpose', 'NonUnit', N, M1, ONE,
 
702
     $            DWORK, M, DWORK( IWB ), N )
 
703
C
 
704
C     Compute At .
 
705
C
 
706
      CALL DLACPY( 'Full', N, N, A, LDA, AK, LDAK )
 
707
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, M1, ONE,
 
708
     $            DWORK( IWB ), N, DWORK( IWL ), M, ONE, AK, LDAK )
 
709
C
 
710
C     Scale Bt1 .
 
711
C
 
712
      CALL DSCAL( N*M1, GAMMA, DWORK( IWB ), 1 )
 
713
C
 
714
C     Compute |Dt11| = |R2 |*inv(V21t) .
 
715
C             |Dt21|   |D21|
 
716
C
 
717
      CALL DLACPY( 'Full', M2, M1, DWORK( IR2 ), M, DWORK( IWD ), 
 
718
     $             M2+NP2 )
 
719
      CALL DLACPY( 'Full', NP2, M1, D( NP1+1, 1 ), LDD, DWORK( IWD+M2 ),
 
720
     $             M2+NP2 )     
 
721
      CALL DTRSM( 'Right', 'Lower', 'Transpose', 'NonUnit', M2+NP2,
 
722
     $            M1, ONE, DWORK, M, DWORK( IWD ), M2+NP2 )
 
723
C
 
724
C     Compute Ct = |Ct1| = |L2| + |Dt11|*inv(V21t')*L_Nabla .
 
725
C                  |Ct2| = |C2| + |Dt21|
 
726
C
 
727
      CALL DLACPY( 'Full', M2, N, DWORK( IWL+M1 ), M, DWORK( IWC ), 
 
728
     $             M2+NP2 )
 
729
      CALL DLACPY( 'Full', NP2, N, C( NP1+1, 1 ), LDC, DWORK( IWC+M2 ),
 
730
     $             M2+NP2 )     
 
731
      CALL DGEMM( 'NoTranspose', 'NoTranspose', M2+NP2, N, M1, ONE,
 
732
     $            DWORK( IWD ), M2+NP2, DWORK( IWL ), M, ONE, 
 
733
     $            DWORK( IWC ), M2+NP2 )
 
734
C
 
735
C     Scale |Dt11| .
 
736
C           |Dt21|
 
737
C
 
738
      CALL DSCAL( ( M2+NP2 )*M1, GAMMA, DWORK( IWD ), 1 )
 
739
C
 
740
C     Workspace usage.
 
741
C
 
742
      IWW  = IWD + ( M2 + NP2 )*M1  
 
743
      IWQ  = IWW + ( M2 + NP2 )*( M2 + NP2 )
 
744
      IWL  = IWQ + N*N
 
745
      IWR  = IWL + N*( M2 + NP2 )
 
746
      IWI  = IWR + 2*N
 
747
      IWH  = IWI + 2*N
 
748
      IWS  = IWH + 2*N
 
749
      IWT  = IWS + ( 2*N + M2 + NP2 )*( 2*N + M2 + NP2 )
 
750
      IWU  = IWT + ( 2*N + M2 + NP2 )*2*N 
 
751
      IWG  = IWU + 4*N*N
 
752
      IWRK = IWG + ( M2 + NP2 )*N
 
753
      IS2  = IWW + ( M2 + NP2 )*M2
 
754
      IS3  = IS2 + M2
 
755
C
 
756
C     Compute S0 = |Dt11||Dt11' Dt21'| -|gamma^2*Im2 0| .
 
757
C                  |Dt21|               |   0        0|
 
758
C
 
759
      CALL DSYRK( 'Upper', 'NoTranspose', M2+NP2, M1, ONE, DWORK( IWD ),
 
760
     $            M2+NP2, ZERO, DWORK( IWW ), M2+NP2 )
 
761
      DO 20 J = IWW, IWW - 1 + ( M2 + NP2 )*M2, M2 + NP2 + 1
 
762
         DWORK( J ) = DWORK( J ) - GAMMA*GAMMA
 
763
   20 CONTINUE
 
764
C
 
765
C     Compute Bt1*Bt1' .
 
766
C
 
767
      CALL DSYRK( 'Upper', 'NoTranspose', N, M1, ONE, DWORK( IWB ), N,
 
768
     $            ZERO, DWORK( IWQ ), N )
 
769
C
 
770
C     Compute Bt1*|Dt11' Dt21'| .
 
771
C
 
772
      CALL DGEMM( 'NoTranspose', 'Transpose', N, M2+NP2, M1, ONE,
 
773
     $            DWORK( IWB ), N, DWORK( IWD ), M2+NP2, ZERO,
 
774
     $            DWORK( IWL ), N )
 
775
C
 
776
C     Transpose At in situ (in AK) .
 
777
C
 
778
      DO 30 J = 2, N
 
779
         CALL DSWAP( J-1, AK( J, 1 ), LDAK, AK( 1, J ), 1 )
 
780
   30 CONTINUE 
 
781
C
 
782
C     Transpose Ct .
 
783
C
 
784
      CALL MA02AD( 'Full', M2+NP2, N, DWORK( IWC ), M2+NP2, 
 
785
     $             DWORK( IWG ), N )
 
786
C
 
787
C     Solution of the Z-Riccati equation.
 
788
C     Workspace: need   13*N*N + M*M + (8*N+M+M2+2*NP2)*(M2+NP2) +
 
789
C                       N*(M+NP2) + 6*N + 
 
790
C                       max(14*N+23,16*N,2*N+M2+NP2,3*(M2+NP2));
 
791
C                prefer larger.
 
792
C
 
793
      CALL SB02OD( 'D', 'B', 'N', 'U', 'N', 'S', N, M2+NP2, NP, AK, 
 
794
     $             LDAK, DWORK( IWG ), N, DWORK( IWQ ), N, DWORK( IWW ),
 
795
     $             M2+NP2, DWORK( IWL ), N, RCOND2, Z, LDZ, DWORK( IWR),
 
796
     $             DWORK( IWI ), DWORK( IWH ), DWORK( IWS ), 2*N+M2+NP2,
 
797
     $             DWORK( IWT ), 2*N+M2+NP2, DWORK( IWU ), 2*N, TOLL,
 
798
     $             IWORK, DWORK( IWRK ), LDWORK-IWRK+1, BWORK, INFO2 )
 
799
      IF( INFO2.GT.0 ) THEN
 
800
         INFO = 7
 
801
         RETURN
 
802
      END IF
 
803
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
804
C
 
805
C     Condition estimation.
 
806
C     Workspace: need   4*N*N + M*M + 2*(M2+NP2)*(M2+NP2)+
 
807
C                       N*(M+2*M2+3*NP2) + (M2+NP2)*M1 +
 
808
C                       max(5*N,max(3,2*N*N)+N*N);
 
809
C                prefer larger.
 
810
C
 
811
      IWS  = IWR
 
812
      IWH  = IWS + ( M2 + NP2 )*( M2 + NP2 )
 
813
      IWT  = IWH + N*( M2 + NP2 )
 
814
      IWU  = IWT + N*N 
 
815
      IWG  = IWU + N*N 
 
816
      IWRK = IWG + N*N
 
817
      CALL DLACPY( 'Upper', M2+NP2, M2+NP2, DWORK( IWW ), M2+NP2,
 
818
     $             DWORK( IWS ), M2+NP2 )
 
819
      CALL DSYTRF( 'Upper', M2+NP2, DWORK( IWS ), M2+NP2, IWORK,
 
820
     $             DWORK( IWRK ), LDWORK-IWRK+1, INFO2 )
 
821
      IF( INFO2.GT.0 ) THEN
 
822
         INFO = 5
 
823
         RETURN
 
824
      END IF
 
825
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
826
C
 
827
      CALL DLACPY( 'Full', M2+NP2, N, DWORK( IWC ), M2+NP2, 
 
828
     $             DWORK( IWH ), M2+NP2 )
 
829
      CALL DSYTRS( 'Upper', M2+NP2, N, DWORK( IWS ), M2+NP2, IWORK,
 
830
     $             DWORK( IWH ), M2+NP2, INFO2 )
 
831
      CALL MB01RX( 'Left', 'Upper', 'Transpose', N, M2+NP2, ZERO, ONE,
 
832
     $             DWORK( IWG ), N, DWORK( IWC ), M2+NP2, DWORK( IWH ),
 
833
     $             M2+NP2, INFO2 )
 
834
      CALL SB02SD( 'C', 'N', 'N', 'U', 'O', N, AK, LDAK, DWORK( IWT ),
 
835
     $             N, DWORK( IWU ), N, DWORK( IWG ), N, DWORK( IWQ ), N,
 
836
     $             Z, LDZ, SEPD, RCOND( 8 ), FERR, IWORK, DWORK( IWRK ),
 
837
     $             LDWORK-IWRK+1, INFO2 )
 
838
      IF( INFO2.GT.0 ) RCOND( 8 ) = ZERO
 
839
      LWAMAX = MAX( INT( DWORK( IWRK ) ) + IWRK - 1, LWAMAX )
 
840
C
 
841
C     Workspace usage.
 
842
C
 
843
      IWRK = IWR
 
844
C
 
845
C     Compute the upper triangle of
 
846
C             |St1   St2| = S0 + |Ct1|*Z*|Ct1' Ct2'| .
 
847
C             |St2'  St3|        |Ct2|
 
848
C
 
849
      CALL MB01RU( 'Upper', 'NoTranspose', M2+NP2, N, ONE, ONE,
 
850
     $             DWORK( IWW ), M2+NP2, DWORK( IWC ), M2+NP2, Z, LDZ, 
 
851
     $             DWORK( IWRK ), (M2+NP2)*N, INFO2 )
 
852
C
 
853
C     Compute the Cholesky factorization of St3, St3 = U12'*U12 .
 
854
C
 
855
      ANORM = DLANSY( '1', 'Upper', NP2, DWORK( IS3 ), M2+NP2, 
 
856
     $                DWORK( IWRK ) )
 
857
      CALL DPOTRF( 'Upper', NP2, DWORK( IS3 ), M2+NP2, INFO2 )
 
858
      IF( INFO2.GT.0 ) THEN
 
859
         INFO = 5
 
860
         RETURN
 
861
      END IF
 
862
      CALL DPOCON( 'Upper', NP2, DWORK( IS3 ), M2+NP2, ANORM,
 
863
     $             RCOND( 4 ), DWORK( IWRK ), IWORK, INFO2 )
 
864
C
 
865
C     Return if the matrix is singular to working precision.
 
866
C
 
867
      IF( RCOND( 4 ).LT.TOLL ) THEN
 
868
         INFO = 5
 
869
         RETURN
 
870
      END IF
 
871
C
 
872
C     Compute St2 <- St2*inv(U12) .
 
873
C 
 
874
      CALL DTRSM( 'Right', 'Upper', 'NoTranspose', 'NonUnit', M2, NP2,
 
875
     $            ONE, DWORK( IS3 ), M2+NP2, DWORK( IS2 ), M2+NP2 )
 
876
C
 
877
C     Check the negative definiteness of St1 - St2*inv(St3)*St2' .
 
878
C
 
879
      CALL DSYRK( 'Upper', 'NoTranspose', M2, NP2, ONE, DWORK( IS2 ), 
 
880
     $            M2+NP2, -ONE, DWORK( IWW ), M2+NP2 )
 
881
      CALL DPOTRF( 'Upper', M2, DWORK( IWW ), M2+NP2, INFO2 )
 
882
      IF( INFO2.GT.0 ) THEN
 
883
         INFO = 5
 
884
         RETURN
 
885
      END IF     
 
886
C
 
887
C     Restore At in situ .
 
888
C
 
889
      DO 40 J = 2, N
 
890
         CALL DSWAP( J-1, AK( J, 1 ), LDAK, AK( 1, J ), 1 )
 
891
   40 CONTINUE 
 
892
C
 
893
C     Compute At*Z .
 
894
C
 
895
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, N, ONE, AK, LDAK,
 
896
     $            Z, LDZ, ZERO, DWORK( IWRK ), N )
 
897
C
 
898
C     Compute Mt2 = Bt1*Dt21' + At*Z*Ct2' in BK .
 
899
C
 
900
      CALL DLACPY( 'Full', N, NP2, DWORK( IWL+N*M2 ), N, BK, LDBK )
 
901
      CALL DGEMM( 'NoTranspose', 'Transpose', N, NP2, N, ONE,
 
902
     $            DWORK( IWRK ), N, DWORK( IWC+M2 ), M2+NP2, ONE,
 
903
     $            BK, LDBK )
 
904
C
 
905
C     Compute St2 <- St2*inv(U12') .
 
906
C
 
907
      CALL DTRSM( 'Right', 'Upper', 'Transpose', 'NonUnit', M2, NP2,
 
908
     $            ONE, DWORK( IS3 ), M2+NP2, DWORK( IS2 ), M2+NP2 )
 
909
C
 
910
C     Compute DKHAT = -inv(V12)*St2 in DK .
 
911
C
 
912
      CALL DLACPY( 'Full', M2, NP2, DWORK( IS2 ), M2+NP2, DK, LDDK ) 
 
913
      CALL DTRSM( 'Left', 'Lower', 'Transpose', 'NonUnit', M2, NP2,
 
914
     $            -ONE, DWORK( IR3 ), M, DK, LDDK )
 
915
C
 
916
C     Compute CKHAT = -inv(V12)*(Ct1 - St2*inv(St3)*Ct2) in CK .
 
917
C
 
918
      CALL DLACPY( 'Full', M2, N, DWORK( IWC ), M2+NP2, CK, LDCK ) 
 
919
      CALL DGEMM( 'NoTranspose', 'NoTranspose', M2, N, NP2, -ONE,
 
920
     $            DWORK( IS2 ), M2+NP2, DWORK( IWC+M2 ), M2+NP2, ONE, 
 
921
     $            CK, LDCK )
 
922
      CALL DTRSM( 'Left', 'Lower', 'Transpose', 'NonUnit', M2, N, -ONE,
 
923
     $            DWORK( IR3 ), M, CK, LDCK )
 
924
C
 
925
C     Compute Mt2*inv(St3) in BK .
 
926
C
 
927
      CALL DTRSM( 'Right', 'Upper', 'NoTranspose', 'NonUnit', N, NP2,
 
928
     $            ONE, DWORK( IS3 ), M2+NP2, BK, LDBK )
 
929
      CALL DTRSM( 'Right', 'Upper', 'Transpose', 'NonUnit', N, NP2,
 
930
     $            ONE, DWORK( IS3 ), M2+NP2, BK, LDBK )
 
931
C
 
932
C     Compute AKHAT in AK .
 
933
C
 
934
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, M2, ONE,
 
935
     $            B( 1, M1+1 ), LDB, CK, LDCK, ONE, AK, LDAK )
 
936
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, NP2, -ONE, BK, 
 
937
     $            LDBK, DWORK( IWC+M2 ), M2+NP2, ONE, AK, LDAK )
 
938
C
 
939
C     Compute BKHAT in BK .
 
940
C
 
941
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, NP2, M2, ONE, 
 
942
     $            B( 1, M1+1 ), LDB, DK, LDDK, ONE, BK, LDBK )
 
943
C
 
944
C     Compute Im2 + DKHAT*D22 .
 
945
C
 
946
      IWRK = M2*M2 + 1
 
947
      CALL DLASET( 'Full', M2, M2, ZERO, ONE, DWORK, M2 )
 
948
      CALL DGEMM( 'NoTranspose', 'NoTranspose', M2, M2, NP2, ONE, DK,
 
949
     $            LDDK, D( NP1+1, M1+1 ), LDD, ONE, DWORK, M2 )
 
950
      ANORM = DLANGE( '1', M2, M2, DWORK, M2, DWORK( IWRK ) )
 
951
      CALL DGETRF( M2, M2, DWORK, M2, IWORK, INFO2 )
 
952
      IF( INFO2.GT.0 ) THEN
 
953
         INFO = 8
 
954
         RETURN
 
955
      END IF
 
956
      CALL DGECON( '1', M2, DWORK, M2, ANORM, RCOND( 6 ), DWORK( IWRK ),
 
957
     $             IWORK( M2+1 ), INFO2 )
 
958
C
 
959
C     Return if the matrix is singular to working precision.
 
960
C
 
961
      IF( RCOND( 6 ).LT.TOLL ) THEN
 
962
         INFO = 8
 
963
         RETURN
 
964
      END IF
 
965
C
 
966
C     Compute CK .
 
967
C
 
968
      CALL DGETRS( 'NoTranspose', M2, N, DWORK, M2, IWORK, CK, LDCK, 
 
969
     $             INFO2 )
 
970
C
 
971
C     Compute DK .
 
972
C
 
973
      CALL DGETRS( 'NoTranspose', M2, NP2, DWORK, M2, IWORK, DK, LDDK,
 
974
     $             INFO2 )
 
975
C
 
976
C     Compute AK .
 
977
C
 
978
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, M2, NP2, ONE, BK, 
 
979
     $            LDBK, D( NP1+1, M1+1 ), LDD, ZERO, DWORK, N )
 
980
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, M2, -ONE, DWORK, 
 
981
     $            N, CK, LDCK, ONE, AK, LDAK )
 
982
C
 
983
C     Compute BK .
 
984
C 
 
985
      CALL DGEMM( 'NoTranspose', 'NoTranspose', N, NP2, M2, -ONE, DWORK,
 
986
     $            N, DK, LDDK, ONE, BK, LDBK )
 
987
C
 
988
      DWORK( 1 ) = DBLE( LWAMAX )
 
989
      RETURN
 
990
C *** Last line of SB10DD ***
 
991
      END