← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/slicot/mb04oy.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/mb04oy.f
 
1
      SUBROUTINE MB04OY( M, N, V, TAU, A, LDA, B, LDB, DWORK )
 
2
C
 
3
C     RELEASE 4.0, WGS COPYRIGHT 1999.
 
4
C
 
5
C     PURPOSE
 
6
C
 
7
C     To apply a real elementary reflector H to a real (m+1)-by-n
 
8
C     matrix C = [ A ], from the left, where A has one row. H is
 
9
C                [ B ]
 
10
C     represented in the form
 
11
C                                        ( 1 )
 
12
C           H = I - tau * u *u',    u  = (   ),
 
13
C                                        ( v )
 
14
C     where tau is a real scalar and v is a real m-vector.
 
15
C
 
16
C     If tau = 0, then H is taken to be the unit matrix.
 
17
C
 
18
C     In-line code is used if H has order < 11.
 
19
C
 
20
C     ARGUMENTS
 
21
C
 
22
C     Input/Output Parameters
 
23
C
 
24
C     M       (input) INTEGER
 
25
C             The number of rows of the matrix B.  M >= 0.
 
26
C
 
27
C     N       (input) INTEGER
 
28
C             The number of columns of the matrices A and B.  N >= 0.
 
29
C
 
30
C     V       (input) DOUBLE PRECISION array, dimension (M)
 
31
C             The vector v in the representation of H.
 
32
C
 
33
C     TAU     (input) DOUBLE PRECISION
 
34
C             The scalar factor of the elementary reflector H.
 
35
C
 
36
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 
37
C             On entry, the leading 1-by-N part of this array must
 
38
C             contain the matrix A.
 
39
C             On exit, the leading 1-by-N part of this array contains
 
40
C             the updated matrix A (the first row of H * C).
 
41
C
 
42
C     LDA     INTEGER
 
43
C             The leading dimension of array A.  LDA >= 1.
 
44
C
 
45
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
 
46
C             On entry, the leading M-by-N part of this array must
 
47
C             contain the matrix B. 
 
48
C             On exit, the leading M-by-N part of this array contains
 
49
C             the updated matrix B (the last m rows of H * C).
 
50
C
 
51
C     LDB     INTEGER
 
52
C             The leading dimension of array B.  LDB >= MAX(1,M).
 
53
C
 
54
C     Workspace
 
55
C
 
56
C     DWORK   DOUBLE PRECISION array, dimension (N)
 
57
C             DWORK is not referenced if H has order less than 11.
 
58
C
 
59
C     METHOD
 
60
C
 
61
C     The routine applies the elementary reflector H, taking the special
 
62
C     structure of C into account.
 
63
C  
 
64
C     NUMERICAL ASPECTS
 
65
C
 
66
C     The algorithm is backward stable.
 
67
C
 
68
C     CONTRIBUTORS
 
69
C
 
70
C     V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1997.
 
71
C     Based on LAPACK routines DLARFX and DLATZM.
 
72
C
 
73
C     REVISIONS
 
74
C
 
75
C     Dec. 1997.
 
76
C
 
77
C     KEYWORDS
 
78
C
 
79
C     Elementary matrix operations, elementary reflector, orthogonal
 
80
C     transformation.
 
81
C
 
82
C     ******************************************************************
 
83
C
 
84
C     .. Parameters ..
 
85
      DOUBLE PRECISION  ZERO, ONE
 
86
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
 
87
C     .. Scalar Arguments ..
 
88
      INTEGER           LDA, LDB, M, N
 
89
      DOUBLE PRECISION  TAU
 
90
C     .. Array Arguments ..
 
91
      DOUBLE PRECISION  A( LDA, * ), B( LDB, * ), DWORK( * ), V( * )
 
92
C     .. Local Scalars ..
 
93
      INTEGER           J
 
94
      DOUBLE PRECISION  SUM, T1, T2, T3, T4, T5, T6, T7, T8, T9, V1, V2,
 
95
     $                  V3, V4, V5, V6, V7, V8, V9
 
96
C     .. External Subroutines ..
 
97
      EXTERNAL          DAXPY, DCOPY, DGEMV, DGER
 
98
C
 
