trunk/blas/dtrsv.f

      SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
*
************************************************************************
*
*     .. Scalar Arguments ..
      INTEGER            INCX, LDA, N
      CHARACTER*1        DIAG, TRANS, UPLO
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, N ), X( * )
*     ..
*
*  Purpose
*  =======
*
*  DTRSV  solves one of the systems of equations
*
*     A*x = b,   or   A'*x = b,
*
*  where b and x are n element vectors and A is an n by n unit, or
*  non-unit, upper or lower triangular matrix.
*
*  No test for singularity or near-singularity is included in this
*  routine. Such tests must be performed before calling this routine.
*
*  Parameters
*  ==========
*
*  UPLO   - CHARACTER*1.
*           On entry, UPLO specifies whether the matrix is an upper or
*           lower triangular matrix as follows:
*
*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
*
*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
*
*           Unchanged on exit.
*
*  TRANS  - CHARACTER*1.
*           On entry, TRANS specifies the equations to be solved as
*           follows:
*
*              TRANS = 'N' or 'n'   A*x = b.
*
*              TRANS = 'T' or 't'   A'*x = b.
*
*              TRANS = 'C' or 'c'   A'*x = b.
*
*           Unchanged on exit.
*
*  DIAG   - CHARACTER*1.
*           On entry, DIAG specifies whether or not A is unit
*           triangular as follows:
*
*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
*
*              DIAG = 'N' or 'n'   A is not assumed to be unit
*                                  triangular.
*
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry, N specifies the order of the matrix A.
*           N must be at least zero.
*           Unchanged on exit.
*
*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
*           Before entry with  UPLO = 'U' or 'u', the leading n by n
*           upper triangular part of the array A must contain the upper
*           triangular matrix and the strictly lower triangular part of
*           A is not referenced.
*           Before entry with UPLO = 'L' or 'l', the leading n by n
*           lower triangular part of the array A must contain the lower
*           triangular matrix and the strictly upper triangular part of
*           A is not referenced.
*           Note that when  DIAG = 'U' or 'u', the diagonal elements of
*           A are not referenced either, but are assumed to be unity.
*           Unchanged on exit.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program. LDA must be at least
*           max( 1, n ).
*           Unchanged on exit.
*
*  X      - DOUBLE PRECISION array of dimension at least
*           ( 1 + ( n - 1 )*abs( INCX ) ).
*           Before entry, the incremented array X must contain the n
*           element right-hand side vector b. On exit, X is overwritten
*           with the solution vector x.
*
*  INCX   - INTEGER.
*           On entry, INCX specifies the increment for the elements of
*           X. INCX must not be zero.
*           Unchanged on exit.
*
*
*  Level 2 Blas routine.
*
*  -- Written on 22-October-1986.
*     Jack Dongarra, Argonne National Lab.
*     Jeremy Du Croz, Nag Central Office.
*     Sven Hammarling, Nag Central Office.
*     Richard Hanson, Sandia National Labs.
*
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER        ( ZERO = 0.0D+0 )
*     .. Local Scalars ..
      DOUBLE PRECISION   TEMP
      INTEGER            I, INFO, IX, J, JX, KX
      LOGICAL            NOUNIT
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. External Subroutines ..
      EXTERNAL           XERBLA
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters.
*
      INFO = 0
      IF     ( .NOT.LSAME( UPLO , 'U' ).AND.
     $         .NOT.LSAME( UPLO , 'L' )      )THEN
     INFO = 1
      ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
     $         .NOT.LSAME( TRANS, 'T' ).AND.
     $         .NOT.LSAME( TRANS, 'C' )      )THEN
     INFO = 2
      ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
     $         .NOT.LSAME( DIAG , 'N' )      )THEN
     INFO = 3
      ELSE IF( N.LT.0 )THEN
     INFO = 4
      ELSE IF( LDA.LT.MAX( 1, N ) )THEN
     INFO = 6
      ELSE IF( INCX.EQ.0 )THEN
     INFO = 8
      END IF
      IF( INFO.NE.0 )THEN
     CALL XERBLA( 'DTRSV ', INFO )
     RETURN
      END IF
*
*     Quick return if possible.
*
      IF( N.EQ.0 )
     $   RETURN
*
      NOUNIT = LSAME( DIAG, 'N' )
*
*     Set up the start point in X if the increment is not unity. This
*     will be  ( N - 1 )*INCX  too small for descending loops.
*
      IF( INCX.LE.0 )THEN
     KX = 1 - ( N - 1 )*INCX
      ELSE IF( INCX.NE.1 )THEN
     KX = 1
      END IF
*
*     Start the operations. In this version the elements of A are
*     accessed sequentially with one pass through A.
*
      IF( LSAME( TRANS, 'N' ) )THEN
