Free cookie consent management tool by TermsFeed Policy Generator

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