← Back to branch summary

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

« back to all changes in this revision

Viewing changes to routines/slicot/mb01ry.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/mb01ry.f
 
1
      SUBROUTINE MB01RY( SIDE, UPLO, TRANS, M, ALPHA, BETA, R, LDR, H,
 
2
     $                   LDH, B, LDB, DWORK, INFO )
 
3
C
 
4
C     RELEASE 4.0, WGS COPYRIGHT 1999.
 
5
C
 
6
C     PURPOSE
 
7
C
 
8
C     To compute either the upper or lower triangular part of one of the
 
9
C     matrix formulas
 
10
C        _
 
11
C        R = alpha*R + beta*op( H )*B,                               (1)
 
12
C        _
 
13
C        R = alpha*R + beta*B*op( H ),                               (2)
 
14
C                                                    _     
 
15
C     where alpha and beta are scalars, H, B, R, and R are m-by-m
 
16
C     matrices, H is an upper Hessenberg matrix, and op( H ) is one of
 
17
C     
 
18
C        op( H ) = H   or   op( H ) = H',  the transpose of H.
 
19
C     
 
20
C     The result is overwritten on R. 
 
21
C
 
22
C     ARGUMENTS
 
23
C
 
24
C     Mode Parameters
 
25
C
 
26
C     SIDE    CHARACTER*1
 
27
C             Specifies whether the Hessenberg matrix H appears on the 
 
28
C             left or right in the matrix product as follows:
 
29
C                     _                           
 
30
C             = 'L':  R = alpha*R + beta*op( H )*B;
 
31
C                     _                           
 
32
C             = 'R':  R = alpha*R + beta*B*op( H ).
 
33
C            
 
34
C     UPLO    CHARACTER*1                               _
 
35
C             Specifies which triangles of the matrices R and R are
 
36
C             computed and given, respectively, as follows:
 
37
C             = 'U':  the upper triangular part;
 
38
C             = 'L':  the lower triangular part.
 
39
C            
 
40
C     TRANS   CHARACTER*1
 
41
C             Specifies the form of op( H ) to be used in the matrix
 
42
C             multiplication as follows:
 
43
C             = 'N':  op( H ) = H;
 
44
C             = 'T':  op( H ) = H';
 
45
C             = 'C':  op( H ) = H'.
 
46
C            
 
47
C     Input/Output Parameters
 
48
C
 
49
C     M       (input) INTEGER           _
 
50
C             The order of the matrices R, R, H and B.  M >= 0.
 
51
C
 
52
C     ALPHA   (input) DOUBLE PRECISION
 
53
C             The scalar alpha. When alpha is zero then R need not be  
 
54
C             set before entry.
 
55
C            
 
56
C     BETA    (input) DOUBLE PRECISION
 
57
C             The scalar beta. When beta is zero then H and B are not
 
58
C             referenced.
 
59
C            
 
60
C     R       (input/output) DOUBLE PRECISION array, dimension (LDR,M)
 
61
C             On entry with UPLO = 'U', the leading M-by-M upper
 
62
C             triangular part of this array must contain the upper
 
63
C             triangular part of the matrix R; the strictly lower
 
64
C             triangular part of the array is not referenced.
 
65
C             On entry with UPLO = 'L', the leading M-by-M lower
 
66
C             triangular part of this array must contain the lower
 
67
C             triangular part of the matrix R; the strictly upper
 
68
C             triangular part of the array is not referenced.
 
69
C             On exit, the leading M-by-M upper triangular part (if 
 
70
C             UPLO = 'U'), or lower triangular part (if UPLO = 'L') of
 
71
C             this array contains the corresponding triangular part of
 
72
C                                 _
 
73
C             the computed matrix R. 
 
74
C
 
75
C     LDR     INTEGER
 
76
C             The leading dimension of array R.  LDR >= MAX(1,M).
 
77
C
 
78
C     H       (input) DOUBLE PRECISION array, dimension (LDH,M)
 
79
C             On entry, the leading M-by-M upper Hessenberg part of
 
