Free cookie consent management tool by TermsFeed Policy Generator

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