← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/slicot/sb03my.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/sb03my.f
 
1
      SUBROUTINE SB03MY( TRANA, N, A, LDA, C, LDC, SCALE, INFO )
 
2
C
 
3
C     RELEASE 4.0, WGS COPYRIGHT 1999.
 
4
C
 
5
C     PURPOSE
 
6
C
 
7
C     To solve the real Lyapunov matrix equation
 
8
C
 
9
C            op(A)'*X + X*op(A) = scale*C
 
10
C
 
11
C     where op(A) = A or A' (A**T), A is upper quasi-triangular and C is
 
12
C     symmetric (C = C'). (A' denotes the transpose of the matrix A.)
 
13
C     A is N-by-N, the right hand side C and the solution X are N-by-N,
 
14
C     and scale is an output scale factor, set less than or equal to 1
 
15
C     to avoid overflow in X. The solution matrix X is overwritten
 
16
C     onto C.
 
17
C
 
18
C     A must be in Schur canonical form (as returned by LAPACK routines
 
19
C     DGEES or DHSEQR), that is, block upper triangular with 1-by-1 and
 
20
C     2-by-2 diagonal blocks; each 2-by-2 diagonal block has its
 
21
C     diagonal elements equal and its off-diagonal elements of opposite
 
22
C     sign.
 
23
C
 
24
C     ARGUMENTS 
 
25
C
 
26
C     Mode Parameters
 
27
C
 
28
C     TRANA   CHARACTER*1
 
29
C             Specifies the form of op(A) to be used, as follows:
 
30
C             = 'N':  op(A) = A    (No transpose);
 
31
C             = 'T':  op(A) = A**T (Transpose);
 
32
C             = 'C':  op(A) = A**T (Conjugate transpose = Transpose).
 
33
C
 
34
C     Input/Output Parameters
 
35
C
 
36
C     N       (input) INTEGER
 
37
C             The order of the matrices A, X, and C.  N >= 0.
 
38
C
 
39
C     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
 
40
C             The leading N-by-N part of this array must contain the
 
41
C             upper quasi-triangular matrix A, in Schur canonical form.
 
42
C             The part of A below the first sub-diagonal is not
 
43
C             referenced.
 
44
C
 
45
C     LDA     INTEGER
 
46
C             The leading dimension of array A.  LDA >= MAX(1,N).
 
47
C
 
48
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 
49
C             On entry, the leading N-by-N part of this array must  
 
50
C             contain the symmetric matrix C.
 
51
C             On exit, if INFO >= 0, the leading N-by-N part of this
 
52
C             array contains the symmetric solution matrix X.
 
53
C
 
54
C     LDC     INTEGER
 
55
C             The leading dimension of array C.  LDC >= MAX(1,N).
 
56
C
 
57
C     SCALE   (output) DOUBLE PRECISION
 
58
C             The scale factor, scale, set less than or equal to 1 to
 
59
C             prevent the solution overflowing.
 
60
C
 
61
C     Error Indicator
 
62
C
 
63
C     INFO    INTEGER
 
64
C             = 0:  successful exit;
 
65
C             < 0:  if INFO = -i, the i-th argument had an illegal
 
66
C                   value;
 
67
C             = 1:  if A and -A have common or very close eigenvalues;
 
68
C                   perturbed values were used to solve the equation
 
69
C                   (but the matrix A is unchanged).
 
70
C
 
71
C     METHOD
 
72
C
 
73
C     Bartels-Stewart algorithm is used. A set of equivalent linear
 
74
C     algebraic systems of equations of order at most four are formed
 
75
C     and solved using Gaussian elimination with complete pivoting.
 
76
C
 
77
C     REFERENCES
 
78
C
 
79
C     [1] Bartels, R.H. and Stewart, G.W.  T
 
80
C         Solution of the matrix equation A X + XB = C.
 
81
C         Comm. A.C.M., 15, pp. 820-826, 1972.
 
82
C
 
83
C     NUMERICAL ASPECTS
 
84
C                               3
 