99
C     .. Executable Statements ..
 
100
C
 
101
      IF( TAU.EQ.ZERO )
 
102
     $   RETURN
 
103
C
 
104
C     Form  H * C, where H has order m+1.
 
105
C
 
106
      GO TO ( 10, 30, 50, 70, 90, 110, 130, 150,
 
107
     $        170, 190 ) M+1
 
108
C
 
109
C     Code for general M. Compute
 
110
C
 
111
C     w := C'*u,  C := C - tau * u * w'.
 
112
C
 
113
      CALL DCOPY( N, A, LDA, DWORK, 1 )
 
114
      CALL DGEMV( 'Transpose', M, N, ONE, B, LDB, V, 1, ONE, DWORK, 1 )
 
115
      CALL DAXPY( N, -TAU, DWORK, 1, A, LDA )
 
116
      CALL DGER( M, N, -TAU, V, 1, DWORK, 1, B, LDB )
 
117
      GO TO 210
 
118
   10 CONTINUE
 
119
C
 
120
C     Special code for 1 x 1 Householder
 
121
C
 
122
      T1 = ONE - TAU
 
123
      DO 20 J = 1, N
 
124
         A( 1, J ) = T1*A( 1, J )
 
125
   20 CONTINUE
 
126
      GO TO 210
 
127
   30 CONTINUE
 
128
C
 
129
C     Special code for 2 x 2 Householder
 
130
C
 
131
      V1 = V( 1 )
 
132
      T1 = TAU*V1
 
133
      DO 40 J = 1, N
 
134
         SUM = A( 1, J ) + V1*B( 1, J )
 
135
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
136
         B( 1, J ) = B( 1, J ) - SUM*T1
 
137
   40 CONTINUE
 
138
      GO TO 210
 
139
   50 CONTINUE
 
140
C
 
141
C     Special code for 3 x 3 Householder
 
142
C
 
143
      V1 = V( 1 )
 
144
      T1 = TAU*V1
 
145
      V2 = V( 2 )
 
146
      T2 = TAU*V2
 
147
      DO 60 J = 1, N
 
148
         SUM = A( 1, J ) + V1*B( 1, J ) + V2*B( 2, J )
 
149
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
150
         B( 1, J ) = B( 1, J ) - SUM*T1
 
151
         B( 2, J ) = B( 2, J ) - SUM*T2
 
152
   60 CONTINUE
 
153
      GO TO 210
 
154
   70 CONTINUE
 
155
C
 
156
C     Special code for 4 x 4 Householder
 
157
C
 
158
      V1 = V( 1 )
 
159
      T1 = TAU*V1
 
160
      V2 = V( 2 )
 
161
      T2 = TAU*V2
 
162
      V3 = V( 3 )
 
163
      T3 = TAU*V3
 
164
      DO 80 J = 1, N
 
165
         SUM = A( 1, J ) + V1*B( 1, J ) + V2*B( 2, J ) + V3*B( 3, J )
 
166
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
167
         B( 1, J ) = B( 1, J ) - SUM*T1
 
168
         B( 2, J ) = B( 2, J ) - SUM*T2
 
169
         B( 3, J ) = B( 3, J ) - SUM*T3
 
170
   80 CONTINUE
 
171
      GO TO 210
 
172
   90 CONTINUE
 
173
C
 
174
C     Special code for 5 x 5 Householder
 
175
C
 
176
      V1 = V( 1 )
 
177
      T1 = TAU*V1
 
178
      V2 = V( 2 )
 
179
      T2 = TAU*V2
 
180
      V3 = V( 3 )
 
181
      T3 = TAU*V3
 
182
      V4 = V( 4 )
 
183
      T4 = TAU*V4
 
184
      DO 100 J = 1, N
 
185
         SUM = A( 1, J ) +  V1*B( 1, J ) + V2*B( 2, J ) + V3*B( 3, J ) +
 
186
     $                      V4*B( 4, J )
 
187
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
188
         B( 1, J ) = B( 1, J ) - SUM*T1
 
189
         B( 2, J ) = B( 2, J ) - SUM*T2
 
190
         B( 3, J ) = B( 3, J ) - SUM*T3
 
191
         B( 4, J ) = B( 4, J ) - SUM*T4
 
