← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/slicot/sb04rx.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/sb04rx.f
 
1
      SUBROUTINE SB04RX( RC, UL, M, A, LDA, LAMBD1, LAMBD2, LAMBD3,
 
2
     $                   LAMBD4, D, TOL, IWORK, DWORK, LDDWOR, INFO )
 
3
C
 
4
C     RELEASE 4.0, WGS COPYRIGHT 2000.
 
5
C
 
6
C     PURPOSE
 
7
C
 
8
C     To solve a system of equations in quasi-Hessenberg form
 
9
C     (Hessenberg form plus two consecutive offdiagonals) with two 
 
10
C     right-hand sides.
 
11
C
 
12
C     ARGUMENTS
 
13
C
 
14
C     Mode Parameters
 
15
C
 
16
C     RC      CHARACTER*1
 
17
C             Indicates processing by columns or rows, as follows:
 
18
C             = 'R':  Row transformations are applied;
 
19
C             = 'C':  Column transformations are applied.
 
20
C
 
21
C     UL      CHARACTER*1
 
22
C             Indicates whether A is upper or lower Hessenberg matrix,
 
23
C             as follows:
 
24
C             = 'U':  A is upper Hessenberg;
 
25
C             = 'L':  A is lower Hessenberg.
 
26
C
 
27
C     Input/Output Parameters
 
28
C
 
29
C     M       (input) INTEGER
 
30
C             The order of the matrix A.  M >= 0.
 
31
C
 
32
C     A       (input) DOUBLE PRECISION array, dimension (LDA,M)
 
33
C             The leading M-by-M part of this array must contain a
 
34
C             matrix A in Hessenberg form.
 
35
C
 
36
C     LDA     INTEGER
 
37
C             The leading dimension of array A.  LDA >= MAX(1,M).
 
38
C
 
39
C     LAMBD1, (input) DOUBLE PRECISION
 
40
C     LAMBD2, These variables must contain the 2-by-2 block to be 
 
41
C     LAMBD3, multiplied to the elements of A.
 
42
C     LAMBD4
 
43
C
 
44
C     D       (input/output) DOUBLE PRECISION array, dimension (2*M)
 
45
C             On entry, this array must contain the two right-hand
 
46
C             side vectors of the quasi-Hessenberg system, stored 
 
47
C             row-wise.
 
48
C             On exit, if INFO = 0, this array contains the two solution
 
49
C             vectors of the quasi-Hessenberg system, stored row-wise.
 
50
C
 
51
C     Tolerances
 
52
C
 
53
C     TOL     DOUBLE PRECISION
 
54
C             The tolerance to be used to test for near singularity of
 
55
C             the triangular factor R of the quasi-Hessenberg matrix.
 
56
C             A matrix whose estimated condition number is less
 
57
C             than 1/TOL is considered to be nonsingular.
 
58
C
 
59
C     Workspace
 
60
C
 
61
C     IWORK   INTEGER array, dimension (2*M)
 
62
C
 
63
C     DWORK   DOUBLE PRECISION array, dimension (LDDWOR,2*M+3)
 
64
C             The leading 2*M-by-2*M part of this array is used for
 
65
C             computing the triangular factor of the QR decomposition 
 
66
C             of the quasi-Hessenberg matrix. The remaining 6*M elements
 
67
C             are used as workspace for the computation of the
 
68
C             reciprocal condition estimate.
 
69
C
 
70
C     LDDWOR  INTEGER
 
71
C             The leading dimension of array DWORK.
 
72
C             LDDWOR >= MAX(1,2*M).
 
73
C
 
74
C     Error Indicator
 
75
C
 
76
C     INFO    INTEGER
 
77
C             = 0:  successful exit;
 
78
C             = 1:  if the quasi-Hessenberg matrix is (numerically)
 
79
C                   singular. That is, its estimated reciprocal
 
80
C                   condition number is less than or equal to TOL.
 
81
C
 
82
C     NUMERICAL ASPECTS
 
83
C
 
84
C     None.
 
85
C
 
86
C     CONTRIBUTORS
 
87
C
 
88
C     D. Sima, University of Bucharest, May 2000.
 