80
C             this array must contain the upper Hessenberg part of the
 
81
C             matrix H.
 
82
C             The elements below the subdiagonal are not referenced,
 
83
C             except possibly for those in the first column, which
 
84
C             could be overwritten, but are restored on exit.
 
85
C
 
86
C     LDH     INTEGER
 
87
C             The leading dimension of array H.  LDH >= MAX(1,M).
 
88
C
 
89
C     B       (input) DOUBLE PRECISION array, dimension (LDB,M)
 
90
C             On entry, the leading M-by-M part of this array must
 
91
C             contain the matrix B.
 
92
C
 
93
C     LDB     INTEGER
 
94
C             The leading dimension of array B.  LDB >= MAX(1,M).
 
95
C
 
96
C     Workspace
 
97
C
 
98
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
 
99
C             LDWORK >= M, if  beta <> 0 and SIDE = 'L';
 
100
C             LDWORK >= 0, if  beta =  0 or  SIDE = 'R'.
 
101
C             This array is not referenced when beta = 0 or SIDE = 'R'.
 
102
C
 
103
C     Error Indicator
 
104
C
 
105
C     INFO    INTEGER
 
106
C             = 0:  successful exit;
 
107
C             < 0:  if INFO = -i, the i-th argument had an illegal
 
108
C                   value.
 
109
C
 
110
C     METHOD
 
111
C
 
112
C     The matrix expression is efficiently evaluated taking the
 
113
C     Hessenberg/triangular structure into account. BLAS 2 operations
 
114
C     are used. A block algorithm can be constructed; it can use BLAS 3
 
115
C     GEMM operations for most computations, and calls of this BLAS 2
 
116
C     algorithm for computing the triangles.
 
117
C
 
118
C     FURTHER COMMENTS
 
119
C
 
120
C     The main application of this routine is when the result should
 
121
C     be a symmetric matrix, e.g., when B = X*op( H )', for (1), or
 
122
C     B = op( H )'*X, for (2), where B is already available and X = X'.
 
123
C     
 
124
C     CONTRIBUTORS
 
125
C
 
126
C     V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1999.
 
127
C
 
128
C     REVISIONS
 
129
C
 
130
C     -
 
131
C
 
132
C     KEYWORDS
 
133
C
 
134
C     Elementary matrix operations, matrix algebra, matrix operations.
 
135
C
 
136
C     ******************************************************************
 
137
C
 
138
C     .. Parameters ..
 
139
      DOUBLE PRECISION  ZERO, ONE
 
140
      PARAMETER         ( ZERO = 0.0D0, ONE = 1.0D0 )
 
141
C     .. Scalar Arguments ..
 
142
      CHARACTER         SIDE, TRANS, UPLO
 
143
      INTEGER           INFO, LDB, LDH, LDR, M
 
144
      DOUBLE PRECISION  ALPHA, BETA
 
145
C     .. Array Arguments ..
 
146
      DOUBLE PRECISION  B(LDB,*), DWORK(*), H(LDH,*), R(LDR,*)
 
147
C     .. Local Scalars ..
 
148
      LOGICAL           LSIDE, LTRANS, LUPLO
 
149
      INTEGER           I, J
 
150
C     .. External Functions ..
 
151
      LOGICAL           LSAME
 
152
      DOUBLE PRECISION  DDOT
 
153
      EXTERNAL          DDOT, LSAME
 
154
C     .. External Subroutines ..
 
155
      EXTERNAL          DCOPY, DGEMV, DLASCL, DLASET, DSCAL, DSWAP,
 
156
     $                  DTRMV, XERBLA
 
157
C     .. Intrinsic Functions ..
 
158
      INTRINSIC         MAX, MIN
 
159
C     .. Executable Statements ..
 
160
C
 
161
C     Test the input scalar arguments.
 
162
C
 
163
      INFO   = 0
 
164
      LSIDE  = LSAME( SIDE,  'L' )
 
165
      LUPLO  = LSAME( UPLO,  'U' )
 