85
C     The algorithm requires 0(N ) operations.
 
86
C
 
87
C     CONTRIBUTOR
 
88
C
 
89
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, May 1997.
 
90
C     Supersedes Release 2.0 routine SB03AY by Control Systems Research
 
91
C     Group, Kingston Polytechnic, United Kingdom, October 1982.
 
92
C     Based on DTRLYP by P. Petkov, Tech. University of Sofia, September
 
93
C     1993.
 
94
C
 
95
C     REVISIONS
 
96
C
 
97
C     V. Sima, Katholieke Univ. Leuven, Belgium, May 1999.
 
98
C
 
99
C     KEYWORDS
 
100
C
 
101
C     Continuous-time system, Lyapunov equation, matrix algebra, real
 
102
C     Schur form.
 
103
C
 
104
C     ******************************************************************
 
105
C
 
106
C     .. Parameters ..
 
107
      DOUBLE PRECISION   ZERO, ONE
 
108
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 
109
C     ..
 
110
C     .. Scalar Arguments ..
 
111
      CHARACTER          TRANA
 
112
      INTEGER            INFO, LDA, LDC, N
 
113
      DOUBLE PRECISION   SCALE
 
114
C     ..
 
115
C     .. Array Arguments ..
 
116
      DOUBLE PRECISION   A( LDA, * ), C( LDC, * )
 
117
C     ..
 
118
C     .. Local Scalars ..
 
119
      LOGICAL            NOTRNA, LUPPER
 
120
      INTEGER            IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT,
 
121
     $                   MINK1N, MINK2N, MINL1N, MINL2N
 
122
      DOUBLE PRECISION   A11, BIGNUM, DA11, DB, EPS, SCALOC, SMIN,
 
123
     $                   SMLNUM, XNORM
 
124
C     ..
 
125
C     .. Local Arrays ..
 
126
      DOUBLE PRECISION   DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
 
127
C     ..
 
128
C     .. External Functions ..
 
129
      LOGICAL            LSAME
 
130
      DOUBLE PRECISION   DDOT, DLAMCH, DLANHS
 
131
      EXTERNAL           DDOT, DLAMCH, DLANHS, LSAME
 
132
C     ..
 
133
C     .. External Subroutines ..
 
134
      EXTERNAL           DLABAD, DLALN2, DLASY2, DSCAL, SB03MW, XERBLA
 
135
C     ..
 
136
C     .. Intrinsic Functions ..
 
137
      INTRINSIC          ABS, DBLE, MAX, MIN
 
138
C     ..
 
139
C     .. Executable Statements ..
 
140
C
 
141
C     Decode and Test input parameters.
 
142
C
 
143
      NOTRNA = LSAME( TRANA, 'N' )
 
144
      LUPPER = .TRUE.
 
145
C
 
146
      INFO = 0
 
147
      IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND.
 
148
     $                      .NOT.LSAME( TRANA, 'C' ) ) THEN
 
149
         INFO = -1
 
150
      ELSE IF( N.LT.0 ) THEN
 
151
         INFO = -2
 
152
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 
153
         INFO = -4
 
154
      ELSE IF( LDC.LT.MAX( 1, N ) ) THEN
 
155
         INFO = -6
 
156
      END IF
 
157
C
 
158
      IF( INFO.NE.0 ) THEN
 
159
         CALL XERBLA( 'SB03MY', -INFO )
 
160
         RETURN
 
161
      END IF
 
162
C
 
163
      SCALE = ONE
 
164
C
 
165
C     Quick return if possible.
 
166
C
 
167
      IF( N.EQ.0 )
 
168
     $   RETURN
 
169
C
 
170
C     Set constants to control overflow.
 
171
C
 
172
      EPS = DLAMCH( 'P' )
 
173
      SMLNUM = DLAMCH( 'S' )
 
174
      BIGNUM = ONE / SMLNUM
 
175
      CALL DLABAD( SMLNUM, BIGNUM )
 
176
      SMLNUM = SMLNUM*DBLE( N*N ) / EPS
 
