Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2789_MathNetNumerics-Exploration/HeuristicLab.Algorithms.DataAnalysis.Experimental/sbart/zsymm.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: 9.4 KB
Line 
1      SUBROUTINE ZSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2*     .. Scalar Arguments ..
3      DOUBLE COMPLEX ALPHA,BETA
4      INTEGER LDA,LDB,LDC,M,N
5      CHARACTER SIDE,UPLO
6*     ..
7*     .. Array Arguments ..
8      DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*)
9*     ..
10*
11*  Purpose
12*  =======
13*
14*  ZSYMM  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 a symmetric matrix and  B and
23*  C are m by n matrices.
24*
25*  Arguments
26*  ==========
27*
28*  SIDE   - CHARACTER*1.
29*           On entry,  SIDE  specifies whether  the  symmetric 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  symmetric  matrix   A  is  to  be
41*           referenced as follows:
42*
43*              UPLO = 'U' or 'u'   Only the upper triangular part of the
44*                                  symmetric matrix is to be referenced.
45*
46*              UPLO = 'L' or 'l'   Only the lower triangular part of the
47*                                  symmetric 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  symmetric 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  symmetric 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  symmetric
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  symmetric 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  symmetric 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  symmetric
85*           matrix and the  strictly upper triangular part of  A  is not
86*           referenced.
87*           Unchanged on exit.
88*
89*  LDA    - INTEGER.
90*           On entry, LDA specifies the first dimension of A as declared
91*           in the  calling (sub) program. When  SIDE = 'L' or 'l'  then
92*           LDA must be at least  max( 1, m ), otherwise  LDA must be at
93*           least max( 1, n ).
94*           Unchanged on exit.
95*
96*  B      - COMPLEX*16       array of DIMENSION ( LDB, n ).
97*           Before entry, the leading  m by n part of the array  B  must
98*           contain the matrix B.
99*           Unchanged on exit.
100*
101*  LDB    - INTEGER.
102*           On entry, LDB specifies the first dimension of B as declared
103*           in  the  calling  (sub)  program.   LDB  must  be  at  least
104*           max( 1, m ).
105*           Unchanged on exit.
106*
107*  BETA   - COMPLEX*16      .
108*           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
109*           supplied as zero then C need not be set on input.
110*           Unchanged on exit.
111*
112*  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
113*           Before entry, the leading  m by n  part of the array  C must
114*           contain the matrix  C,  except when  beta  is zero, in which
115*           case C need not be set on entry.
116*           On exit, the array  C  is overwritten by the  m by n updated
117*           matrix.
118*
119*  LDC    - INTEGER.
120*           On entry, LDC specifies the first dimension of C as declared
121*           in  the  calling  (sub)  program.   LDC  must  be  at  least
122*           max( 1, m ).
123*           Unchanged on exit.
124*
125*
126*  Level 3 Blas routine.
127*
128*  -- Written on 8-February-1989.
129*     Jack Dongarra, Argonne National Laboratory.
130*     Iain Duff, AERE Harwell.
131*     Jeremy Du Croz, Numerical Algorithms Group Ltd.
132*     Sven Hammarling, Numerical Algorithms Group Ltd.
133*
134*
135*     .. External Functions ..
136      LOGICAL LSAME
137      EXTERNAL LSAME
138*     ..
139*     .. External Subroutines ..
140      EXTERNAL XERBLA
141*     ..
142*     .. Intrinsic Functions ..
143      INTRINSIC MAX
144*     ..
145*     .. Local Scalars ..
146      DOUBLE COMPLEX TEMP1,TEMP2
147      INTEGER I,INFO,J,K,NROWA
148      LOGICAL UPPER
149*     ..
150*     .. Parameters ..
151      DOUBLE COMPLEX ONE
152      PARAMETER (ONE= (1.0D+0,0.0D+0))
153      DOUBLE COMPLEX ZERO
154      PARAMETER (ZERO= (0.0D+0,0.0D+0))
155*     ..
156*
157*     Set NROWA as the number of rows of A.
158*
159      IF (LSAME(SIDE,'L')) THEN
160          NROWA = M
161      ELSE
162          NROWA = N
163      END IF
164      UPPER = LSAME(UPLO,'U')
165*
166*     Test the input parameters.
167*
168      INFO = 0
169      IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
170          INFO = 1
171      ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
172          INFO = 2
173      ELSE IF (M.LT.0) THEN
174          INFO = 3
175      ELSE IF (N.LT.0) THEN
176          INFO = 4
177      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
178          INFO = 7
179      ELSE IF (LDB.LT.MAX(1,M)) THEN
180          INFO = 9
181      ELSE IF (LDC.LT.MAX(1,M)) THEN
182          INFO = 12
183      END IF
184      IF (INFO.NE.0) THEN
185          CALL XERBLA('ZSYMM ',INFO)
186          RETURN
187      END IF
188*
189*     Quick return if possible.
190*
191      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
192     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
193*
194*     And when  alpha.eq.zero.
195*
196      IF (ALPHA.EQ.ZERO) THEN
197          IF (BETA.EQ.ZERO) THEN
198              DO 20 J = 1,N
199                  DO 10 I = 1,M
200                      C(I,J) = ZERO
201   10             CONTINUE
202   20         CONTINUE
203          ELSE
204              DO 40 J = 1,N
205                  DO 30 I = 1,M
206                      C(I,J) = BETA*C(I,J)
207   30             CONTINUE
208   40         CONTINUE
209          END IF
210          RETURN
211      END IF
212*
213*     Start the operations.
214*
215      IF (LSAME(SIDE,'L')) THEN
216*
217*        Form  C := alpha*A*B + beta*C.
218*
219          IF (UPPER) THEN
220              DO 70 J = 1,N
221                  DO 60 I = 1,M
222                      TEMP1 = ALPHA*B(I,J)
223                      TEMP2 = ZERO
224                      DO 50 K = 1,I - 1
225                          C(K,J) = C(K,J) + TEMP1*A(K,I)
226                          TEMP2 = TEMP2 + B(K,J)*A(K,I)
227   50                 CONTINUE
228                      IF (BETA.EQ.ZERO) THEN
229                          C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
230                      ELSE
231                          C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
232     +                             ALPHA*TEMP2
233                      END IF
234   60             CONTINUE
235   70         CONTINUE
236          ELSE
237              DO 100 J = 1,N
238                  DO 90 I = M,1,-1
239                      TEMP1 = ALPHA*B(I,J)
240                      TEMP2 = ZERO
241                      DO 80 K = I + 1,M
242                          C(K,J) = C(K,J) + TEMP1*A(K,I)
243                          TEMP2 = TEMP2 + B(K,J)*A(K,I)
244   80                 CONTINUE
245                      IF (BETA.EQ.ZERO) THEN
246                          C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
247                      ELSE
248                          C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
249     +                             ALPHA*TEMP2
250                      END IF
251   90             CONTINUE
252  100         CONTINUE
253          END IF
254      ELSE
255*
256*        Form  C := alpha*B*A + beta*C.
257*
258          DO 170 J = 1,N
259              TEMP1 = ALPHA*A(J,J)
260              IF (BETA.EQ.ZERO) THEN
261                  DO 110 I = 1,M
262                      C(I,J) = TEMP1*B(I,J)
263  110             CONTINUE
264              ELSE
265                  DO 120 I = 1,M
266                      C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
267  120             CONTINUE
268              END IF
269              DO 140 K = 1,J - 1
270                  IF (UPPER) THEN
271                      TEMP1 = ALPHA*A(K,J)
272                  ELSE
273                      TEMP1 = ALPHA*A(J,K)
274                  END IF
275                  DO 130 I = 1,M
276                      C(I,J) = C(I,J) + TEMP1*B(I,K)
277  130             CONTINUE
278  140         CONTINUE
279              DO 160 K = J + 1,N
280                  IF (UPPER) THEN
281                      TEMP1 = ALPHA*A(J,K)
282                  ELSE
283                      TEMP1 = ALPHA*A(K,J)
284                  END IF
285                  DO 150 I = 1,M
286                      C(I,J) = C(I,J) + TEMP1*B(I,K)
287  150             CONTINUE
288  160         CONTINUE
289  170     CONTINUE
290      END IF
291*
292      RETURN
293*
294*     End of ZSYMM .
295*
296      END
Note: See TracBrowser for help on using the repository browser.