← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/linalg/dsyr2k.f

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/linalg/dsyr2k.f
 
1
*> \brief \b DSYR2K
 
2
*
 
3
*  =========== DOCUMENTATION ===========
 
4
*
 
5
* Online html documentation available at 
 
6
*            http://www.netlib.org/lapack/explore-html/ 
 
7
*
 
8
*  Definition:
 
9
*  ===========
 
10
*
 
11
*       SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
 
12
* 
 
13
*       .. Scalar Arguments ..
 
14
*       DOUBLE PRECISION ALPHA,BETA
 
15
*       INTEGER K,LDA,LDB,LDC,N
 
16
*       CHARACTER TRANS,UPLO
 
17
*       ..
 
18
*       .. Array Arguments ..
 
19
*       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
 
20
*       ..
 
21
*  
 
22
*
 
23
*> \par Purpose:
 
24
*  =============
 
25
*>
 
26
*> \verbatim
 
27
*>
 
28
*> DSYR2K  performs one of the symmetric rank 2k operations
 
29
*>
 
30
*>    C := alpha*A*B**T + alpha*B*A**T + beta*C,
 
31
*>
 
32
*> or
 
33
*>
 
34
*>    C := alpha*A**T*B + alpha*B**T*A + beta*C,
 
35
*>
 
36
*> where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
 
37
*> and  A and B  are  n by k  matrices  in the  first  case  and  k by n
 
38
*> matrices in the second case.
 
39
*> \endverbatim
 
40
*
 
41
*  Arguments:
 
42
*  ==========
 
43
*
 
44
*> \param[in] UPLO
 
45
*> \verbatim
 
46
*>          UPLO is CHARACTER*1
 
47
*>           On  entry,   UPLO  specifies  whether  the  upper  or  lower
 
48
*>           triangular  part  of the  array  C  is to be  referenced  as
 
49
*>           follows:
 
50
*>
 
51
*>              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
 
52
*>                                  is to be referenced.
 
53
*>
 
54
*>              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
 
55
*>                                  is to be referenced.
 
56
*> \endverbatim
 
57
*>
 
58
*> \param[in] TRANS
 
59
*> \verbatim
 
60
*>          TRANS is CHARACTER*1
 
61
*>           On entry,  TRANS  specifies the operation to be performed as
 
62
*>           follows:
 
63
*>
 
64
*>              TRANS = 'N' or 'n'   C := alpha*A*B**T + alpha*B*A**T +
 
65
*>                                        beta*C.
 
66
*>
 
67
*>              TRANS = 'T' or 't'   C := alpha*A**T*B + alpha*B**T*A +
 
68
*>                                        beta*C.
 
69
*>
 
70
*>              TRANS = 'C' or 'c'   C := alpha*A**T*B + alpha*B**T*A +
 
71
*>                                        beta*C.
 
72
*> \endverbatim
 
73
*>
 
74
*> \param[in] N
 
75
*> \verbatim
 
76
*>          N is INTEGER
 
77
*>           On entry,  N specifies the order of the matrix C.  N must be
 
78
*>           at least zero.
 
79
*> \endverbatim
 
80
*>
 
81
*> \param[in] K
 
82
*> \verbatim
 
83
*>          K is INTEGER
 
84
*>           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
 
85
*>           of  columns  of the  matrices  A and B,  and on  entry  with
 
86
*>           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
 
87
*>           of rows of the matrices  A and B.  K must be at least  zero.
 
88
*> \endverbatim
 
89
*>
 
90
*> \param[in] ALPHA
 
91
*> \verbatim
 
92
*>          ALPHA is DOUBLE PRECISION.
 
93
*>           On entry, ALPHA specifies the scalar alpha.
 
94
*> \endverbatim
 
95
*>
 
96
*> \param[in] A
 
97
*> \verbatim
 
98
*>          A is DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
 
99
*>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
 
100
*>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
 
101
*>           part of the array  A  must contain the matrix  A,  otherwise
 