*
*        Form  x := inv( A )*x.
*
     IF( LSAME( UPLO, 'U' ) )THEN
        IF( INCX.EQ.1 )THEN
           DO 20, J = N, 1, -1
          IF( X( J ).NE.ZERO )THEN
             IF( NOUNIT )
     $                  X( J ) = X( J )/A( J, J )
             TEMP = X( J )
             DO 10, I = J - 1, 1, -1
            X( I ) = X( I ) - TEMP*A( I, J )
   10                CONTINUE
          END IF
   20          CONTINUE
        ELSE
           JX = KX + ( N - 1 )*INCX
           DO 40, J = N, 1, -1
          IF( X( JX ).NE.ZERO )THEN
             IF( NOUNIT )
     $                  X( JX ) = X( JX )/A( J, J )
             TEMP = X( JX )
             IX   = JX
             DO 30, I = J - 1, 1, -1
            IX      = IX      - INCX
            X( IX ) = X( IX ) - TEMP*A( I, J )
   30                CONTINUE
          END IF
          JX = JX - INCX
   40          CONTINUE
        END IF
     ELSE
        IF( INCX.EQ.1 )THEN
           DO 60, J = 1, N
          IF( X( J ).NE.ZERO )THEN
             IF( NOUNIT )
     $                  X( J ) = X( J )/A( J, J )
             TEMP = X( J )
             DO 50, I = J + 1, N
            X( I ) = X( I ) - TEMP*A( I, J )
   50                CONTINUE
          END IF
   60          CONTINUE
        ELSE
           JX = KX
           DO 80, J = 1, N
          IF( X( JX ).NE.ZERO )THEN
             IF( NOUNIT )
     $                  X( JX ) = X( JX )/A( J, J )
             TEMP = X( JX )
             IX   = JX
             DO 70, I = J + 1, N
            IX      = IX      + INCX
            X( IX ) = X( IX ) - TEMP*A( I, J )
   70                CONTINUE
          END IF
          JX = JX + INCX
   80          CONTINUE
        END IF
     END IF
      ELSE
