← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/lapack/dlasy2.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/dlasy2.f
 
1
      SUBROUTINE DLASY2( LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR,
 
2
     $                   LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO )
 
3
*
 
4
*  -- LAPACK auxiliary routine (version 2.0) --
 
5
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
6
*     Courant Institute, Argonne National Lab, and Rice University
 
7
*     October 31, 1992
 
8
*
 
9
*     .. Scalar Arguments ..
 
10
      LOGICAL            LTRANL, LTRANR
 
11
      INTEGER            INFO, ISGN, LDB, LDTL, LDTR, LDX, N1, N2
 
12
      DOUBLE PRECISION   SCALE, XNORM
 
13
*     ..
 
14
*     .. Array Arguments ..
 
15
      DOUBLE PRECISION   B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
 
16
     $                   X( LDX, * )
 
17
*     ..
 
18
*
 
19
*  Purpose
 
20
*  =======
 
21
*
 
22
*  DLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in
 
23
*
 
24
*         op(TL)*X + ISGN*X*op(TR) = SCALE*B,
 
25
*
 
26
*  where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
 
27
*  -1.  op(T) = T or T', where T' denotes the transpose of T.
 
28
*
 
29
*  Arguments
 
30
*  =========
 
31
*
 
32
*  LTRANL  (input) LOGICAL
 
33
*          On entry, LTRANL specifies the op(TL):
 
34
*             = .FALSE., op(TL) = TL,
 
35
*             = .TRUE., op(TL) = TL'.
 
36
*
 
37
*  LTRANR  (input) LOGICAL
 
38
*          On entry, LTRANR specifies the op(TR):
 
39
*            = .FALSE., op(TR) = TR,
 
40
*            = .TRUE., op(TR) = TR'.
 
41
*
 
42
*  ISGN    (input) INTEGER
 
43
*          On entry, ISGN specifies the sign of the equation
 
44
*          as described before. ISGN may only be 1 or -1.
 
45
*
 
46
*  N1      (input) INTEGER
 
47
*          On entry, N1 specifies the order of matrix TL.
 
48
*          N1 may only be 0, 1 or 2.
 
49
*
 
50
*  N2      (input) INTEGER
 
51
*          On entry, N2 specifies the order of matrix TR.
 
52
*          N2 may only be 0, 1 or 2.
 
53
*
 
54
*  TL      (input) DOUBLE PRECISION array, dimension (LDTL,2)
 
55
*          On entry, TL contains an N1 by N1 matrix.
 
56
*
 
57
*  LDTL    (input) INTEGER
 
58
*          The leading dimension of the matrix TL. LDTL >= max(1,N1).
 
59
*
 
60
*  TR      (input) DOUBLE PRECISION array, dimension (LDTR,2)
 
61
*          On entry, TR contains an N2 by N2 matrix.
 
62
*
 
63
*  LDTR    (input) INTEGER
 
64
*          The leading dimension of the matrix TR. LDTR >= max(1,N2).
 
65
*
 
66
*  B       (input) DOUBLE PRECISION array, dimension (LDB,2)
 
67
*          On entry, the N1 by N2 matrix B contains the right-hand
 
68
*          side of the equation.
 
69
*
 
70
*  LDB     (input) INTEGER
 
71
*          The leading dimension of the matrix B. LDB >= max(1,N1).
 
72
*
 
73
*  SCALE   (output) DOUBLE PRECISION
 
74
*          On exit, SCALE contains the scale factor. SCALE is chosen
 
75
*          less than or equal to 1 to prevent the solution overflowing.
 
76
*
 
77
*  X       (output) DOUBLE PRECISION array, dimension (LDX,2)
 
78
*          On exit, X contains the N1 by N2 solution.
 
79
*
 
80
*  LDX     (input) INTEGER
 
81
*          The leading dimension of the matrix X. LDX >= max(1,N1).
 
82
*
 
83
*  XNORM   (output) DOUBLE PRECISION
 
84
*          On exit, XNORM is the infinity-norm of the solution.
 
85
*
 
86
*  INFO    (output) INTEGER
 
87
*          On exit, INFO is set to
 
88
*             0: successful exit.
 