192
  100 CONTINUE
 
193
      GO TO 210
 
194
  110 CONTINUE
 
195
C
 
196
C     Special code for 6 x 6 Householder
 
197
C
 
198
      V1 = V( 1 )
 
199
      T1 = TAU*V1
 
200
      V2 = V( 2 )
 
201
      T2 = TAU*V2
 
202
      V3 = V( 3 )
 
203
      T3 = TAU*V3
 
204
      V4 = V( 4 )
 
205
      T4 = TAU*V4
 
206
      V5 = V( 5 )
 
207
      T5 = TAU*V5
 
208
      DO 120 J = 1, N
 
209
         SUM = A( 1, J ) +  V1*B( 1, J ) + V2*B( 2, J ) + V3*B( 3, J ) +
 
210
     $                      V4*B( 4, J ) + V5*B( 5, J )
 
211
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
212
         B( 1, J ) = B( 1, J ) - SUM*T1
 
213
         B( 2, J ) = B( 2, J ) - SUM*T2
 
214
         B( 3, J ) = B( 3, J ) - SUM*T3
 
215
         B( 4, J ) = B( 4, J ) - SUM*T4
 
216
         B( 5, J ) = B( 5, J ) - SUM*T5
 
217
  120 CONTINUE
 
218
      GO TO 210
 
219
  130 CONTINUE
 
220
C
 
221
C     Special code for 7 x 7 Householder
 
222
C
 
223
      V1 = V( 1 )
 
224
      T1 = TAU*V1
 
225
      V2 = V( 2 )
 
226
      T2 = TAU*V2
 
227
      V3 = V( 3 )
 
228
      T3 = TAU*V3
 
229
      V4 = V( 4 )
 
230
      T4 = TAU*V4
 
231
      V5 = V( 5 )
 
232
      T5 = TAU*V5
 
233
      V6 = V( 6 )
 
234
      T6 = TAU*V6
 
235
      DO 140 J = 1, N
 
236
         SUM = A( 1, J ) +  V1*B( 1, J ) + V2*B( 2, J ) + V3*B( 3, J ) +
 
237
     $                      V4*B( 4, J ) + V5*B( 5, J ) + V6*B( 6, J )
 
238
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
239
         B( 1, J ) = B( 1, J ) - SUM*T1
 
240
         B( 2, J ) = B( 2, J ) - SUM*T2
 
241
         B( 3, J ) = B( 3, J ) - SUM*T3
 
242
         B( 4, J ) = B( 4, J ) - SUM*T4
 
243
         B( 5, J ) = B( 5, J ) - SUM*T5
 
244
         B( 6, J ) = B( 6, J ) - SUM*T6
 
245
  140 CONTINUE
 
246
      GO TO 210
 
247
  150 CONTINUE
 
248
C
 
249
C     Special code for 8 x 8 Householder
 
250
C
 
251
      V1 = V( 1 )
 
252
      T1 = TAU*V1
 
253
      V2 = V( 2 )
 
254
      T2 = TAU*V2
 
255
      V3 = V( 3 )
 
256
      T3 = TAU*V3
 
257
      V4 = V( 4 )
 
258
      T4 = TAU*V4
 
259
      V5 = V( 5 )
 
260
      T5 = TAU*V5
 
261
      V6 = V( 6 )
 
262
      T6 = TAU*V6
 
263
      V7 = V( 7 )
 
264
      T7 = TAU*V7
 
265
      DO 160 J = 1, N
 
266
         SUM = A( 1, J ) +  V1*B( 1, J ) + V2*B( 2, J ) + V3*B( 3, J ) +
 
267
     $                      V4*B( 4, J ) + V5*B( 5, J ) + V6*B( 6, J ) +
 
268
     $                      V7*B( 7, J )
 
269
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
270
         B( 1, J ) = B( 1, J ) - SUM*T1
 
271
         B( 2, J ) = B( 2, J ) - SUM*T2
 
272
         B( 3, J ) = B( 3, J ) - SUM*T3
 
273
         B( 4, J ) = B( 4, J ) - SUM*T4
 
274
         B( 5, J ) = B( 5, J ) - SUM*T5
 
