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	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