177
      BIGNUM = ONE / SMLNUM
 
178
C
 
179
      SMIN = MAX( SMLNUM, EPS*DLANHS( 'Max', N, A, LDA, DUM ) )
 
180
C
 
181
      IF( NOTRNA ) THEN
 
182
C
 
183
C        Solve    A'*X + X*A = scale*C.
 
184
C
 
185
C        The (K,L)th block of X is determined starting from
 
186
C        upper-left corner column by column by
 
187
C
 
188
C          A(K,K)'*X(K,L) + X(K,L)*A(L,L) = C(K,L) - R(K,L),
 
189
C
 
190
C        where
 
191
C                   K-1                    L-1
 
192
C          R(K,L) = SUM [A(I,K)'*X(I,L)] + SUM [X(K,J)*A(J,L)].
 
193
C                   I=1                    J=1
 
194
C
 
195
C        Start column loop (index = L).
 
196
C        L1 (L2): column index of the first (last) row of X(K,L).
 
197
C
 
198
         LNEXT = 1
 
199
C
 
200
         DO 60 L = 1, N
 
201
            IF( L.LT.LNEXT )
 
202
     $         GO TO 60
 
203
            L1 = L
 
204
            L2 = L
 
205
            IF( L.LT.N ) THEN
 
206
               IF( A( L+1, L ).NE.ZERO )
 
207
     $            L2 = L2 + 1
 
208
               LNEXT = L2 + 1
 
209
            END IF
 
210
C
 
211
C           Start row loop (index = K).
 
212
C           K1 (K2): row index of the first (last) row of X(K,L).
 
213
C
 
214
            KNEXT = L
 
215
C
 
216
            DO 50 K = L, N
 
217
               IF( K.LT.KNEXT )
 
218
     $            GO TO 50
 
219
               K1 = K
 
220
               K2 = K
 
221
               IF( K.LT.N ) THEN
 
222
                  IF( A( K+1, K ).NE.ZERO )
 
223
     $               K2 = K2 + 1
 
224
                  KNEXT = K2 + 1
 
225
               END IF
 
226
C
 
227
               IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
 
228
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
229
     $               ( DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) +
 
230
     $                 DDOT( L1-1, C( K1, 1 ), LDC, A( 1, L1 ), 1 ) )
 
231
                  SCALOC = ONE
 
232
C
 
233
                  A11 = A( K1, K1 ) + A( L1, L1 )
 
234
                  DA11 = ABS( A11 )
 
235
                  IF( DA11.LE.SMIN ) THEN
 
236
                     A11 = SMIN
 
237
                     DA11 = SMIN
 
238
                     INFO = 1
 
239
                  END IF
 
240
                  DB = ABS( VEC( 1, 1 ) )
 
241
                  IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
 
242
                     IF( DB.GT.BIGNUM*DA11 )
 
243
     $                  SCALOC = ONE / DB
 
244
                  END IF
 
245
                  X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
 
246
C
 
247
                  IF( SCALOC.NE.ONE ) THEN
 
248
C
 
249
                     DO 10 J = 1, N
 
250
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
251
   10                CONTINUE
 
252
C
 
253
                     SCALE = SCALE*SCALOC
 
254
                  END IF
 
255
                  C( K1, L1 ) = X( 1, 1 )
 
256
                  IF( K1.NE.L1 ) THEN
 
257
                     C( L1, K1 ) = X( 1, 1 )
 
258
                  END IF  
 
259
C
 
260
               ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
 
261
C
 
262
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
263
     $               ( DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) +
 
264
     $                 DDOT( L1-1, C( K1, 1 ), LDC, A( 1, L1 ), 1 ) )
 
265
C
 
266
                  VEC( 2, 1 ) = C( K2, L1 ) -
 
267
     $               ( DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) +
 
268
     $                 DDOT( L1-1, C( K2, 1 ), LDC, A( 1, L1 ), 1 ) )
 
269
C
 
270
                  CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
 