89
C
 
90
C     REVISIONS
 
91
C
 
92
C     -
 
93
C
 
94
C     Note that RC, UL, M, LDA, and LDDWOR must be such that the value
 
95
C     of the LOGICAL variable OK in the following statement is true.
 
96
C
 
97
C      OK = ( ( UL.EQ.'U' ) .OR. ( UL.EQ.'u' ) .OR.
 
98
C             ( UL.EQ.'L' ) .OR. ( UL.EQ.'l' ) )
 
99
C           .AND.
 
100
C           ( ( RC.EQ.'R' ) .OR. ( RC.EQ.'r' ) .OR.
 
101
C             ( RC.EQ.'C' ) .OR. ( RC.EQ.'c' ) )
 
102
C           .AND.
 
103
C           ( M.GE.0 )
 
104
C           .AND.
 
105
C           ( LDA.GE.MAX( 1, M ) )
 
106
C           .AND.
 
107
C           ( LDDWOR.GE.MAX( 1, 2*M ) )
 
108
C
 
109
C     These conditions are not checked by the routine.
 
110
C
 
111
C     KEYWORDS
 
112
C
 
113
C     Hessenberg form, orthogonal transformation, real Schur form,
 
114
C     Sylvester equation.
 
115
C
 
116
C     ******************************************************************
 
117
C
 
118
      DOUBLE PRECISION  ONE, ZERO
 
119
      PARAMETER         ( ONE = 1.0D0, ZERO = 0.0D0 )
 
120
C     .. Scalar Arguments ..
 
121
      CHARACTER         RC, UL
 
122
      INTEGER           INFO, LDA, LDDWOR, M
 
123
      DOUBLE PRECISION  LAMBD1, LAMBD2, LAMBD3, LAMBD4, TOL
 
124
C     .. Array Arguments ..
 
125
      INTEGER           IWORK(*)
 
126
      DOUBLE PRECISION  A(LDA,*), D(*), DWORK(LDDWOR,*)
 
127
C     .. Local Scalars ..
 
128
      CHARACTER         TRANS
 
129
      INTEGER           J, J1, J2, M2, MJ, ML
 
130
      DOUBLE PRECISION  C, R, RCOND, S
 
131
C     .. External Functions ..
 
132
      LOGICAL           LSAME
 
133
      EXTERNAL          LSAME
 
134
C     .. External Subroutines ..
 
135
      EXTERNAL          DCOPY, DLARTG, DLASET, DROT, DSCAL, DTRCON, 
 
136
     $                  DTRSV
 
137
C     .. Intrinsic Functions ..
 
138
      INTRINSIC         MAX, MIN, MOD
 
139
C     .. Executable Statements ..
 
140
C
 
141
      INFO = 0
 
142
C
 
143
C     For speed, no tests on the input scalar arguments are made.
 
144
C     Quick return if possible.
 
145
C
 
146
      IF ( M.EQ.0 )
 
147
     $   RETURN
 
148
C
 
149
      M2 = M*2
 
150
      IF ( LSAME( UL, 'U' ) ) THEN
 
151
C
 
152
         DO 20 J = 1, M
 
153
            J2 = J*2
 
154
            ML = MIN( M, J + 1 )
 
155
            CALL DLASET( 'Full', M2, 2, ZERO, ZERO, DWORK(1,J2-1),
 
156
     $                   LDDWOR )
 
157
            CALL DCOPY( ML, A(1,J), 1, DWORK(1,J2-1), 2 )
 
158
            CALL DSCAL( ML, LAMBD1, DWORK(1,J2-1), 2 )
 
159
            CALL DCOPY( ML, A(1,J), 1, DWORK(2,J2-1), 2 )
 
160
            CALL DSCAL( ML, LAMBD3, DWORK(2,J2-1), 2 )
 
161
            CALL DCOPY( ML, A(1,J), 1, DWORK(1,J2), 2 )
 
162
            CALL DSCAL( ML, LAMBD2, DWORK(1,J2), 2 )
 
163
            CALL DCOPY( ML, A(1,J), 1, DWORK(2,J2), 2 )
 
