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

« back to all changes in this revision

Viewing changes to routines/blas/zhemm.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
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
      SUBROUTINE ZHEMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB,
 
2
     $                   BETA, C, LDC )
 
3
*     .. Scalar Arguments ..
 
4
      CHARACTER*1        SIDE, UPLO
 
5
      INTEGER            M, N, LDA, LDB, LDC
 
6
      COMPLEX*16         ALPHA, BETA
 
7
*     .. Array Arguments ..
 
8
      COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
 
9
*     ..
 
10
*
 
11
*  Purpose
 
12
*  =======
 
13
*
 
14
*  ZHEMM  performs one of the matrix-matrix operations
 
15
*
 
16
*     C := alpha*A*B + beta*C,
 
17
*
 
18
*  or
 
19
*
 
20
*     C := alpha*B*A + beta*C,
 
21
*
 
22
*  where alpha and beta are scalars, A is an hermitian matrix and  B and
 
23
*  C are m by n matrices.
 
24
*
 
25
*  Parameters
 
26
*  ==========
 
27
*
 
28
*  SIDE   - CHARACTER*1.
 
29
*           On entry,  SIDE  specifies whether  the  hermitian matrix  A
 
30
*           appears on the  left or right  in the  operation as follows:
 
31
*
 
32
*              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
 
33
*
 
34
*              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
 
35
*
 
36
*           Unchanged on exit.
 
37
*
 
38
*  UPLO   - CHARACTER*1.
 
39
*           On  entry,   UPLO  specifies  whether  the  upper  or  lower
 
40
*           triangular  part  of  the  hermitian  matrix   A  is  to  be
 
41
*           referenced as follows:
 
42
*
 
43
*              UPLO = 'U' or 'u'   Only the upper triangular part of the
 
44
*                                  hermitian matrix is to be referenced.
 
45
*
 
46
*              UPLO = 'L' or 'l'   Only the lower triangular part of the
 
47
*                                  hermitian matrix is to be referenced.
 
48
*
 
49
*           Unchanged on exit.
 
50
*
 
51
*  M      - INTEGER.
 
52
*           On entry,  M  specifies the number of rows of the matrix  C.
 
53
*           M  must be at least zero.
 
54
*           Unchanged on exit.
 
55
*
 
56
*  N      - INTEGER.
 
57
*           On entry, N specifies the number of columns of the matrix C.
 
58
*           N  must be at least zero.
 
59
*           Unchanged on exit.
 
60
*
 
61
*  ALPHA  - COMPLEX*16      .
 
62
*           On entry, ALPHA specifies the scalar alpha.
 
63
*           Unchanged on exit.
 
64
*
 
65
*  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
 
66
*           m  when  SIDE = 'L' or 'l'  and is n  otherwise.
 
67
*           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
 
68
*           the array  A  must contain the  hermitian matrix,  such that
 
69
*           when  UPLO = 'U' or 'u', the leading m by m upper triangular
 
70
*           part of the array  A  must contain the upper triangular part
 
71
*           of the  hermitian matrix and the  strictly  lower triangular
 
72
*           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
 
73
*           the leading  m by m  lower triangular part  of the  array  A
 
74
*           must  contain  the  lower triangular part  of the  hermitian
 
75
*           matrix and the  strictly upper triangular part of  A  is not
 
76
*           referenced.
 
77
*           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
 
78
*           the array  A  must contain the  hermitian matrix,  such that
 
79
*           when  UPLO = 'U' or 'u', the leading n by n upper triangular
 
80
*           part of the array  A  must contain the upper triangular part
 
81
*           of the  hermitian matrix and the  strictly  lower triangular
 
82
*           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
 
83
*           the leading  n by n  lower triangular part  of the  array  A
 
84
*           must  contain  the  lower triangular part  of the  hermitian
 
85
*           matrix and the  strictly upper triangular part of  A  is not
 
86
*           referenced.
 
87
*           Note that the imaginary parts  of the diagonal elements need
 
88
*           not be set, they are assumed to be zero.
 
89
*           Unchanged on exit.
 
90
*
 
91
*  LDA    - INTEGER.
 
92
*           On entry, LDA specifies the first dimension of A as declared
 
93
*           in the  calling (sub) program. When  SIDE = 'L' or 'l'  then
 
94
*           LDA must be at least  max( 1, m ), otherwise  LDA must be at
 
95
*           least max( 1, n ).
 
96
*           Unchanged on exit.
 