271
     $                         LDA, ONE, ONE, VEC, 2, -A( L1, L1 ),
 
272
     $                         ZERO, X, 2, SCALOC, XNORM, IERR )
 
273
                  IF( IERR.NE.0 )
 
274
     $               INFO = 1
 
275
C
 
276
                  IF( SCALOC.NE.ONE ) THEN
 
277
C
 
278
                     DO 20 J = 1, N
 
279
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
280
   20                CONTINUE
 
281
C
 
282
                     SCALE = SCALE*SCALOC
 
283
                  END IF
 
284
                  C( K1, L1 ) = X( 1, 1 )
 
285
                  C( K2, L1 ) = X( 2, 1 )
 
286
                  C( L1, K1 ) = X( 1, 1 )
 
287
                  C( L1, K2 ) = X( 2, 1 )  
 
288
C
 
289
               ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
 
290
C
 
291
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
292
     $               ( DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) +
 
293
     $                 DDOT( L1-1, C( K1, 1 ), LDC, A( 1, L1 ), 1 ) )
 
294
C
 
295
                  VEC( 2, 1 ) = C( K1, L2 ) -
 
296
     $               ( DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) +
 
297
     $                 DDOT( L1-1, C( K1, 1 ), LDC, A( 1, L2 ), 1 ) )
 
298
C
 
299
                  CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( L1, L1 ),
 
300
     $                         LDA, ONE, ONE, VEC, 2, -A( K1, K1 ),
 
301
     $                         ZERO, X, 2, SCALOC, XNORM, IERR )
 
302
                  IF( IERR.NE.0 )
 
303
     $               INFO = 1
 
304
C
 
305
                  IF( SCALOC.NE.ONE ) THEN
 
306
C
 
307
                     DO 30 J = 1, N
 
308
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
309
   30                CONTINUE
 
310
C
 
311
                     SCALE = SCALE*SCALOC
 
312
                  END IF
 
313
                  C( K1, L1 ) = X( 1, 1 )
 
314
                  C( K1, L2 ) = X( 2, 1 )
 
315
                  C( L1, K1 ) = X( 1, 1 )
 
316
                  C( L2, K1 ) = X( 2, 1 )
 
317
C
 
318
               ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
 
319
C
 
320
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
321
     $               ( DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) +
 
322
     $                 DDOT( L1-1, C( K1, 1 ), LDC, A( 1, L1 ), 1 ) )
 
323
C
 
324
                  VEC( 1, 2 ) = C( K1, L2 ) -
 
325
     $               ( DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) +
 
326
     $                 DDOT( L1-1, C( K1, 1 ), LDC, A( 1, L2 ), 1 ) )
 
327
C                 
 
328
                  VEC( 2, 1 ) = C( K2, L1 ) -
 
329
     $               ( DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) +
 
330
     $                 DDOT( L1-1, C( K2, 1 ), LDC, A( 1, L1 ), 1 ) )
 
331
C
 
332
                  VEC( 2, 2 ) = C( K2, L2 ) -
 
333
     $               ( DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) +
 
334
     $                 DDOT( L1-1, C( K2, 1 ), LDC, A( 1, L2 ), 1 ) )
 
335
C
 
336
                  IF( K1.EQ.L1 ) THEN
 
337
                     CALL SB03MW( .FALSE., LUPPER, A( K1, K1 ), LDA,
 
338
     $                            VEC, 2, SCALOC, X, 2, XNORM, IERR )
 
339
                     IF( LUPPER ) THEN
 
340
                        X( 2, 1 ) = X( 1, 2 ) 
 
341
                     ELSE
 
342
                        X( 1, 2 ) = X( 2, 1 ) 
 
343
                     END IF
 
344
                  ELSE    
 
345
                     CALL DLASY2( .TRUE., .FALSE., 1, 2, 2, A( K1, K1 ),
 
346
     $                            LDA, A( L1, L1 ), LDA, VEC, 2, SCALOC,
 
347
     $                            X, 2, XNORM, IERR )
 