164
            CALL DSCAL( ML, LAMBD4, DWORK(2,J2), 2 )
 
165
C
 
166
            DWORK(J2-1,J2-1) = DWORK(J2-1,J2-1) + ONE
 
167
            DWORK(J2,J2) = DWORK(J2,J2) + ONE
 
168
   20    CONTINUE
 
169
C
 
170
         IF ( LSAME( RC, 'R' ) ) THEN
 
171
            TRANS = 'N'
 
172
C
 
173
C           A is an upper Hessenberg matrix, row transformations.
 
174
C
 
175
            DO 40 J = 1, M2 - 1
 
176
               MJ = M2 - J
 
177
               IF ( MOD(J,2).EQ.1 .AND. J.LT.M2-2 ) THEN
 
178
                  IF ( DWORK(J+3,J).NE.ZERO ) THEN
 
179
                     CALL DLARTG( DWORK(J+2,J), DWORK(J+3,J), C, S, R )
 
180
                     DWORK(J+2,J) = R
 
181
                     DWORK(J+3,J) = ZERO
 
182
                     CALL DROT( MJ, DWORK(J+2,J+1), LDDWOR,
 
183
     $                          DWORK(J+3,J+1), LDDWOR, C, S )
 
184
                     CALL DROT( 1, D(J+2), 1, D(J+3), 1, C, S )
 
185
                  END IF
 
186
               END IF                
 
187
               IF ( J.LT.M2-1 ) THEN
 
188
                  IF ( DWORK(J+2,J).NE.ZERO ) THEN
 
189
                     CALL DLARTG( DWORK(J+1,J), DWORK(J+2,J), C, S, R )
 
190
                     DWORK(J+1,J) = R
 
191
                     DWORK(J+2,J) = ZERO
 
192
                     CALL DROT( MJ, DWORK(J+1,J+1), LDDWOR,
 
193
     $                          DWORK(J+2,J+1), LDDWOR, C, S )
 
194
                     CALL DROT( 1, D(J+1), 1, D(J+2), 1, C, S )
 
195
                  END IF
 
196
               END IF
 
197
               IF ( DWORK(J+1,J).NE.ZERO ) THEN
 
198
                  CALL DLARTG( DWORK(J,J), DWORK(J+1,J), C, S, R )
 
199
                  DWORK(J,J)   = R
 
200
                  DWORK(J+1,J) = ZERO
 
201
                  CALL DROT( MJ, DWORK(J,J+1), LDDWOR, DWORK(J+1,J+1),
 
202
     $                       LDDWOR, C, S )
 
203
                  CALL DROT( 1, D(J), 1, D(J+1), 1, C, S )
 
204
               END IF
 
205
   40       CONTINUE
 
206
C
 
207
         ELSE
 
208
            TRANS = 'T'
 
209
C
 
210
C           A is an upper Hessenberg matrix, column transformations.
 
211
C
 
212
            DO 60 J = 1, M2 - 1
 
213
               MJ = M2 - J
 
214
               IF ( MOD(J,2).EQ.1 .AND. J.LT.M2-2 ) THEN
 
215
                  IF ( DWORK(MJ+1,MJ-2).NE.ZERO ) THEN
 
216
                     CALL DLARTG( DWORK(MJ+1,MJ-1), DWORK(MJ+1,MJ-2), C,
 
217
     $                            S, R )
 
218
                     DWORK(MJ+1,MJ-1) = R
 
219
                     DWORK(MJ+1,MJ-2) = ZERO
 
220
                     CALL DROT( MJ, DWORK(1,MJ-1), 1, DWORK(1,MJ-2), 1, 
 
221
     $                          C, S )
 
222
                     CALL DROT( 1, D(MJ-1), 1, D(MJ-2), 1, C, S )
 
223
                  END IF
 
224
               END IF
 
225
               IF ( J.LT.M2-1 ) THEN
 
226
                  IF ( DWORK(MJ+1,MJ-1).NE.ZERO ) THEN
 
227
                     CALL DLARTG( DWORK(MJ+1,MJ), DWORK(MJ+1,MJ-1), C,
 
228
     $                            S, R )
 
229
                     DWORK(MJ+1,MJ)   = R
 
