← Back to branch summary

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

« back to all changes in this revision

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