89
*             1: TL and TR have too close eigenvalues, so TL or
 
90
*                TR is perturbed to get a nonsingular equation.
 
91
*          NOTE: In the interests of speed, this routine does not
 
92
*                check the inputs for errors.
 
93
*
 
94
* =====================================================================
 
95
*
 
96
*     .. Parameters ..
 
97
      DOUBLE PRECISION   ZERO, ONE
 
98
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 
99
      DOUBLE PRECISION   TWO, HALF, EIGHT
 
100
      PARAMETER          ( TWO = 2.0D+0, HALF = 0.5D+0, EIGHT = 8.0D+0 )
 
101
*     ..
 
102
*     .. Local Scalars ..
 
103
      LOGICAL            BSWAP, XSWAP
 
104
      INTEGER            I, IP, IPIV, IPSV, J, JP, JPSV, K
 
105
      DOUBLE PRECISION   BET, EPS, GAM, L21, SGN, SMIN, SMLNUM, TAU1,
 
106
     $                   TEMP, U11, U12, U22, XMAX
 
107
*     ..
 
108
*     .. Local Arrays ..
 
109
      LOGICAL            BSWPIV( 4 ), XSWPIV( 4 )
 
110
      INTEGER            JPIV( 4 ), LOCL21( 4 ), LOCU12( 4 ),
 
111
     $                   LOCU22( 4 )
 
112
      DOUBLE PRECISION   BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
 
113
*     ..
 
114
*     .. External Functions ..
 
115
      INTEGER            IDAMAX
 
116
      DOUBLE PRECISION   DLAMCH
 
117
      EXTERNAL           IDAMAX, DLAMCH
 
118
*     ..
 
119
*     .. External Subroutines ..
 
120
      EXTERNAL           DCOPY, DSWAP
 
121
*     ..
 
122
*     .. Intrinsic Functions ..
 
123
      INTRINSIC          ABS, MAX
 
124
*     ..
 
125
*     .. Data statements ..
 
126
      DATA               LOCU12 / 3, 4, 1, 2 / , LOCL21 / 2, 1, 4, 3 / ,
 
127
     $                   LOCU22 / 4, 3, 2, 1 /
 
128
      DATA               XSWPIV / .FALSE., .FALSE., .TRUE., .TRUE. /
 
129
      DATA               BSWPIV / .FALSE., .TRUE., .FALSE., .TRUE. /
 
130
*     ..
 
131
*     .. Executable Statements ..
 
132
*
 
133
*     Do not check the input parameters for errors
 
134
*
 
135
      INFO = 0
 
136
*
 
137
*     Quick return if possible
 
138
*
 
139
      IF( N1.EQ.0 .OR. N2.EQ.0 )
 
140
     $   RETURN
 
141
*
 
142
*     Set constants to control overflow
 
143
*
 
144
      EPS = DLAMCH( 'P' )
 
145
      SMLNUM = DLAMCH( 'S' ) / EPS
 
146
      SGN = ISGN
 
147
*
 
148
      K = N1 + N1 + N2 - 2
 
149
      GO TO ( 10, 20, 30, 50 )K
 
150
*
 
151
*     1 by 1: TL11*X + SGN*X*TR11 = B11
 
152
*
 
153
   10 CONTINUE
 
154
      TAU1 = TL( 1, 1 ) + SGN*TR( 1, 1 )
 
155
      BET = ABS( TAU1 )
 
156
      IF( BET.LE.SMLNUM ) THEN
 
157
         TAU1 = SMLNUM
 
158
         BET = SMLNUM
 
159
         INFO = 1
 
160
      END IF
 
161
*
 
162
      SCALE = ONE
 
163
      GAM = ABS( B( 1, 1 ) )
 
164
      IF( SMLNUM*GAM.GT.BET )
 
165
     $   SCALE = ONE / GAM
 
166
*
 
167
      X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / TAU1
 
168
      XNORM = ABS( X( 1, 1 ) )
 
169
      RETURN
 
170
*
 
171
*     1 by 2:
 
172
*     TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12]  = [B11 B12]
 
173
*                                       [TR21 TR22]
 
174
*
 
175
   20 CONTINUE
 
176
*
 