275
         B( 6, J ) = B( 6, J ) - SUM*T6
 
276
         B( 7, J ) = B( 7, J ) - SUM*T7
 
277
  160 CONTINUE
 
278
      GO TO 210
 
279
  170 CONTINUE
 
280
C
 
281
C     Special code for 9 x 9 Householder
 
282
C
 
283
      V1 = V( 1 )
 
284
      T1 = TAU*V1
 
285
      V2 = V( 2 )
 
286
      T2 = TAU*V2
 
287
      V3 = V( 3 )
 
288
      T3 = TAU*V3
 
289
      V4 = V( 4 )
 
290
      T4 = TAU*V4
 
291
      V5 = V( 5 )
 
292
      T5 = TAU*V5
 
293
      V6 = V( 6 )
 
294
      T6 = TAU*V6
 
295
      V7 = V( 7 )
 
296
      T7 = TAU*V7
 
297
      V8 = V( 8 )
 
298
      T8 = TAU*V8
 
299
      DO 180 J = 1, N
 
300
         SUM = A( 1, J ) +  V1*B( 1, J ) + V2*B( 2, J ) + V3*B( 3, J ) +
 
301
     $                      V4*B( 4, J ) + V5*B( 5, J ) + V6*B( 6, J ) +
 
302
     $                      V7*B( 7, J ) + V8*B( 8, J )
 
303
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
304
         B( 1, J ) = B( 1, J ) - SUM*T1
 
305
         B( 2, J ) = B( 2, J ) - SUM*T2
 
306
         B( 3, J ) = B( 3, J ) - SUM*T3
 
307
         B( 4, J ) = B( 4, J ) - SUM*T4
 
308
         B( 5, J ) = B( 5, J ) - SUM*T5
 
309
         B( 6, J ) = B( 6, J ) - SUM*T6
 
310
         B( 7, J ) = B( 7, J ) - SUM*T7
 
311
         B( 8, J ) = B( 8, J ) - SUM*T8
 
312
  180 CONTINUE
 
313
      GO TO 210
 
314
  190 CONTINUE
 
315
C
 
316
C     Special code for 10 x 10 Householder
 
317
C
 
318
      V1 = V( 1 )
 
319
      T1 = TAU*V1
 
320
      V2 = V( 2 )
 
321
      T2 = TAU*V2
 
322
      V3 = V( 3 )
 
323
      T3 = TAU*V3
 
324
      V4 = V( 4 )
 
325
      T4 = TAU*V4
 
326
      V5 = V( 5 )
 
327
      T5 = TAU*V5
 
328
      V6 = V( 6 )
 
329
      T6 = TAU*V6
 
330
      V7 = V( 7 )
 
331
      T7 = TAU*V7
 
332
      V8 = V( 8 )
 
333
      T8 = TAU*V8
 
334
      V9 = V( 9 )
 
335
      T9 = TAU*V9
 
336
      DO 200 J = 1, N
 
337
         SUM = A( 1, J ) +  V1*B( 1, J ) + V2*B( 2, J ) + V3*B( 3, J ) +
 
338
     $                      V4*B( 4, J ) + V5*B( 5, J ) + V6*B( 6, J ) +
 
339
     $                      V7*B( 7, J ) + V8*B( 8, J ) + V9*B( 9, J )
 
340
         A( 1, J ) = A( 1, J ) - SUM*TAU
 
341
         B( 1, J ) = B( 1, J ) - SUM*T1
 
342
         B( 2, J ) = B( 2, J ) - SUM*T2
 
343
         B( 3, J ) = B( 3, J ) - SUM*T3
 
344
         B( 4, J ) = B( 4, J ) - SUM*T4
 
345
         B( 5, J ) = B( 5, J ) - SUM*T5
 
346
         B( 6, J ) = B( 6, J ) - SUM*T6
 
347
         B( 7, J ) = B( 7, J ) - SUM*T7
 
348
         B( 8, J ) = B( 8, J ) - SUM*T8
 
349
         B( 9, J ) = B( 9, J ) - SUM*T9
 
350
  200 CONTINUE
 
351
  210 CONTINUE
 
352
      RETURN
 
353
C *** Last line of MB04OY ***
 
354
      END