102
*>           the leading  k by n  part of the array  A  must contain  the
 
103
*>           matrix A.
 
104
*> \endverbatim
 
105
*>
 
106
*> \param[in] LDA
 
107
*> \verbatim
 
108
*>          LDA is INTEGER
 
109
*>           On entry, LDA specifies the first dimension of A as declared
 
110
*>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
 
111
*>           then  LDA must be at least  max( 1, n ), otherwise  LDA must
 
112
*>           be at least  max( 1, k ).
 
113
*> \endverbatim
 
114
*>
 
115
*> \param[in] B
 
116
*> \verbatim
 
117
*>          B is DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
 
118
*>           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
 
119
*>           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
 
120
*>           part of the array  B  must contain the matrix  B,  otherwise
 
121
*>           the leading  k by n  part of the array  B  must contain  the
 
122
*>           matrix B.
 
123
*> \endverbatim
 
124
*>
 
125
*> \param[in] LDB
 
126
*> \verbatim
 
127
*>          LDB is INTEGER
 
128
*>           On entry, LDB specifies the first dimension of B as declared
 
129
*>           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
 
130
*>           then  LDB must be at least  max( 1, n ), otherwise  LDB must
 
131
*>           be at least  max( 1, k ).
 
132
*> \endverbatim
 
133
*>
 
134
*> \param[in] BETA
 
135
*> \verbatim
 
136
*>          BETA is DOUBLE PRECISION.
 
137
*>           On entry, BETA specifies the scalar beta.
 
138
*> \endverbatim
 
139
*>
 
140
*> \param[in,out] C
 
141
*> \verbatim
 
142
*>          C is DOUBLE PRECISION array of DIMENSION ( LDC, n ).
 
143
*>           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
 
144
*>           upper triangular part of the array C must contain the upper
 
145
*>           triangular part  of the  symmetric matrix  and the strictly
 
146
*>           lower triangular part of C is not referenced.  On exit, the
 
147
*>           upper triangular part of the array  C is overwritten by the
 
148
*>           upper triangular part of the updated matrix.
 
149
*>           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
 
150
*>           lower triangular part of the array C must contain the lower
 
151
*>           triangular part  of the  symmetric matrix  and the strictly
 
152
*>           upper triangular part of C is not referenced.  On exit, the
 
153
*>           lower triangular part of the array  C is overwritten by the
 
154
*>           lower triangular part of the updated matrix.
 
155
*> \endverbatim
 
156
*>
 
157
*> \param[in] LDC
 
158
*> \verbatim
 
159
*>          LDC is INTEGER
 
160
*>           On entry, LDC specifies the first dimension of C as declared
 
161
*>           in  the  calling  (sub)  program.   LDC  must  be  at  least
 
162
*>           max( 1, n ).
 
163
*> \endverbatim
 
164
*
 
165
*  Authors:
 
166
*  ========
 
167
*
 
168
*> \author Univ. of Tennessee 
 
169
*> \author Univ. of California Berkeley 
 
170
*> \author Univ. of Colorado Denver 
 
171
*> \author NAG Ltd. 
 
172
*
 
173
*> \date November 2011
 
174
*
 
175
*> \ingroup double_blas_level3
 
176
*
 
177
*> \par Further Details:
 
178
*  =====================
 
179
*>
 
180
*> \verbatim
 
181
*>
 
182
*>  Level 3 Blas routine.
 
183
*>
 
184
*>
 
185
*>  -- Written on 8-February-1989.
 
186
*>     Jack Dongarra, Argonne National Laboratory.
 
187
*>     Iain Duff, AERE Harwell.
 
188
*>     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 
189
*>     Sven Hammarling, Numerical Algorithms Group Ltd.
 
190
*> \endverbatim
 
191
*>
 
192
*  =====================================================================
 
193
      SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
 
194
*
 
195
*  -- Reference BLAS level3 routine (version 3.4.0) --
 
196
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 
197
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 
198
*     November 2011
 
199
*
 