230
                     DWORK(MJ+1,MJ-1) = ZERO
 
231
                     CALL DROT( MJ, DWORK(1,MJ), 1, DWORK(1,MJ-1), 1, C,
 
232
     $                          S )
 
233
                     CALL DROT( 1, D(MJ), 1, D(MJ-1), 1, C, S )
 
234
                  END IF
 
235
               END IF
 
236
               IF ( DWORK(MJ+1,MJ).NE.ZERO ) THEN
 
237
                  CALL DLARTG( DWORK(MJ+1,MJ+1), DWORK(MJ+1,MJ), C, S,
 
238
     $                         R )
 
239
                  DWORK(MJ+1,MJ+1) = R
 
240
                  DWORK(MJ+1,MJ)   = ZERO
 
241
                  CALL DROT( MJ, DWORK(1,MJ+1), 1, DWORK(1,MJ), 1, C,
 
242
     $                       S )
 
243
                  CALL DROT( 1, D(MJ+1), 1, D(MJ), 1, C, S )
 
244
               END IF
 
245
   60       CONTINUE
 
246
C
 
247
         END IF
 
248
      ELSE
 
249
C
 
250
         DO 80 J = 1, M
 
251
            J2 = J*2
 
252
            J1 = MAX( J - 1, 1 )
 
253
            ML = MIN( M - J + 2, M )
 
254
            CALL DLASET( 'Full', M2, 2, ZERO, ZERO, DWORK(1,J2-1),
 
255
     $                   LDDWOR )
 
256
            CALL DCOPY( ML, A(J1,J), 1, DWORK(J1*2-1,J2-1), 2 )
 
257
            CALL DSCAL( ML, LAMBD1, DWORK(J1*2-1,J2-1), 2 )
 
258
            CALL DCOPY( ML, A(J1,J), 1, DWORK(J1*2,J2-1), 2 )
 
259
            CALL DSCAL( ML, LAMBD3, DWORK(J1*2,J2-1), 2 )
 
260
            CALL DCOPY( ML, A(J1,J), 1, DWORK(J1*2-1,J2), 2 )
 
261
            CALL DSCAL( ML, LAMBD2, DWORK(J1*2-1,J2), 2 )
 
262
            CALL DCOPY( ML, A(J1,J), 1, DWORK(J1*2,J2), 2 )
 
263
            CALL DSCAL( ML, LAMBD4, DWORK(J1*2,J2), 2 )
 
264
C
 
265
            DWORK(J2-1,J2-1) = DWORK(J2-1,J2-1) + ONE
 
266
            DWORK(J2,J2) = DWORK(J2,J2) + ONE
 
267
   80    CONTINUE
 
268
C
 
269
         IF ( LSAME( RC, 'R' ) ) THEN
 
270
            TRANS = 'N'
 
271
C
 
272
C           A is a lower Hessenberg matrix, row transformations.
 
273
C
 
274
            DO 100 J = 1, M2 - 1
 
275
               MJ = M2 - J
 
276
               IF ( MOD(J,2).EQ.1 .AND. J.LT.M2-2 ) THEN
 
277
                  IF ( DWORK(MJ-2,MJ+1).NE.ZERO ) THEN
 
278
                     CALL DLARTG( DWORK(MJ-1,MJ+1), DWORK(MJ-2,MJ+1), C,
 
279
     $                            S, R )
 
280
                     DWORK(MJ-1,MJ+1) = R
 
281
                     DWORK(MJ-2,MJ+1) = ZERO
 
282
                     CALL DROT( MJ, DWORK(MJ-1,1), LDDWOR, 
 
283
     $                          DWORK(MJ-2,1), LDDWOR, C, S )
 
284
                     CALL DROT( 1, D(MJ-1), 1, D(MJ-2), 1, C, S )
 
285
                  END IF
 
286
               END IF
 
287
               IF ( J.LT.M2-1 ) THEN
 
288
                  IF ( DWORK(MJ-1,MJ+1).NE.ZERO ) THEN
 
289
                     CALL DLARTG( DWORK(MJ,MJ+1), DWORK(MJ-1,MJ+1), C,
 
290
     $                           S, R )
 
