Free cookie consent management tool by TermsFeed Policy Generator

source: branches/MathNetNumerics-Exploration-2789/HeuristicLab.Algorithms.DataAnalysis.Experimental/sbart/dsymm.f @ 15457

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

#2789 added Finbarr O'Sullivan smoothing spline code

File size: 9.4 KB
Line 
1      SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2*     .. Scalar Arguments ..
3      DOUBLE PRECISION ALPHA,BETA
4      INTEGER LDA,LDB,LDC,M,N
5      CHARACTER SIDE,UPLO
6*     ..
7*     .. Array Arguments ..
8      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
9*     ..
10*
11*  Purpose
12*  =======
13*
14*  DSYMM  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  - DOUBLE PRECISION.
62*           On entry, ALPHA specifies the scalar alpha.
63*           Unchanged on exit.
64*
65*  A      - DOUBLE PRECISION 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      - DOUBLE PRECISION 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   - DOUBLE PRECISION.
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      - DOUBLE PRECISION 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 PRECISION TEMP1,TEMP2
147      INTEGER I,INFO,J,K,NROWA
148      LOGICAL UPPER
149*     ..
150*     .. Parameters ..
151      DOUBLE PRECISION ONE,ZERO
152      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
153*     ..
154*
155*     Set NROWA as the number of rows of A.
156*
157      IF (LSAME(SIDE,'L')) THEN
158          NROWA = M
159      ELSE
160          NROWA = N
161      END IF
162      UPPER = LSAME(UPLO,'U')
163*
164*     Test the input parameters.
165*
166      INFO = 0
167      IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
168          INFO = 1
169      ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
170          INFO = 2
171      ELSE IF (M.LT.0) THEN
172          INFO = 3
173      ELSE IF (N.LT.0) THEN
174          INFO = 4
175      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
176          INFO = 7
177      ELSE IF (LDB.LT.MAX(1,M)) THEN
178          INFO = 9
179      ELSE IF (LDC.LT.MAX(1,M)) THEN
180          INFO = 12
181      END IF
182      IF (INFO.NE.0) THEN
183          CALL XERBLA('DSYMM ',INFO)
184          RETURN
185      END IF
186*
187*     Quick return if possible.
188*
189      IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
190     +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
191*
192*     And when  alpha.eq.zero.
193*
194      IF (ALPHA.EQ.ZERO) THEN
195          IF (BETA.EQ.ZERO) THEN
196              DO 20 J = 1,N
197                  DO 10 I = 1,M
198                      C(I,J) = ZERO
199   10             CONTINUE
200   20         CONTINUE
201          ELSE
202              DO 40 J = 1,N
203                  DO 30 I = 1,M
204                      C(I,J) = BETA*C(I,J)
205   30             CONTINUE
206   40         CONTINUE
207          END IF
208          RETURN
209      END IF
210*
211*     Start the operations.
212*
213      IF (LSAME(SIDE,'L')) THEN
214*
215*        Form  C := alpha*A*B + beta*C.
216*
217          IF (UPPER) THEN
218              DO 70 J = 1,N
219                  DO 60 I = 1,M
220                      TEMP1 = ALPHA*B(I,J)
221                      TEMP2 = ZERO
222                      DO 50 K = 1,I - 1
223                          C(K,J) = C(K,J) + TEMP1*A(K,I)
224                          TEMP2 = TEMP2 + B(K,J)*A(K,I)
225   50                 CONTINUE
226                      IF (BETA.EQ.ZERO) THEN
227                          C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
228                      ELSE
229                          C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
230     +                             ALPHA*TEMP2
231                      END IF
232   60             CONTINUE
233   70         CONTINUE
234          ELSE
235              DO 100 J = 1,N
236                  DO 90 I = M,1,-1
237                      TEMP1 = ALPHA*B(I,J)
238                      TEMP2 = ZERO
239                      DO 80 K = I + 1,M
240                          C(K,J) = C(K,J) + TEMP1*A(K,I)
241                          TEMP2 = TEMP2 + B(K,J)*A(K,I)
242   80                 CONTINUE
243                      IF (BETA.EQ.ZERO) THEN
244                          C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
245                      ELSE
246                          C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
247     +                             ALPHA*TEMP2
248                      END IF
249   90             CONTINUE
250  100         CONTINUE
251          END IF
252      ELSE
253*
254*        Form  C := alpha*B*A + beta*C.
255*
256          DO 170 J = 1,N
257              TEMP1 = ALPHA*A(J,J)
258              IF (BETA.EQ.ZERO) THEN
259                  DO 110 I = 1,M
260                      C(I,J) = TEMP1*B(I,J)
261  110             CONTINUE
262              ELSE
263                  DO 120 I = 1,M
264                      C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
265  120             CONTINUE
266              END IF
267              DO 140 K = 1,J - 1
268                  IF (UPPER) THEN
269                      TEMP1 = ALPHA*A(K,J)
270                  ELSE
271                      TEMP1 = ALPHA*A(J,K)
272                  END IF
273                  DO 130 I = 1,M
274                      C(I,J) = C(I,J) + TEMP1*B(I,K)
275  130             CONTINUE
276  140         CONTINUE
277              DO 160 K = J + 1,N
278                  IF (UPPER) THEN
279                      TEMP1 = ALPHA*A(J,K)
280                  ELSE
281                      TEMP1 = ALPHA*A(K,J)
282                  END IF
283                  DO 150 I = 1,M
284                      C(I,J) = C(I,J) + TEMP1*B(I,K)
285  150             CONTINUE
286  160         CONTINUE
287  170     CONTINUE
288      END IF
289*
290      RETURN
291*
292*     End of DSYMM .
293*
294      END
Note: See TracBrowser for help on using the repository browser.