166
      LTRANS = LSAME( TRANS, 'T' ) .OR. LSAME( TRANS, 'C' )
 
167
C
 
168
      IF(      ( .NOT.LSIDE  ).AND.( .NOT.LSAME( SIDE,  'R' ) ) )THEN
 
169
         INFO = -1
 
170
      ELSE IF( ( .NOT.LUPLO  ).AND.( .NOT.LSAME( UPLO,  'L' ) ) )THEN
 
171
         INFO = -2
 
172
      ELSE IF( ( .NOT.LTRANS ).AND.( .NOT.LSAME( TRANS, 'N' ) ) )THEN
 
173
         INFO = -3
 
174
      ELSE IF( M.LT.0 ) THEN
 
175
         INFO = -4
 
176
      ELSE IF( LDR.LT.MAX( 1, M ) ) THEN
 
177
         INFO = -8
 
178
      ELSE IF( LDH.LT.MAX( 1, M ) ) THEN
 
179
         INFO = -10    
 
180
      ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
 
181
         INFO = -12
 
182
      END IF
 
183
C
 
184
      IF ( INFO.NE.0 ) THEN
 
185
C
 
186
C        Error return.
 
187
C
 
188
         CALL XERBLA( 'MB01RY', -INFO )
 
189
         RETURN
 
190
      END IF
 
191
C
 
192
C     Quick return if possible.
 
193
C
 
194
      IF ( M.EQ.0 ) 
 
195
     $   RETURN
 
196
C
 
197
      IF ( BETA.EQ.ZERO ) THEN
 
198
         IF ( ALPHA.EQ.ZERO ) THEN
 
199
C
 
200
C           Special case when both alpha = 0 and beta = 0.
 
201
C
 
202
            CALL DLASET( UPLO, M, M, ZERO, ZERO, R, LDR )
 
203
         ELSE
 
204
C
 
205
C           Special case beta = 0.
 
206
C
 
207
            IF ( ALPHA.NE.ONE ) 
 
208
     $         CALL DLASCL( UPLO, 0, 0, ONE, ALPHA, M, M, R, LDR, INFO )
 
209
         END IF
 
210
         RETURN
 
211
      END IF
 
212
C     
 
213
C     General case: beta <> 0.
 
214
C     Compute the required triangle of (1) or (2) using BLAS 2 
 
215
C     operations. 
 
216
C
 
217
      IF( LSIDE ) THEN
 
218
C
 
219
C        To avoid repeated references to the subdiagonal elements of H, 
 
220
C        these are swapped with the corresponding elements of H in the
 
221
C        first column, and are finally restored.
 
222
C 
 
223
         IF( M.GT.2 ) 
 
224
     $      CALL DSWAP( M-2, H( 3, 2 ), LDH+1, H( 3, 1 ), 1 )
 
225
C
 
226
         IF( LUPLO ) THEN
 
227
            IF ( LTRANS ) THEN
 
228
C
 
229
               DO 20 J = 1, M
 
230
C              
 
231
C                 Multiply the transposed upper triangle of the leading
 
232
C                 j-by-j submatrix of H by the leading part of the j-th
 
233
C                 column of B.
 
234
C
 
235
                  CALL DCOPY( J, B( 1, J ), 1, DWORK, 1 )
 
236
                  CALL DTRMV( 'Upper', TRANS, 'Non-unit', J, H, LDH,
 
237
     $                        DWORK, 1 )
 
238
C              
 
239
C                 Add the contribution of the subdiagonal of H to
 
240
C                 the j-th column of the product.
 
241
C              
 
242
                  DO 10 I = 1, MIN( J, M - 1 )
 
243
                     R( I, J ) = ALPHA*R( I, J ) + BETA*( DWORK( I ) +
 
244
     $                           H( I+1, 1 )*B( I+1, J ) )
 
245
   10             CONTINUE
 
246
C
 
247
   20          CONTINUE
 
248
C
 
249
               R( M, M ) = ALPHA*R( M, M ) + BETA*DWORK( M )
 