291
                     DWORK(MJ,MJ+1)   = R
 
292
                     DWORK(MJ-1,MJ+1) = ZERO
 
293
                     CALL DROT( MJ, DWORK(MJ,1), LDDWOR, DWORK(MJ-1,1),
 
294
     $                          LDDWOR, C, S )
 
295
                     CALL DROT( 1, D(MJ), 1, D(MJ-1), 1, C, S )
 
296
                  END IF
 
297
               END IF
 
298
               IF ( DWORK(MJ,MJ+1).NE.ZERO ) THEN
 
299
                  CALL DLARTG( DWORK(MJ+1,MJ+1), DWORK(MJ,MJ+1), C, S,
 
300
     $                        R )
 
301
                  DWORK(MJ+1,MJ+1) = R
 
302
                  DWORK(MJ,MJ+1)   = ZERO
 
303
                  CALL DROT( MJ, DWORK(MJ+1,1), LDDWOR, DWORK(MJ,1),
 
304
     $                       LDDWOR, C, S)
 
305
                  CALL DROT( 1, D(MJ+1), 1, D(MJ), 1, C, S )
 
306
               END IF
 
307
  100       CONTINUE
 
308
C
 
309
         ELSE
 
310
            TRANS = 'T'
 
311
C
 
312
C           A is a lower Hessenberg matrix, column transformations.
 
313
C
 
314
            DO 120 J = 1, M2 - 1
 
315
               MJ = M2 - J
 
316
               IF ( MOD(J,2).EQ.1 .AND. J.LT.M2-2 ) THEN
 
317
                  IF ( DWORK(J,J+3).NE.ZERO ) THEN
 
318
                     CALL DLARTG( DWORK(J,J+2), DWORK(J,J+3), C, S, R )
 
319
                     DWORK(J,J+2) = R
 
320
                     DWORK(J,J+3) = ZERO
 
321
                     CALL DROT( MJ, DWORK(J+1,J+2), 1, DWORK(J+1,J+3),
 
322
     $                          1, C, S )
 
323
                     CALL DROT( 1, D(J+2), 1, D(J+3), 1, C, S )
 
324
                  END IF
 
325
               END IF
 
326
               IF ( J.LT.M2-1 ) THEN
 
327
                  IF ( DWORK(J,J+2).NE.ZERO ) THEN
 
328
                     CALL DLARTG( DWORK(J,J+1), DWORK(J,J+2), C, S, R )
 
329
                     DWORK(J,J+1) = R
 
330
                     DWORK(J,J+2) = ZERO
 
331
                     CALL DROT( MJ, DWORK(J+1,J+1), 1, DWORK(J+1,J+2),
 
332
     $                          1, C, S )
 
333
                     CALL DROT( 1, D(J+1), 1, D(J+2), 1, C, S )
 
334
                  END IF
 
335
               END IF
 
336
               IF ( DWORK(J,J+1).NE.ZERO ) THEN
 
337
                  CALL DLARTG( DWORK(J,J), DWORK(J,J+1), C, S, R )
 
338
                  DWORK(J,J)   = R
 
339
                  DWORK(J,J+1) = ZERO
 
340
                  CALL DROT( MJ, DWORK(J+1,J), 1, DWORK(J+1,J+1), 1, C,
 
341
     $                       S )
 
342
                  CALL DROT( 1, D(J), 1, D(J+1), 1, C, S )
 
343
               END IF
 
344
  120       CONTINUE
 
345
C
 
346
         END IF
 
347
      END IF
 
348
C
 
349
      CALL DTRCON( '1-norm', UL, 'Non-unit', M2, DWORK, LDDWOR, RCOND,
 
350
     $             DWORK(1,M2+1), IWORK, INFO )
 
351
      IF ( RCOND.LE.TOL ) THEN
 
352
         INFO = 1
 
353
      ELSE
 
354
         CALL DTRSV( UL, TRANS, 'Non-unit', M2, DWORK, LDDWOR, D, 1 )
 
355
      END IF
 
356
C
 
357
      RETURN
 
358
C *** Last line of SB04RX ***
 
359
      END