200
*     .. Scalar Arguments ..
 
201
      DOUBLE PRECISION ALPHA,BETA
 
202
      INTEGER K,LDA,LDB,LDC,N
 
203
      CHARACTER TRANS,UPLO
 
204
*     ..
 
205
*     .. Array Arguments ..
 
206
      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
 
207
*     ..
 
208
*
 
209
*  =====================================================================
 
210
*
 
211
*     .. External Functions ..
 
212
      LOGICAL LSAME
 
213
      EXTERNAL LSAME
 
214
*     ..
 
215
*     .. External Subroutines ..
 
216
      EXTERNAL XERBLA
 
217
*     ..
 
218
*     .. Intrinsic Functions ..
 
219
      INTRINSIC MAX
 
220
*     ..
 
221
*     .. Local Scalars ..
 
222
      DOUBLE PRECISION TEMP1,TEMP2
 
223
      INTEGER I,INFO,J,L,NROWA
 
224
      LOGICAL UPPER
 
225
*     ..
 
226
*     .. Parameters ..
 
227
      DOUBLE PRECISION ONE,ZERO
 
228
      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
 
229
*     ..
 
230
*
 
231
*     Test the input parameters.
 
232
*
 
233
      IF (LSAME(TRANS,'N')) THEN
 
234
          NROWA = N
 
235
      ELSE
 
236
          NROWA = K
 
237
      END IF
 
238
      UPPER = LSAME(UPLO,'U')
 
239
*
 
240
      INFO = 0
 
241
      IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
 
242
          INFO = 1
 
243
      ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
 
244
     +         (.NOT.LSAME(TRANS,'T')) .AND.
 
245
     +         (.NOT.LSAME(TRANS,'C'))) THEN
 
246
          INFO = 2
 
247
      ELSE IF (N.LT.0) THEN
 
248
          INFO = 3
 
249
      ELSE IF (K.LT.0) THEN
 
250
          INFO = 4
 
251
      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
 
252
          INFO = 7
 
253
      ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
 
254
          INFO = 9
 
255
      ELSE IF (LDC.LT.MAX(1,N)) THEN
 
256
          INFO = 12
 
257
      END IF
 
258
      IF (INFO.NE.0) THEN
 
259
          CALL XERBLA('DSYR2K',INFO)
 
260
          RETURN
 
261
      END IF
 
262
*
 
263
*     Quick return if possible.
 
264
*
 
265
      IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
 
266
     +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
 
267
*
 
268
*     And when  alpha.eq.zero.
 
269
*
 
270
      IF (ALPHA.EQ.ZERO) THEN
 
271
          IF (UPPER) THEN
 
272
              IF (BETA.EQ.ZERO) THEN
 
273
                  DO 20 J = 1,N
 
274
                      DO 10 I = 1,J
 
275
                          C(I,J) = ZERO
 
276
   10                 CONTINUE
 
277
   20             CONTINUE
 
278
              ELSE
 
279
                  DO 40 J = 1,N
 
280
                      DO 30 I = 1,J
 
281
                          C(I,J) = BETA*C(I,J)
 
282
   30                 CONTINUE
 
283
   40             CONTINUE
 
284
              END IF
 
285
          ELSE
 
286
              IF (BETA.EQ.ZERO) THEN
 
287
                  DO 60 J = 1,N
 
288
                      DO 50 I = J,N
 
289
                          C(I,J) = ZERO
 
290
   50                 CONTINUE
 
291
   60             CONTINUE
 
292
              ELSE
 
293
                  DO 80 J = 1,N
 
294
                      DO 70 I = J,N
 
295
                          C(I,J) = BETA*C(I,J)
 
296
   70                 CONTINUE
 
297
   80             CONTINUE
 
298
              END IF
 
299
          END IF
 
300
          RETURN
 
301
      END IF
 
302
*
 
303
*     Start the operations.
 
304
*
 
305
      IF (LSAME(TRANS,'N')) THEN
 
306
*
 