348
                  END IF 
 
349
                  IF( IERR.NE.0 )
 
350
     $               INFO = 1
 
351
C
 
352
                  IF( SCALOC.NE.ONE ) THEN
 
353
C
 
354
                     DO 40 J = 1, N
 
355
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
356
   40                CONTINUE
 
357
C
 
358
                     SCALE = SCALE*SCALOC
 
359
                  END IF
 
360
                  C( K1, L1 ) = X( 1, 1 )
 
361
                  C( K1, L2 ) = X( 1, 2 )
 
362
                  C( K2, L1 ) = X( 2, 1 )
 
363
                  C( K2, L2 ) = X( 2, 2 )
 
364
                  IF( K1.NE.L1 ) THEN
 
365
                     C( L1, K1 ) = X( 1, 1 )
 
366
                     C( L2, K1 ) = X( 1, 2 )
 
367
                     C( L1, K2 ) = X( 2, 1 )
 
368
                     C( L2, K2 ) = X( 2, 2 )
 
369
                  END IF   
 
370
               END IF
 
371
C
 
372
   50       CONTINUE
 
373
C
 
374
   60    CONTINUE
 
375
C
 
376
      ELSE
 
377
C
 
378
C        Solve    A*X + X*A' = scale*C.
 
379
C
 
380
C        The (K,L)th block of X is determined starting from
 
381
C        bottom-right corner column by column by
 
382
C
 
383
C            A(K,K)*X(K,L) + X(K,L)*A(L,L)' = C(K,L) - R(K,L),
 
384
C
 
385
C        where
 
386
C                      N                     N
 
387
C            R(K,L) = SUM [A(K,I)*X(I,L)] + SUM [X(K,J)*A(L,J)'].
 
388
C                    I=K+1                 J=L+1
 
389
C
 
390
C        Start column loop (index = L).
 
391
C        L1 (L2): column index of the first (last) row of X(K,L).
 
392
C
 
393
         LNEXT = N
 
394
C
 
395
         DO 120 L = N, 1, -1
 
396
            IF( L.GT.LNEXT )
 
397
     $         GO TO 120
 
398
            L1 = L
 
399
            L2 = L
 
400
            IF( L.GT.1 ) THEN
 
401
               IF( A( L, L-1 ).NE.ZERO )
 
402
     $            L1 = L1 - 1
 
403
               LNEXT = L1 - 1
 
404
            END IF
 
405
            MINL1N = MIN( L1+1, N )
 
406
            MINL2N = MIN( L2+1, N )
 
407
C
 
408
C           Start row loop (index = K).
 
409
C           K1 (K2): row index of the first (last) row of X(K,L).
 
410
C
 
411
            KNEXT = L
 
412
C
 
413
            DO 110 K = L, 1, -1
 
414
               IF( K.GT.KNEXT )
 
415
     $            GO TO 110
 
416
               K1 = K
 
417
               K2 = K
 
418
               IF( K.GT.1 ) THEN
 
419
                  IF( A( K, K-1 ).NE.ZERO )
 
420
     $               K1 = K1 - 1
 
421
                  KNEXT = K1 - 1
 
422
               END IF
 
423
               MINK1N = MIN( K1+1, N )
 
424
               MINK2N = MIN( K2+1, N )
 
425
C
 
426
               IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
 
427
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
428
     $               ( DDOT( N-K1, A( K1, MINK1N ), LDA,
 
429
     $                       C( MINK1N, L1 ), 1 ) +
 
430
     $                 DDOT( N-L1, C( K1, MINL1N ), LDC,
 
431
     $                       A( L1, MINL1N ), LDA ) )
 
432
                  SCALOC = ONE
 
433
C
 
434
                  A11 = A( K1, K1 ) + A( L1, L1 )
 
435
                  DA11 = ABS( A11 )
 
436
                  IF( DA11.LE.SMIN ) THEN
 
437
                     A11 = SMIN
 
438
                     DA11 = SMIN
 
439
                     INFO = 1
 
440
                  END IF
 
441
                  DB = ABS( VEC( 1, 1 ) )
 
442
                  IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
 
443
                     IF( DB.GT.BIGNUM*DA11 )
 
444
     $                  SCALOC = ONE / DB
 
445
                  END IF
 
446
                  X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
 
447
C
 
448
                  IF( SCALOC.NE.ONE ) THEN
 
449
C
 
450
                     DO 70 J = 1, N
 
451
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
452
   70                CONTINUE
 
453
C
 
454
                     SCALE = SCALE*SCALOC
 
455
                  END IF
 
456
                  C( K1, L1 ) = X( 1, 1 )
 
457
                  IF( K1.NE.L1 ) THEN
 
458
                     C( L1, K1 ) = X( 1, 1 )
 
459
                  END IF   
 
460
C
 
461
               ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
 
462
C
 
463
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
464
     $               ( DDOT( N-K2, A( K1, MINK2N ), LDA,
 
465
     $                       C( MINK2N, L1 ), 1 ) +
 
466
     $                 DDOT( N-L2, C( K1, MINL2N ), LDC,
 
467
     $                       A( L1, MINL2N ), LDA ) )
 
468
C
 
469
                  VEC( 2, 1 ) = C( K2, L1 ) -
 
470
     $               ( DDOT( N-K2, A( K2, MINK2N ), LDA,
 
471
     $                     C( MINK2N, L1 ), 1 ) +
 
472
     $                 DDOT( N-L2, C( K2, MINL2N ), LDC,
 
473
     $                     A( L1, MINL2N ), LDA ) )
 
474
C
 
475
                  CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
 
476
     $                         LDA, ONE, ONE, VEC, 2, -A( L1, L1 ),
 
477
     $                         ZERO, X, 2, SCALOC, XNORM, IERR )
 