177
      SMIN = MAX( EPS*MAX( ABS( TL( 1, 1 ) ), ABS( TR( 1, 1 ) ),
 
178
     $       ABS( TR( 1, 2 ) ), ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) ),
 
179
     $       SMLNUM )
 
180
      TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
 
181
      TMP( 4 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
 
182
      IF( LTRANR ) THEN
 
183
         TMP( 2 ) = SGN*TR( 2, 1 )
 
184
         TMP( 3 ) = SGN*TR( 1, 2 )
 
185
      ELSE
 
186
         TMP( 2 ) = SGN*TR( 1, 2 )
 
187
         TMP( 3 ) = SGN*TR( 2, 1 )
 
188
      END IF
 
189
      BTMP( 1 ) = B( 1, 1 )
 
190
      BTMP( 2 ) = B( 1, 2 )
 
191
      GO TO 40
 
192
*
 
193
*     2 by 1:
 
194
*          op[TL11 TL12]*[X11] + ISGN* [X11]*TR11  = [B11]
 
195
*            [TL21 TL22] [X21]         [X21]         [B21]
 
196
*
 
197
   30 CONTINUE
 
198
      SMIN = MAX( EPS*MAX( ABS( TR( 1, 1 ) ), ABS( TL( 1, 1 ) ),
 
199
     $       ABS( TL( 1, 2 ) ), ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) ),
 
200
     $       SMLNUM )
 