250
C     
 
251
            ELSE
 
252
C     
 
253
               DO 40 J = 1, M
 
254
C              
 
255
C                 Multiply the upper triangle of the leading j-by-j
 
256
C                 submatrix of H by the leading part of the j-th column
 
257
C                 of B.
 
258
C
 
259
                  CALL DCOPY( J, B( 1, J ), 1, DWORK, 1 )
 
260
                  CALL DTRMV( 'Upper', TRANS, 'Non-unit', J, H, LDH,
 
261
     $                        DWORK, 1 )
 
262
                  IF( J.LT.M ) THEN
 
263
C
 
264
C                    Multiply the remaining right part of the leading
 
265
C                    j-by-M submatrix of H by the trailing part of the
 
266
C                    j-th column of B.
 
267
C                    
 
268
                     CALL DGEMV( TRANS, J, M-J, BETA, H( 1, J+1 ), LDH,
 
269
     $                           B( J+1, J ), 1, ALPHA, R( 1, J ), 1 )
 
270
                  ELSE
 
271
                     CALL DSCAL( M, ALPHA, R( 1, M ), 1 ) 
 
272
                  END IF
 
273
C              
 
274
C                 Add the contribution of the subdiagonal of H to
 
275
C                 the j-th column of the product.
 
276
C
 
277
                  R( 1, J ) = R( 1, J ) + BETA*DWORK( 1 )
 
278
C              
 
279
                  DO 30 I = 2, J
 
280
                     R( I, J ) = R( I, J ) + BETA*( DWORK( I ) +
 
281
     $                           H( I, 1 )*B( I-1, J ) )
 
282
   30             CONTINUE
 
283
C
 
284
   40          CONTINUE
 
285
C
 
286
            END IF
 
287
C
 
288
         ELSE
 
289
C
 
290
            IF ( LTRANS ) THEN
 
291
C
 
292
               DO 60 J = M, 1, -1
 
293
C              
 
294
C                 Multiply the transposed upper triangle of the trailing
 
295
C                 (M-j+1)-by-(M-j+1) submatrix of H by the trailing part
 
296
C                 of the j-th column of B.
 
297
C
 
298
                  CALL DCOPY( M-J+1, B( J, J ), 1, DWORK( J ), 1 )
 
299
                  CALL DTRMV( 'Upper', TRANS, 'Non-unit', M-J+1,
 
300
     $                        H( J, J ), LDH, DWORK( J ), 1 )
 
301
                  IF( J.GT.1 ) THEN
 
302
C
 
303
C                    Multiply the remaining left part of the trailing 
 
304
C                    (M-j+1)-by-(j-1) submatrix of H' by the leading
 
305
C                    part of the j-th column of B.
 
306
C                    
 
307
                     CALL DGEMV( TRANS, J-1, M-J+1, BETA, H( 1, J ),
 
308
     $                           LDH, B( 1, J ), 1, ALPHA, R( J, J ), 
 
309
     $                           1 )
 
310
                  ELSE
 
311
                     CALL DSCAL( M, ALPHA, R( 1, 1 ), 1 ) 
 
312
                  END IF
 
313
C              
 
314
C                 Add the contribution of the subdiagonal of H to
 
315
C                 the j-th column of the product.
 
316
C
 
317
                  DO 50 I = J, M - 1
 
318
                     R( I, J ) = R( I, J ) + BETA*( DWORK( I ) +
 
319
     $                           H( I+1, 1 )*B( I+1, J ) )
 
320
   50             CONTINUE
 
321
C
 
322
                  R( M, J ) = R( M, J ) + BETA*DWORK( M )
 
323
   60          CONTINUE
 
324
C              
 
325
            ELSE
 
326
C
 
327
               DO 80 J = M, 1, -1
 
328
C              
 
329
C                 Multiply the upper triangle of the trailing
 
330
C                 (M-j+1)-by-(M-j+1) submatrix of H by the trailing
 
331
C                 part of the j-th column of B.
 
332
C
 