*
*        Form  x := inv( A' )*x.
*
     IF( LSAME( UPLO, 'U' ) )THEN
        IF( INCX.EQ.1 )THEN
           DO 100, J = 1, N
          TEMP = X( J )
          DO 90, I = 1, J - 1
             TEMP = TEMP - A( I, J )*X( I )
   90             CONTINUE
          IF( NOUNIT )
     $               TEMP = TEMP/A( J, J )
          X( J ) = TEMP
  100          CONTINUE
        ELSE
           JX = KX
           DO 120, J = 1, N
          TEMP = X( JX )
          IX   = KX
          DO 110, I = 1, J - 1
             TEMP = TEMP - A( I, J )*X( IX )
             IX   = IX   + INCX
  110             CONTINUE
          IF( NOUNIT )
     $               TEMP = TEMP/A( J, J )
          X( JX ) = TEMP
          JX      = JX   + INCX
  120          CONTINUE
        END IF
     ELSE
        IF( INCX.EQ.1 )THEN
           DO 140, J = N, 1, -1
          TEMP = X( J )
          DO 130, I = N, J + 1, -1
             TEMP = TEMP - A( I, J )*X( I )
  130             CONTINUE
          IF( NOUNIT )
     $               TEMP = TEMP/A( J, J )
          X( J ) = TEMP
  140          CONTINUE
        ELSE
           KX = KX + ( N - 1 )*INCX
           JX = KX
           DO 160, J = N, 1, -1
          TEMP = X( JX )
          IX   = KX
          DO 150, I = N, J + 1, -1
             TEMP = TEMP - A( I, J )*X( IX )
             IX   = IX   - INCX
  150             CONTINUE
          IF( NOUNIT )
     $               TEMP = TEMP/A( J, J )
          X( JX ) = TEMP
          JX      = JX   - INCX
  160          CONTINUE
        END IF
     END IF
      END IF
*
      RETURN
*
*     End of DTRSV .
*
      END
1	aw0a	1	SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
2			*
3			************************************************************************
4			*
5			* .. Scalar Arguments ..
6			INTEGER INCX, LDA, N
7			CHARACTER*1 DIAG, TRANS, UPLO
8			* .. Array Arguments ..
9			DOUBLE PRECISION A( LDA, N ), X( * )
10			* ..
11			*
12			* Purpose
13			* =======
14			*
15			* DTRSV solves one of the systems of equations
16			*
17			* Ax = b, or A'x = b,
18			*
19			* where b and x are n element vectors and A is an n by n unit, or
20			* non-unit, upper or lower triangular matrix.
21			*
22			* No test for singularity or near-singularity is included in this
23			* routine. Such tests must be performed before calling this routine.
24			*
25			* Parameters
26			* ==========
27			*
28			* UPLO - CHARACTER*1.
29			* On entry, UPLO specifies whether the matrix is an upper or
30			* lower triangular matrix as follows:
31			*
32			* UPLO = 'U' or 'u' A is an upper triangular matrix.
33			*
34			* UPLO = 'L' or 'l' A is a lower triangular matrix.
35			*
36			* Unchanged on exit.
37			*
38			* TRANS - CHARACTER*1.
39			* On entry, TRANS specifies the equations to be solved as
40			* follows:
41			*
42			* TRANS = 'N' or 'n' A*x = b.
43			*
44			* TRANS = 'T' or 't' A'*x = b.
45			*
46			* TRANS = 'C' or 'c' A'*x = b.
47			*
48			* Unchanged on exit.
49			*
50			* DIAG - CHARACTER*1.
51			* On entry, DIAG specifies whether or not A is unit
52			* triangular as follows:
53			*
54			* DIAG = 'U' or 'u' A is assumed to be unit triangular.
55			*
56			* DIAG = 'N' or 'n' A is not assumed to be unit
57			* triangular.
58			*
59			* Unchanged on exit.
60			*
61			* N - INTEGER.
62			* On entry, N specifies the order of the matrix A.
63			* N must be at least zero.
64			* Unchanged on exit.
65			*
66			* A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
67			* Before entry with UPLO = 'U' or 'u', the leading n by n
68			* upper triangular part of the array A must contain the upper
69			* triangular matrix and the strictly lower triangular part of
70			* A is not referenced.
71			* Before entry with UPLO = 'L' or 'l', the leading n by n
72			* lower triangular part of the array A must contain the lower
73			* triangular matrix and the strictly upper triangular part of
74			* A is not referenced.
75			* Note that when DIAG = 'U' or 'u', the diagonal elements of
76			* A are not referenced either, but are assumed to be unity.
77			* Unchanged on exit.
78			*
79			* LDA - INTEGER.
80			* On entry, LDA specifies the first dimension of A as declared
81			* in the calling (sub) program. LDA must be at least
82			* max( 1, n ).
83			* Unchanged on exit.
84			*
85			* X - DOUBLE PRECISION array of dimension at least
86			* ( 1 + ( n - 1 )*abs( INCX ) ).
87			* Before entry, the incremented array X must contain the n
88			* element right-hand side vector b. On exit, X is overwritten
89			* with the solution vector x.
90			*
91			* INCX - INTEGER.
92			* On entry, INCX specifies the increment for the elements of
93			* X. INCX must not be zero.
94			* Unchanged on exit.
95			*
96			*
97			* Level 2 Blas routine.
98			*
99			* -- Written on 22-October-1986.
100			* Jack Dongarra, Argonne National Lab.
101			* Jeremy Du Croz, Nag Central Office.
102			* Sven Hammarling, Nag Central Office.
103			* Richard Hanson, Sandia National Labs.
104			*
105			*
106			* .. Parameters ..
107			DOUBLE PRECISION ZERO
108			PARAMETER ( ZERO = 0.0D+0 )
109			* .. Local Scalars ..
110			DOUBLE PRECISION TEMP
111			INTEGER I, INFO, IX, J, JX, KX
112			LOGICAL NOUNIT
113			* .. External Functions ..
114			LOGICAL LSAME
115			EXTERNAL LSAME
116			* .. External Subroutines ..
117			EXTERNAL XERBLA
118			* .. Intrinsic Functions ..
119			INTRINSIC MAX
120			* ..
121			* .. Executable Statements ..
122			*
123			* Test the input parameters.
124			*
125			INFO = 0
126			IF ( .NOT.LSAME( UPLO , 'U' ).AND.
127			$ .NOT.LSAME( UPLO , 'L' ) )THEN
128			INFO = 1
129			ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.
130			$ .NOT.LSAME( TRANS, 'T' ).AND.
131			$ .NOT.LSAME( TRANS, 'C' ) )THEN
132			INFO = 2
133			ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.
134			$ .NOT.LSAME( DIAG , 'N' ) )THEN
135			INFO = 3
136			ELSE IF( N.LT.0 )THEN
137			INFO = 4
138			ELSE IF( LDA.LT.MAX( 1, N ) )THEN
139			INFO = 6
140			ELSE IF( INCX.EQ.0 )THEN
141			INFO = 8
142			END IF
143			IF( INFO.NE.0 )THEN
144			CALL XERBLA( 'DTRSV ', INFO )
145			RETURN
146			END IF
147			*
148			* Quick return if possible.
149			*
150			IF( N.EQ.0 )
151			$ RETURN
152			*
153			NOUNIT = LSAME( DIAG, 'N' )
154			*
155			* Set up the start point in X if the increment is not unity. This
156			* will be ( N - 1 )*INCX too small for descending loops.
157			*
158			IF( INCX.LE.0 )THEN
159			KX = 1 - ( N - 1 )*INCX
160			ELSE IF( INCX.NE.1 )THEN
161			KX = 1
162			END IF
163			*
164			* Start the operations. In this version the elements of A are
165			* accessed sequentially with one pass through A.
166			*
167			IF( LSAME( TRANS, 'N' ) )THEN
168			*
169			* Form x := inv( A )*x.
170			*
171			IF( LSAME( UPLO, 'U' ) )THEN
172			IF( INCX.EQ.1 )THEN
173			DO 20, J = N, 1, -1
174			IF( X( J ).NE.ZERO )THEN
175			IF( NOUNIT )
176			$ X( J ) = X( J )/A( J, J )
177			TEMP = X( J )
178			DO 10, I = J - 1, 1, -1
179			X( I ) = X( I ) - TEMP*A( I, J )
180			10 CONTINUE
181			END IF
182			20 CONTINUE
183			ELSE
184			JX = KX + ( N - 1 )*INCX
185			DO 40, J = N, 1, -1
186			IF( X( JX ).NE.ZERO )THEN
187			IF( NOUNIT )
188			$ X( JX ) = X( JX )/A( J, J )
189			TEMP = X( JX )
190			IX = JX
191			DO 30, I = J - 1, 1, -1
192			IX = IX - INCX
193			X( IX ) = X( IX ) - TEMP*A( I, J )
194			30 CONTINUE
195			END IF
196			JX = JX - INCX
197			40 CONTINUE
198			END IF
199			ELSE
200			IF( INCX.EQ.1 )THEN
201			DO 60, J = 1, N
202			IF( X( J ).NE.ZERO )THEN
203			IF( NOUNIT )
204			$ X( J ) = X( J )/A( J, J )
205			TEMP = X( J )
206			DO 50, I = J + 1, N
207			X( I ) = X( I ) - TEMP*A( I, J )
208			50 CONTINUE
209			END IF
210			60 CONTINUE
211			ELSE
212			JX = KX
213			DO 80, J = 1, N
214			IF( X( JX ).NE.ZERO )THEN
215			IF( NOUNIT )
216			$ X( JX ) = X( JX )/A( J, J )
217			TEMP = X( JX )
218			IX = JX
219			DO 70, I = J + 1, N
220			IX = IX + INCX
221			X( IX ) = X( IX ) - TEMP*A( I, J )
222			70 CONTINUE
223			END IF
224			JX = JX + INCX
225			80 CONTINUE
226			END IF
227			END IF
228			ELSE
229			*
230			* Form x := inv( A' )*x.
231			*
232			IF( LSAME( UPLO, 'U' ) )THEN
233			IF( INCX.EQ.1 )THEN
234			DO 100, J = 1, N
235			TEMP = X( J )
236			DO 90, I = 1, J - 1
237			TEMP = TEMP - A( I, J )*X( I )
238			90 CONTINUE
239			IF( NOUNIT )
240			$ TEMP = TEMP/A( J, J )
241			X( J ) = TEMP
242			100 CONTINUE
243			ELSE
244			JX = KX
245			DO 120, J = 1, N
246			TEMP = X( JX )
247			IX = KX
248			DO 110, I = 1, J - 1
249			TEMP = TEMP - A( I, J )*X( IX )
250			IX = IX + INCX
251			110 CONTINUE
252			IF( NOUNIT )
253			$ TEMP = TEMP/A( J, J )
254			X( JX ) = TEMP
255			JX = JX + INCX
256			120 CONTINUE
257			END IF
258			ELSE
259			IF( INCX.EQ.1 )THEN
260			DO 140, J = N, 1, -1
261			TEMP = X( J )
262			DO 130, I = N, J + 1, -1
263			TEMP = TEMP - A( I, J )*X( I )
264			130 CONTINUE
265			IF( NOUNIT )
266			$ TEMP = TEMP/A( J, J )
267			X( J ) = TEMP
268			140 CONTINUE
269			ELSE
270			KX = KX + ( N - 1 )*INCX
271			JX = KX
272			DO 160, J = N, 1, -1
273			TEMP = X( JX )
274			IX = KX
275			DO 150, I = N, J + 1, -1
276			TEMP = TEMP - A( I, J )*X( IX )
277			IX = IX - INCX
278			150 CONTINUE
279			IF( NOUNIT )
280			$ TEMP = TEMP/A( J, J )
281			X( JX ) = TEMP
282			JX = JX - INCX
283			160 CONTINUE
284			END IF
285			END IF
286			END IF
287			*
288			RETURN
289			*
290			* End of DTRSV .
291			*
292			END