97
*
 
98
*  B      - COMPLEX*16       array of DIMENSION ( LDB, n ).
 
99
*           Before entry, the leading  m by n part of the array  B  must
 
100
*           contain the matrix B.
 
101
*           Unchanged on exit.
 
102
*
 
103
*  LDB    - INTEGER.
 
104
*           On entry, LDB specifies the first dimension of B as declared
 
105
*           in  the  calling  (sub)  program.   LDB  must  be  at  least
 
106
*           max( 1, m ).
 
107
*           Unchanged on exit.
 
108
*
 
109
*  BETA   - COMPLEX*16      .
 
110
*           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
 
111
*           supplied as zero then C need not be set on input.
 
112
*           Unchanged on exit.
 
113
*
 
114
*  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
 
115
*           Before entry, the leading  m by n  part of the array  C must
 
116
*           contain the matrix  C,  except when  beta  is zero, in which
 
117
*           case C need not be set on entry.
 
118
*           On exit, the array  C  is overwritten by the  m by n updated
 
119
*           matrix.
 
120
*
 
121
*  LDC    - INTEGER.
 
122
*           On entry, LDC specifies the first dimension of C as declared
 
123
*           in  the  calling  (sub)  program.   LDC  must  be  at  least
 
124
*           max( 1, m ).
 
125
*           Unchanged on exit.
 
126
*
 
127
*
 
128
*  Level 3 Blas routine.
 
129
*
 
130
*  -- Written on 8-February-1989.
 
131
*     Jack Dongarra, Argonne National Laboratory.
 
132
*     Iain Duff, AERE Harwell.
 
133
*     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 
134
*     Sven Hammarling, Numerical Algorithms Group Ltd.
 
135
*
 
136
*
 
137
*     .. External Functions ..
 
138
      LOGICAL            LSAME
 
139
      EXTERNAL           LSAME
 
140
*     .. External Subroutines ..
 
141
      EXTERNAL           XERBLA
 
142
*     .. Intrinsic Functions ..
 
143
      INTRINSIC          DCONJG, MAX, DBLE
 
144
*     .. Local Scalars ..
 
145
      LOGICAL            UPPER
 
146
      INTEGER            I, INFO, J, K, NROWA
 
147
      COMPLEX*16         TEMP1, TEMP2
 
148
*     .. Parameters ..
 
149
      COMPLEX*16         ONE
 
150
      PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
 
151
      COMPLEX*16         ZERO
 
152
      PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 
153
*     ..
 
154
*     .. Executable Statements ..
 
155
*
 
156
*     Set NROWA as the number of rows of A.
 
157
*
 
158
      IF( LSAME( SIDE, 'L' ) )THEN
 
159
         NROWA = M
 
160
      ELSE
 
161
         NROWA = N
 
162
      END IF
 
163
      UPPER = LSAME( UPLO, 'U' )
 
164
*
 
165
*     Test the input parameters.
 
166
*
 
167
      INFO = 0
 
168
      IF(      ( .NOT.LSAME( SIDE, 'L' ) ).AND.
 
169
     $         ( .NOT.LSAME( SIDE, 'R' ) )      )THEN
 
170
         INFO = 1
 
171
      ELSE IF( ( .NOT.UPPER              ).AND.
 
172
     $         ( .NOT.LSAME( UPLO, 'L' ) )      )THEN
 
173
         INFO = 2
 
174
      ELSE IF( M  .LT.0               )THEN
 
175
         INFO = 3
 
176
      ELSE IF( N  .LT.0               )THEN
 
177
         INFO = 4
 
178
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
 
179
         INFO = 7
 
180
      ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
 
181
         INFO = 9
 
182
      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
 
183
         INFO = 12
 
184
      END IF
 
185
      IF( INFO.NE.0 )THEN
 
186
         CALL XERBLA( 'ZHEMM ', INFO )
 
187
         RETURN
 
188
      END IF
 
189
*
 
190
*     Quick return if possible.
 
191
*
 
192
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
 
193
     $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
 
194
     $   RETURN
 
195
*
 
196
*     And when  alpha.eq.zero.
 
197
*
 
198
      IF( ALPHA.EQ.ZERO )THEN
 
199
         IF( BETA.EQ.ZERO )THEN
 
200
            DO 20, J = 1, N
 
201
               DO 10, I = 1, M
 
202
                  C( I, J ) = ZERO
 
203
   10          CONTINUE
 
204
   20       CONTINUE
 