333
                  CALL DCOPY( M-J+1, B( J, J ), 1, DWORK( J ), 1 )
 
334
                  CALL DTRMV( 'Upper', TRANS, 'Non-unit', M-J+1,
 
335
     $                        H( J, J ), LDH, DWORK( J ), 1 )
 
336
C              
 
337
C                 Add the contribution of the subdiagonal of H to
 
338
C                 the j-th column of the product.
 
339
C                 
 
340
                  DO 70 I = MAX( J, 2 ), M
 
341
                     R( I, J ) = ALPHA*R( I, J ) + BETA*( DWORK( I )
 
342
     $                               + H( I, 1 )*B( I-1, J ) )
 
343
   70             CONTINUE
 
344
C
 
345
   80          CONTINUE
 
346
C
 
347
               R( 1, 1 ) = ALPHA*R( 1, 1 ) + BETA*DWORK( 1 )
 
348
C
 
349
            END IF
 
350
         END IF
 
351
C     
 
352
         IF( M.GT.2 ) 
 
353
     $      CALL DSWAP( M-2, H( 3, 2 ), LDH+1, H( 3, 1 ), 1 )
 
354
C
 
355
      ELSE
 
356
C     
 
357
C        Row-wise calculations are used for H, if SIDE = 'R' and 
 
358
C        TRANS = 'T'.
 
359
C           
 
360
         IF( LUPLO ) THEN
 
361
            IF( LTRANS ) THEN
 
362
               R( 1, 1 ) = ALPHA*R( 1, 1 ) + 
 
363
     $                     BETA*DDOT( M, B, LDB, H, LDH )
 
364
C
 
365
               DO 90 J = 2, M
 
366
                  CALL DGEMV( 'NoTranspose', J, M-J+2, BETA,
 
367
     $                        B( 1, J-1 ), LDB, H( J, J-1 ), LDH,
 
368
     $                        ALPHA, R( 1, J ), 1 )
 
369
   90          CONTINUE
 
370
C
 
371
            ELSE
 
372
C
 
373
               DO 100 J = 1, M - 1
 
374
                  CALL DGEMV( 'NoTranspose', J, J+1, BETA, B, LDB,
 
375
     $                        H( 1, J ), 1, ALPHA, R( 1, J ), 1 )
 
376
  100          CONTINUE
 
377
C
 
378
               CALL DGEMV( 'NoTranspose', M, M, BETA, B, LDB,
 
379
     $                     H( 1, M ), 1, ALPHA, R( 1, M ), 1 )
 
380
C     
 
381
            END IF
 
382
C
 
383
         ELSE
 
384
C
 
385
            IF( LTRANS ) THEN
 
386
C
 
387
               CALL DGEMV( 'NoTranspose', M, M, BETA, B, LDB, H, LDH,
 
388
     $                     ALPHA, R( 1, 1 ), 1 )
 
389
C
 
390
               DO 110 J = 2, M
 
391
                  CALL DGEMV( 'NoTranspose', M-J+1, M-J+2, BETA,
 
392
     $                        B( J, J-1 ), LDB, H( J, J-1 ), LDH, ALPHA,
 
393
     $                        R( J, J ), 1 )
 
394
  110          CONTINUE
 
395
C     
 
396
            ELSE
 
397
C     
 
398
               DO 120 J = 1, M - 1
 
399
                  CALL DGEMV( 'NoTranspose', M-J+1, J+1, BETA,
 
400
     $                        B( J, 1 ), LDB, H( 1, J ), 1, ALPHA,
 
401
     $                        R( J, J ), 1 )
 
402
  120          CONTINUE
 
403
C
 
404
               R( M, M ) = ALPHA*R( M, M ) + 
 
405
     $                     BETA*DDOT( M, B( M, 1 ), LDB, H( 1, M ), 1 )
 
406
C     
 
407
            END IF
 
408
         END IF
 
409
      END IF
 
410
C
 
411
      RETURN
 
412
C *** Last line of MB01RY ***
 
413
      END