478
                  IF( IERR.NE.0 )
 
479
     $               INFO = 1
 
480
C
 
481
                  IF( SCALOC.NE.ONE ) THEN
 
482
C
 
483
                     DO 80 J = 1, N
 
484
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
485
   80                CONTINUE
 
486
C
 
487
                     SCALE = SCALE*SCALOC
 
488
                  END IF
 
489
                  C( K1, L1 ) = X( 1, 1 )
 
490
                  C( K2, L1 ) = X( 2, 1 )
 
491
                  C( L1, K1 ) = X( 1, 1 )
 
492
                  C( L1, K2 ) = X( 2, 1 ) 
 
493
C
 
494
               ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
 
495
C
 
496
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
497
     $               ( DDOT( N-K1, A( K1, MINK1N ), LDA,
 
498
     $                       C( MINK1N, L1 ), 1 ) +
 
499
     $                 DDOT( N-L2, C( K1, MINL2N ), LDC,
 
500
     $                       A( L1, MINL2N ), LDA ) )
 
501
C
 
502
                  VEC( 2, 1 ) = C( K1, L2 ) -
 
503
     $               ( DDOT( N-K1, A( K1, MINK1N ), LDA,
 
504
     $                       C( MINK1N, L2 ), 1 ) +
 
505
     $                 DDOT( N-L2, C( K1, MINL2N ), LDC,
 
506
     $                       A( L2, MINL2N ), LDA ) )
 
507
C
 
508
                  CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( L1, L1 ),
 
509
     $                         LDA, ONE, ONE, VEC, 2, -A( K1, K1 ),
 
510
     $                         ZERO, X, 2, SCALOC, XNORM, IERR )
 
511
                  IF( IERR.NE.0 )
 
512
     $               INFO = 1
 
513
C
 
514
                  IF( SCALOC.NE.ONE ) THEN
 
515
C
 
516
                     DO 90 J = 1, N
 
517
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
518
   90                CONTINUE
 
519
C
 
520
                     SCALE = SCALE*SCALOC
 
521
                  END IF
 
522
                  C( K1, L1 ) = X( 1, 1 )
 
523
                  C( K1, L2 ) = X( 2, 1 )
 