307
*        Form  C := alpha*A*B**T + alpha*B*A**T + C.
 
308
*
 
309
          IF (UPPER) THEN
 
310
              DO 130 J = 1,N
 
311
                  IF (BETA.EQ.ZERO) THEN
 
312
                      DO 90 I = 1,J
 
313
                          C(I,J) = ZERO
 
314
   90                 CONTINUE
 
315
                  ELSE IF (BETA.NE.ONE) THEN
 
316
                      DO 100 I = 1,J
 
317
                          C(I,J) = BETA*C(I,J)
 
318
  100                 CONTINUE
 
319
                  END IF
 
320
                  DO 120 L = 1,K
 
321
                      IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
 
322
                          TEMP1 = ALPHA*B(J,L)
 
323
                          TEMP2 = ALPHA*A(J,L)
 
324
                          DO 110 I = 1,J
 
325
                              C(I,J) = C(I,J) + A(I,L)*TEMP1 +
 
326
     +                                 B(I,L)*TEMP2
 
327
  110                     CONTINUE
 
328
                      END IF
 
329
  120             CONTINUE
 
330
  130         CONTINUE
 
331
          ELSE
 
332
              DO 180 J = 1,N
 
333
                  IF (BETA.EQ.ZERO) THEN
 
334
                      DO 140 I = J,N
 
335
                          C(I,J) = ZERO
 
336
  140                 CONTINUE
 
337
                  ELSE IF (BETA.NE.ONE) THEN
 
338
                      DO 150 I = J,N
 
339
                          C(I,J) = BETA*C(I,J)
 
340
  150                 CONTINUE
 
341
                  END IF
 
342
                  DO 170 L = 1,K
 
343
                      IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
 
344
                          TEMP1 = ALPHA*B(J,L)
 
345
                          TEMP2 = ALPHA*A(J,L)
 
346
                          DO 160 I = J,N
 
347
                              C(I,J) = C(I,J) + A(I,L)*TEMP1 +
 
348
     +                                 B(I,L)*TEMP2
 
349
  160                     CONTINUE
 
350
                      END IF
 
351
  170             CONTINUE
 
352
  180         CONTINUE
 
353
          END IF
 
354
      ELSE
 
355
*
 
356
*        Form  C := alpha*A**T*B + alpha*B**T*A + C.
 
357
*
 
358
          IF (UPPER) THEN
 
359
              DO 210 J = 1,N
 
360
                  DO 200 I = 1,J
 
361
                      TEMP1 = ZERO
 
362
                      TEMP2 = ZERO
 
363
                      DO 190 L = 1,K
 
364
                          TEMP1 = TEMP1 + A(L,I)*B(L,J)
 
365
                          TEMP2 = TEMP2 + B(L,I)*A(L,J)
 
366
  190                 CONTINUE
 
367
                      IF (BETA.EQ.ZERO) THEN
 
368
                          C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
 
369
                      ELSE
 
370
                          C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
 
371
     +                             ALPHA*TEMP2
 
372
                      END IF
 
373
  200             CONTINUE
 
374
  210         CONTINUE
 
375
          ELSE
 
376
              DO 240 J = 1,N
 
377
                  DO 230 I = J,N
 
378
                      TEMP1 = ZERO
 
379
                      TEMP2 = ZERO
 
380
                      DO 220 L = 1,K
 
381
                          TEMP1 = TEMP1 + A(L,I)*B(L,J)
 
382
                          TEMP2 = TEMP2 + B(L,I)*A(L,J)
 
383
  220                 CONTINUE
 
384
                      IF (BETA.EQ.ZERO) THEN
 
385
                          C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
 
386
                      ELSE
 
387
                          C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
 
388
     +                             ALPHA*TEMP2
 
389
                      END IF
 
390
  230             CONTINUE
 
391
  240         CONTINUE
 
392
          END IF
 
393
      END IF
 
394
*
 
395
      RETURN
 
396
*
 
397
*     End of DSYR2K.
 
398
*
 
399
      END