Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2789_MathNetNumerics-Exploration/HeuristicLab.Algorithms.DataAnalysis.Experimental/sbart/zherk.f @ 17185

Last change on this file since 17185 was 15457, checked in by gkronber, 7 years ago

#2789 added Finbarr O'Sullivan smoothing spline code

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