524
                  C( L1, K1 ) = X( 1, 1 )
 
525
                  C( L2, K1 ) = X( 2, 1 )
 
526
C
 
527
               ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
 
528
C
 
529
                  VEC( 1, 1 ) = C( K1, L1 ) -
 
530
     $               ( DDOT( N-K2, A( K1, MINK2N ), LDA,
 
531
     $                       C( MINK2N, L1 ), 1 ) +
 
532
     $                 DDOT( N-L2, C( K1, MINL2N ), LDC,
 
533
     $                       A( L1, MINL2N ), LDA ) )
 
534
C
 
535
                  VEC( 1, 2 ) = C( K1, L2 ) -
 
536
     $               ( DDOT( N-K2, A( K1, MINK2N ), LDA,
 
537
     $                       C( MINK2N, L2 ), 1 ) +
 
538
     $                 DDOT( N-L2, C( K1, MINL2N ), LDC,
 
539
     $                       A( L2, MINL2N ), LDA ) )
 
540
C                 
 
541
                  VEC( 2, 1 ) = C( K2, L1 ) -
 
542
     $               ( DDOT( N-K2, A( K2, MINK2N ), LDA,
 
543
     $                       C( MINK2N, L1 ), 1 ) +
 
544
     $                 DDOT( N-L2, C( K2, MINL2N ), LDC,
 
545
     $                       A( L1, MINL2N ), LDA ) )
 
546
C
 
547
                  VEC( 2, 2 ) = C( K2, L2 ) -
 
548
     $               ( DDOT( N-K2, A( K2, MINK2N ), LDA,
 
549
     $                       C( MINK2N, L2 ), 1 ) +
 
550
     $                 DDOT( N-L2, C( K2, MINL2N ), LDC,
 
551
     $                       A( L2, MINL2N ), LDA ) )
 
552
C
 
553
                  IF( K1.EQ.L1 ) THEN
 
554
                     CALL SB03MW( .TRUE., LUPPER, A( K1, K1 ), LDA, VEC,
 
555
     $                            2, SCALOC, X, 2, XNORM, IERR )
 
556
                     IF( LUPPER ) THEN
 
557
                        X( 2, 1 ) = X( 1, 2 ) 
 
558
                     ELSE
 
559
                        X( 1, 2 ) = X( 2, 1 ) 
 
560
                     END IF
 
561
                  ELSE  
 
562
                     CALL DLASY2( .FALSE., .TRUE., 1, 2, 2, A( K1, K1 ),
 
563
     $                            LDA, A( L1, L1 ), LDA, VEC, 2, SCALOC,
 
564
     $                            X, 2, XNORM, IERR )
 
565
                  END IF 
 
566
                  IF( IERR.NE.0 )
 
567
     $               INFO = 1
 
568
C
 
569
                  IF( SCALOC.NE.ONE ) THEN
 
570
C
 
571
                     DO 100 J = 1, N
 
572
                        CALL DSCAL( N, SCALOC, C( 1, J ), 1 )
 
573
  100                CONTINUE
 
574
C
 
575
                     SCALE = SCALE*SCALOC
 
576
                  END IF
 
577
                  C( K1, L1 ) = X( 1, 1 )
 
578
                  C( K1, L2 ) = X( 1, 2 )
 
579
                  C( K2, L1 ) = X( 2, 1 )
 
580
                  C( K2, L2 ) = X( 2, 2 )
 
581
                  IF( K1.NE.L1 ) THEN
 
582
                     C( L1, K1 ) = X( 1, 1 )
 
583
                     C( L2, K1 ) = X( 1, 2 )
 
584
                     C( L1, K2 ) = X( 2, 1 )
 
585
                     C( L2, K2 ) = X( 2, 2 )
 
586
                  END IF  
 
587
               END IF
 
588
C
 
589
  110       CONTINUE
 
590
C
 
591
  120    CONTINUE
 
592
C
 
593
      END IF
 
594
C
 
595
      RETURN
 
596
C *** Last line of SB03MY ***
 
597
      END