205
         ELSE
 
206
            DO 40, J = 1, N
 
207
               DO 30, I = 1, M
 
208
                  C( I, J ) = BETA*C( I, J )
 
209
   30          CONTINUE
 
210
   40       CONTINUE
 
211
         END IF
 
212
         RETURN
 
213
      END IF
 
214
*
 
215
*     Start the operations.
 
216
*
 
217
      IF( LSAME( SIDE, 'L' ) )THEN
 
218
*
 
219
*        Form  C := alpha*A*B + beta*C.
 
220
*
 
221
         IF( UPPER )THEN
 
222
            DO 70, J = 1, N
 
223
               DO 60, I = 1, M
 
224
                  TEMP1 = ALPHA*B( I, J )
 
225
                  TEMP2 = ZERO
 
226
                  DO 50, K = 1, I - 1
 
227
                     C( K, J ) = C( K, J ) + TEMP1*A( K, I )
 
228
                     TEMP2     = TEMP2     +
 
229
     $                           B( K, J )*DCONJG( A( K, I ) )
 
230
   50             CONTINUE
 
231
                  IF( BETA.EQ.ZERO )THEN
 
232
                     C( I, J ) = TEMP1*DBLE( A( I, I ) ) +
 
233
     $                           ALPHA*TEMP2
 
234
                  ELSE
 
235
                     C( I, J ) = BETA *C( I, J )         +
 
236
     $                           TEMP1*DBLE( A( I, I ) ) +
 
237
     $                           ALPHA*TEMP2
 
238
                  END IF
 
239
   60          CONTINUE
 
240
   70       CONTINUE
 
241
         ELSE
 
242
            DO 100, J = 1, N
 
243
               DO 90, I = M, 1, -1
 
244
                  TEMP1 = ALPHA*B( I, J )
 
245
                  TEMP2 = ZERO
 
246
                  DO 80, K = I + 1, M
 
247
                     C( K, J ) = C( K, J ) + TEMP1*A( K, I )
 
248
                     TEMP2     = TEMP2     +
 
249
     $                           B( K, J )*DCONJG( A( K, I ) )
 
250
   80             CONTINUE
 
251
                  IF( BETA.EQ.ZERO )THEN
 
252
                     C( I, J ) = TEMP1*DBLE( A( I, I ) ) +
 
253
     $                           ALPHA*TEMP2
 
254
                  ELSE
 
255
                     C( I, J ) = BETA *C( I, J )         +
 
256
     $                           TEMP1*DBLE( A( I, I ) ) +
 
257
     $                           ALPHA*TEMP2
 
258
                  END IF
 
259
   90          CONTINUE
 
260
  100       CONTINUE
 
261
         END IF
 
262
      ELSE
 
263
*
 
264
*        Form  C := alpha*B*A + beta*C.
 
265
*
 
266
         DO 170, J = 1, N
 
267
            TEMP1 = ALPHA*DBLE( A( J, J ) )
 
268
            IF( BETA.EQ.ZERO )THEN
 
269
               DO 110, I = 1, M
 
270
                  C( I, J ) = TEMP1*B( I, J )
 
271
  110          CONTINUE
 
272
            ELSE
 
273
               DO 120, I = 1, M
 
274
                  C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J )
 
275
  120          CONTINUE
 
276
            END IF
 
277
            DO 140, K = 1, J - 1
 
278
               IF( UPPER )THEN
 
279
                  TEMP1 = ALPHA*A( K, J )
 
280
               ELSE
 
281
                  TEMP1 = ALPHA*DCONJG( A( J, K ) )
 
282
               END IF
 
283
               DO 130, I = 1, M
 
284
                  C( I, J ) = C( I, J ) + TEMP1*B( I, K )
 
285
  130          CONTINUE
 
286
  140       CONTINUE
 
287
            DO 160, K = J + 1, N
 
288
               IF( UPPER )THEN
 
289
                  TEMP1 = ALPHA*DCONJG( A( J, K ) )
 
290
               ELSE
 
291
                  TEMP1 = ALPHA*A( K, J )
 
292
               END IF
 
293
               DO 150, I = 1, M
 
294
                  C( I, J ) = C( I, J ) + TEMP1*B( I, K )
 
295
  150          CONTINUE
 
296
  160       CONTINUE
 
297
  170    CONTINUE
 
298
      END IF
 
299
*
 
300
      RETURN
 
301
*
 
302
*     End of ZHEMM .
 
303
*
 
304
      END