source: trunk/SOURCES/BLAS/strsm.f

Last change on this file was 4, checked in by dumas, 10 years ago

initial import GRISLI trunk
File size: 12.0 KB
Line 
1      SUBROUTINE STRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
2     $                   B, LDB )
3*     .. Scalar Arguments ..
4      CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
5      INTEGER            M, N, LDA, LDB
6      REAL               ALPHA
7*     .. Array Arguments ..
8      REAL               A( LDA, * ), B( LDB, * )
9*     ..
10*
11*  Purpose
12*  =======
13*
14*  STRSM  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*  Parameters
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  - REAL            .
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      - REAL             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      - REAL             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*     .. External Subroutines ..
134      EXTERNAL           XERBLA
135*     .. Intrinsic Functions ..
136      INTRINSIC          MAX
137*     .. Local Scalars ..
138      LOGICAL            LSIDE, NOUNIT, UPPER
139      INTEGER            I, INFO, J, K, NROWA
140      REAL               TEMP
141*     .. Parameters ..
142      REAL               ONE         , ZERO
143      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
144*     ..
145*     .. Executable Statements ..
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.
160     $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
161         INFO = 1
162      ELSE IF( ( .NOT.UPPER                ).AND.
163     $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
164         INFO = 2
165      ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
166     $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
167     $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
168         INFO = 3
169      ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
170     $         ( .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( 'STRSM ', INFO )
183         RETURN
184      END IF
185*
186*     Quick return if possible.
187*
188      IF( N.EQ.0 )
189     $   RETURN
190*
191*     And when  alpha.eq.zero.
192*
193      IF( ALPHA.EQ.ZERO )THEN
194         DO 20, J = 1, N
195            DO 10, I = 1, M
196               B( I, J ) = ZERO
197   10       CONTINUE
198   20    CONTINUE
199         RETURN
200      END IF
201*
202*     Start the operations.
203*
204      IF( LSIDE )THEN
205         IF( LSAME( TRANSA, 'N' ) )THEN
206*
207*           Form  B := alpha*inv( A )*B.
208*
209            IF( UPPER )THEN
210               DO 60, J = 1, N
211                  IF( ALPHA.NE.ONE )THEN
212                     DO 30, I = 1, M
213                        B( I, J ) = ALPHA*B( I, J )
214   30                CONTINUE
215                  END IF
216                  DO 50, K = M, 1, -1
217                     IF( B( K, J ).NE.ZERO )THEN
218                        IF( NOUNIT )
219     $                     B( K, J ) = B( K, J )/A( K, K )
220                        DO 40, I = 1, K - 1
221                           B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
222   40                   CONTINUE
223                     END IF
224   50             CONTINUE
225   60          CONTINUE
226            ELSE
227               DO 100, J = 1, N
228                  IF( ALPHA.NE.ONE )THEN
229                     DO 70, I = 1, M
230                        B( I, J ) = ALPHA*B( I, J )
231   70                CONTINUE
232                  END IF
233                  DO 90 K = 1, M
234                     IF( B( K, J ).NE.ZERO )THEN
235                        IF( NOUNIT )
236     $                     B( K, J ) = B( K, J )/A( K, K )
237                        DO 80, I = K + 1, M
238                           B( I, J ) = B( I, J ) - B( K, J )*A( I, K )
239   80                   CONTINUE
240                     END IF
241   90             CONTINUE
242  100          CONTINUE
243            END IF
244         ELSE
245*
246*           Form  B := alpha*inv( A' )*B.
247*
248            IF( UPPER )THEN
249               DO 130, J = 1, N
250                  DO 120, I = 1, M
251                     TEMP = ALPHA*B( I, J )
252                     DO 110, K = 1, I - 1
253                        TEMP = TEMP - A( K, I )*B( K, J )
254  110                CONTINUE
255                     IF( NOUNIT )
256     $                  TEMP = TEMP/A( I, I )
257                     B( I, J ) = TEMP
258  120             CONTINUE
259  130          CONTINUE
260            ELSE
261               DO 160, J = 1, N
262                  DO 150, I = M, 1, -1
263                     TEMP = ALPHA*B( I, J )
264                     DO 140, K = I + 1, M
265                        TEMP = TEMP - A( K, I )*B( K, J )
266  140                CONTINUE
267                     IF( NOUNIT )
268     $                  TEMP = TEMP/A( I, I )
269                     B( I, J ) = TEMP
270  150             CONTINUE
271  160          CONTINUE
272            END IF
273         END IF
274      ELSE
275         IF( LSAME( TRANSA, 'N' ) )THEN
276*
277*           Form  B := alpha*B*inv( A ).
278*
279            IF( UPPER )THEN
280               DO 210, J = 1, N
281                  IF( ALPHA.NE.ONE )THEN
282                     DO 170, I = 1, M
283                        B( I, J ) = ALPHA*B( I, J )
284  170                CONTINUE
285                  END IF
286                  DO 190, K = 1, J - 1
287                     IF( A( K, J ).NE.ZERO )THEN
288                        DO 180, I = 1, M
289                           B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
290  180                   CONTINUE
291                     END IF
292  190             CONTINUE
293                  IF( NOUNIT )THEN
294                     TEMP = ONE/A( J, J )
295                     DO 200, I = 1, M
296                        B( I, J ) = TEMP*B( I, J )
297  200                CONTINUE
298                  END IF
299  210          CONTINUE
300            ELSE
301               DO 260, J = N, 1, -1
302                  IF( ALPHA.NE.ONE )THEN
303                     DO 220, I = 1, M
304                        B( I, J ) = ALPHA*B( I, J )
305  220                CONTINUE
306                  END IF
307                  DO 240, K = J + 1, N
308                     IF( A( K, J ).NE.ZERO )THEN
309                        DO 230, I = 1, M
310                           B( I, J ) = B( I, J ) - A( K, J )*B( I, K )
311  230                   CONTINUE
312                     END IF
313  240             CONTINUE
314                  IF( NOUNIT )THEN
315                     TEMP = ONE/A( J, J )
316                     DO 250, I = 1, M
317                       B( I, J ) = TEMP*B( I, J )
318  250                CONTINUE
319                  END IF
320  260          CONTINUE
321            END IF
322         ELSE
323*
324*           Form  B := alpha*B*inv( A' ).
325*
326            IF( UPPER )THEN
327               DO 310, K = N, 1, -1
328                  IF( NOUNIT )THEN
329                     TEMP = ONE/A( K, K )
330                     DO 270, I = 1, M
331                        B( I, K ) = TEMP*B( I, K )
332  270                CONTINUE
333                  END IF
334                  DO 290, J = 1, K - 1
335                     IF( A( J, K ).NE.ZERO )THEN
336                        TEMP = A( J, K )
337                        DO 280, I = 1, M
338                           B( I, J ) = B( I, J ) - TEMP*B( I, K )
339  280                   CONTINUE
340                     END IF
341  290             CONTINUE
342                  IF( ALPHA.NE.ONE )THEN
343                     DO 300, I = 1, M
344                        B( I, K ) = ALPHA*B( I, K )
345  300                CONTINUE
346                  END IF
347  310          CONTINUE
348            ELSE
349               DO 360, K = 1, N
350                  IF( NOUNIT )THEN
351                     TEMP = ONE/A( K, K )
352                     DO 320, I = 1, M
353                        B( I, K ) = TEMP*B( I, K )
354  320                CONTINUE
355                  END IF
356                  DO 340, J = K + 1, N
357                     IF( A( J, K ).NE.ZERO )THEN
358                        TEMP = A( J, K )
359                        DO 330, I = 1, M
360                           B( I, J ) = B( I, J ) - TEMP*B( I, K )
361  330                   CONTINUE
362                     END IF
363  340             CONTINUE
364                  IF( ALPHA.NE.ONE )THEN
365                     DO 350, I = 1, M
366                        B( I, K ) = ALPHA*B( I, K )
367  350                CONTINUE
368                  END IF
369  360          CONTINUE
370            END IF
371         END IF
372      END IF
373*
374      RETURN
375*
376*     End of STRSM .
377*
378      END
Note: See TracBrowser for help on using the repository browser.