201
      TMP( 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
 
202
      TMP( 4 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
 
203
      IF( LTRANL ) THEN
 
204
         TMP( 2 ) = TL( 1, 2 )
 
205
         TMP( 3 ) = TL( 2, 1 )
 
206
      ELSE
 
207
         TMP( 2 ) = TL( 2, 1 )
 
208
         TMP( 3 ) = TL( 1, 2 )
 
209
      END IF
 
210
      BTMP( 1 ) = B( 1, 1 )
 
211
      BTMP( 2 ) = B( 2, 1 )
 
212
   40 CONTINUE
 
213
*
 
214
*     Solve 2 by 2 system using complete pivoting.
 
215
*     Set pivots less than SMIN to SMIN.
 
216
*
 
217
      IPIV = IDAMAX( 4, TMP, 1 )
 
218
      U11 = TMP( IPIV )
 
219
      IF( ABS( U11 ).LE.SMIN ) THEN
 
220
         INFO = 1
 
221
         U11 = SMIN
 
222
      END IF
 
223
      U12 = TMP( LOCU12( IPIV ) )
 
224
      L21 = TMP( LOCL21( IPIV ) ) / U11
 
225
      U22 = TMP( LOCU22( IPIV ) ) - U12*L21
 
226
      XSWAP = XSWPIV( IPIV )
 
227
      BSWAP = BSWPIV( IPIV )
 
228
      IF( ABS( U22 ).LE.SMIN ) THEN
 
229
         INFO = 1
 
230
         U22 = SMIN
 
231
      END IF
 
232
      IF( BSWAP ) THEN
 
233
         TEMP = BTMP( 2 )
 
234
         BTMP( 2 ) = BTMP( 1 ) - L21*TEMP
 
235
         BTMP( 1 ) = TEMP
 
236
      ELSE
 
237
         BTMP( 2 ) = BTMP( 2 ) - L21*BTMP( 1 )
 
238
      END IF
 
239
      SCALE = ONE
 
240
      IF( ( TWO*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( U22 ) .OR.
 
241
     $    ( TWO*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( U11 ) ) THEN
 
242
         SCALE = HALF / MAX( ABS( BTMP( 1 ) ), ABS( BTMP( 2 ) ) )
 
243
         BTMP( 1 ) = BTMP( 1 )*SCALE
 
244
         BTMP( 2 ) = BTMP( 2 )*SCALE
 
245
      END IF
 
246
      X2( 2 ) = BTMP( 2 ) / U22
 
247
      X2( 1 ) = BTMP( 1 ) / U11 - ( U12 / U11 )*X2( 2 )
 
248
      IF( XSWAP ) THEN
 
249
         TEMP = X2( 2 )
 
250
         X2( 2 ) = X2( 1 )
 
251
         X2( 1 ) = TEMP
 
252
      END IF
 
253
      X( 1, 1 ) = X2( 1 )
 
254
      IF( N1.EQ.1 ) THEN
 
255
         X( 1, 2 ) = X2( 2 )
 
256
         XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
 
257
      ELSE
 
258
         X( 2, 1 ) = X2( 2 )
 
259
         XNORM = MAX( ABS( X( 1, 1 ) ), ABS( X( 2, 1 ) ) )
 
260
      END IF
 
261
      RETURN
 
262
*
 
263
*     2 by 2:
 
264
*     op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
 
265
*       [TL21 TL22] [X21 X22]        [X21 X22]   [TR21 TR22]   [B21 B22]
 
266
*
 
267
*     Solve equivalent 4 by 4 system using complete pivoting.
 
268
*     Set pivots less than SMIN to SMIN.
 
269
*
 
270
   50 CONTINUE
 
271
      SMIN = MAX( ABS( TR( 1, 1 ) ), ABS( TR( 1, 2 ) ),
 
272
     $       ABS( TR( 2, 1 ) ), ABS( TR( 2, 2 ) ) )
 
273
      SMIN = MAX( SMIN, ABS( TL( 1, 1 ) ), ABS( TL( 1, 2 ) ),
 
274
     $       ABS( TL( 2, 1 ) ), ABS( TL( 2, 2 ) ) )
 
275
      SMIN = MAX( EPS*SMIN, SMLNUM )
 
276
      BTMP( 1 ) = ZERO
 
277
      CALL DCOPY( 16, BTMP, 0, T16, 1 )
 
278
      T16( 1, 1 ) = TL( 1, 1 ) + SGN*TR( 1, 1 )
 
279
      T16( 2, 2 ) = TL( 2, 2 ) + SGN*TR( 1, 1 )
 
280
      T16( 3, 3 ) = TL( 1, 1 ) + SGN*TR( 2, 2 )
 
281
      T16( 4, 4 ) = TL( 2, 2 ) + SGN*TR( 2, 2 )
 
282
      IF( LTRANL ) THEN
 
283
         T16( 1, 2 ) = TL( 2, 1 )
 
284
         T16( 2, 1 ) = TL( 1, 2 )
 
285
         T16( 3, 4 ) = TL( 2, 1 )
 
286
         T16( 4, 3 ) = TL( 1, 2 )
 
287
      ELSE
 
288
         T16( 1, 2 ) = TL( 1, 2 )
 
289
         T16( 2, 1 ) = TL( 2, 1 )
 
290
         T16( 3, 4 ) = TL( 1, 2 )
 
291
         T16( 4, 3 ) = TL( 2, 1 )
 
292
      END IF
 
293
      IF( LTRANR ) THEN
 
294
         T16( 1, 3 ) = SGN*TR( 1, 2 )
 
295
         T16( 2, 4 ) = SGN*TR( 1, 2 )
 
296
         T16( 3, 1 ) = SGN*TR( 2, 1 )
 
297
         T16( 4, 2 ) = SGN*TR( 2, 1 )
 
298
      ELSE
 
299
         T16( 1, 3 ) = SGN*TR( 2, 1 )
 
300
         T16( 2, 4 ) = SGN*TR( 2, 1 )
 
301
         T16( 3, 1 ) = SGN*TR( 1, 2 )
 
302
         T16( 4, 2 ) = SGN*TR( 1, 2 )
 
303
      END IF
 
304
      BTMP( 1 ) = B( 1, 1 )
 
305
      BTMP( 2 ) = B( 2, 1 )
 
306
      BTMP( 3 ) = B( 1, 2 )
 
307
      BTMP( 4 ) = B( 2, 2 )
 
308
*
 
309
*     Perform elimination
 
310
*
 
311
      DO 100 I = 1, 3
 
312
         XMAX = ZERO
 
313
         DO 70 IP = I, 4
 
314
            DO 60 JP = I, 4
 
315
               IF( ABS( T16( IP, JP ) ).GE.XMAX ) THEN
 
316
                  XMAX = ABS( T16( IP, JP ) )
 
317
                  IPSV = IP
 
318
                  JPSV = JP
 
319
               END IF
 
320
   60       CONTINUE
 
321
   70    CONTINUE
 
322
         IF( IPSV.NE.I ) THEN
 
323
            CALL DSWAP( 4, T16( IPSV, 1 ), 4, T16( I, 1 ), 4 )
 
324
            TEMP = BTMP( I )
 
325
            BTMP( I ) = BTMP( IPSV )
 
326
            BTMP( IPSV ) = TEMP
 
327
         END IF
 
328
         IF( JPSV.NE.I )
 
329
     $      CALL DSWAP( 4, T16( 1, JPSV ), 1, T16( 1, I ), 1 )
 
330
         JPIV( I ) = JPSV
 
331
         IF( ABS( T16( I, I ) ).LT.SMIN ) THEN
 
332
            INFO = 1
 
333
            T16( I, I ) = SMIN
 
334
         END IF
 
335
         DO 90 J = I + 1, 4
 
336
            T16( J, I ) = T16( J, I ) / T16( I, I )
 
337
            BTMP( J ) = BTMP( J ) - T16( J, I )*BTMP( I )
 
338
            DO 80 K = I + 1, 4
 
339
               T16( J, K ) = T16( J, K ) - T16( J, I )*T16( I, K )
 
340
   80       CONTINUE
 
341
   90    CONTINUE
 
342
  100 CONTINUE
 
343
      IF( ABS( T16( 4, 4 ) ).LT.SMIN )
 
344
     $   T16( 4, 4 ) = SMIN
 
345
      SCALE = ONE
 
346
      IF( ( EIGHT*SMLNUM )*ABS( BTMP( 1 ) ).GT.ABS( T16( 1, 1 ) ) .OR.
 
347
     $    ( EIGHT*SMLNUM )*ABS( BTMP( 2 ) ).GT.ABS( T16( 2, 2 ) ) .OR.
 
348
     $    ( EIGHT*SMLNUM )*ABS( BTMP( 3 ) ).GT.ABS( T16( 3, 3 ) ) .OR.
 
349
     $    ( EIGHT*SMLNUM )*ABS( BTMP( 4 ) ).GT.ABS( T16( 4, 4 ) ) ) THEN
 
350
         SCALE = ( ONE / EIGHT ) / MAX( ABS( BTMP( 1 ) ),
 
351
     $           ABS( BTMP( 2 ) ), ABS( BTMP( 3 ) ), ABS( BTMP( 4 ) ) )
 
352
         BTMP( 1 ) = BTMP( 1 )*SCALE
 
353
         BTMP( 2 ) = BTMP( 2 )*SCALE
 
354
         BTMP( 3 ) = BTMP( 3 )*SCALE
 
355
         BTMP( 4 ) = BTMP( 4 )*SCALE
 
356
      END IF
 
357
      DO 120 I = 1, 4
 
358
         K = 5 - I
 
359
         TEMP = ONE / T16( K, K )
 
360
         TMP( K ) = BTMP( K )*TEMP
 
361
         DO 110 J = K + 1, 4
 
362
            TMP( K ) = TMP( K ) - ( TEMP*T16( K, J ) )*TMP( J )
 
363
  110    CONTINUE
 
364
  120 CONTINUE
 
365
      DO 130 I = 1, 3
 
366
         IF( JPIV( 4-I ).NE.4-I ) THEN
 
367
            TEMP = TMP( 4-I )
 
368
            TMP( 4-I ) = TMP( JPIV( 4-I ) )
 
369
            TMP( JPIV( 4-I ) ) = TEMP
 
370
         END IF
 
371
  130 CONTINUE
 
372
      X( 1, 1 ) = TMP( 1 )
 
373
      X( 2, 1 ) = TMP( 2 )
 
374
      X( 1, 2 ) = TMP( 3 )
 
375
      X( 2, 2 ) = TMP( 4 )
 
376
      XNORM = MAX( ABS( TMP( 1 ) )+ABS( TMP( 3 ) ),
 
377
     $        ABS( TMP( 2 ) )+ABS( TMP( 4 ) ) )
 
378
      RETURN
 
379
*
 
380
*     End of DLASY2
 
381
*
 
382
      END