1 : C DTRSV SOURCE CHAT 06/03/29 21:18:59 5360 2 : SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8 (A-H,O-Z) 5 : -INC CCOPTIO 6 : * .. Scalar Arguments .. 7 : INTEGER INCX, LDA, N 8 : CHARACTER*1 DIAG, TRANS, UPLO 9 : * .. Array Arguments .. 10 : REAL*8 A( LDA, * ), X( * ) 11 : * .. 12 : * 13 : * Purpose 14 : * ======= 15 : * 16 : * DTRSV solves one of the systems of equations 17 : * 18 : * A*x = b, or A'*x = b, 19 : * 20 : * where b and x are n element vectors and A is an n by n unit, or 21 : * non-unit, upper or lower triangular matrix. 22 : * 23 : * 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 80 : * in the calling (sub) program. LDA must be at least 81 : * max( 1, n ). 82 : * Unchanged on exit. 83 : * 84 : * X - DOUBLE PRECISION array of dimension at least 85 : * ( 1 + ( n - 1 )*abs( INCX ) ). 86 : * Before entry, the incremented array X must contain the n 87 : * element right-hand side vector b. On exit, X is overwritten 88 : * with the solution vector x. 89 : * 90 : * INCX - INTEGER. 91 : * On entry, INCX specifies the increment for the elements of 92 : * X. INCX must not be zero. 93 : * Unchanged on exit. 94 : * 95 : * 96 : * Level 2 Blas routine. 97 : * 98 : * -- Written on 22-October-1986. 99 : * Jack Dongarra, Argonne National Lab. 100 : * Jeremy Du Croz, Nag Central Office. 101 : * Sven Hammarling, Nag Central Office. 102 : * Richard Hanson, Sandia National Labs. 103 : C modified 5/11/98 : double precision -> real*8 104 : * commented INTRINSIC (not esope compatible) 105 : C added error handling 106 : * 107 : * 108 : * .. Parameters .. 109 : REAL*8 ZERO 110 : PARAMETER ( ZERO = 0.0D+0 ) 111 : * .. Local Scalars .. 112 : REAL*8 TEMP 113 : INTEGER I, INFO, IX, J, JX, KX 114 : LOGICAL NOUNIT 115 : * .. External Functions .. 116 : LOGICAL LSAME 117 : EXTERNAL LSAME 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 : GOTO 9999 145 : END IF 146 : * 147 : * Quick return if possible. 148 : * 149 : IF( N.EQ.0 ) 150 : $ RETURN 151 : * 152 : NOUNIT = LSAME( DIAG, 'N' ) 153 : * 154 : * Set up the start point in X if the increment is not unity. This 155 : * will be ( N - 1 )*INCX too small for descending loops. 156 : * 157 : IF( INCX.LE.0 )THEN 158 : KX = 1 - ( N - 1 )*INCX 159 : ELSE IF( INCX.NE.1 )THEN 160 : KX = 1 161 : END IF 162 : * 163 : * Start the operations. In this version the elements of A are 164 : * accessed sequentially with one pass through A. 165 : * 166 : IF( LSAME( TRANS, 'N' ) )THEN 167 : * 168 : * Form x := inv( A )*x. 169 : * 170 : IF( LSAME( UPLO, 'U' ) )THEN 171 : IF( INCX.EQ.1 )THEN 172 : DO 20, J = N, 1, -1 173 : IF( X( J ).NE.ZERO )THEN 174 : IF( NOUNIT ) 175 : $ X( J ) = X( J )/A( J, J ) 176 : TEMP = X( J ) 177 : DO 10, I = J - 1, 1, -1 178 : X( I ) = X( I ) - TEMP*A( I, J ) 179 : 10 CONTINUE 180 : END IF 181 : 20 CONTINUE 182 : ELSE 183 : JX = KX + ( N - 1 )*INCX 184 : DO 40, J = N, 1, -1 185 : IF( X( JX ).NE.ZERO )THEN 186 : IF( NOUNIT ) 187 : $ X( JX ) = X( JX )/A( J, J ) 188 : TEMP = X( JX ) 189 : IX = JX 190 : DO 30, I = J - 1, 1, -1 191 : IX = IX - INCX 192 : X( IX ) = X( IX ) - TEMP*A( I, J ) 193 : 30 CONTINUE 194 : END IF 195 : JX = JX - INCX 196 : 40 CONTINUE 197 : END IF 198 : ELSE 199 : IF( INCX.EQ.1 )THEN 200 : DO 60, J = 1, N 201 : IF( X( J ).NE.ZERO )THEN 202 : IF( NOUNIT ) 203 : $ X( J ) = X( J )/A( J, J ) 204 : TEMP = X( J ) 205 : DO 50, I = J + 1, N 206 : X( I ) = X( I ) - TEMP*A( I, J ) 207 : 50 CONTINUE 208 : END IF 209 : 60 CONTINUE 210 : ELSE 211 : JX = KX 212 : DO 80, J = 1, N 213 : IF( X( JX ).NE.ZERO )THEN 214 : IF( NOUNIT ) 215 : $ X( JX ) = X( JX )/A( J, J ) 216 : TEMP = X( JX ) 217 : IX = JX 218 : DO 70, I = J + 1, N 219 : IX = IX + INCX 220 : X( IX ) = X( IX ) - TEMP*A( I, J ) 221 : 70 CONTINUE 222 : END IF 223 : JX = JX + INCX 224 : 80 CONTINUE 225 : END IF 226 : END IF 227 : ELSE 228 : * 229 : * Form x := inv( A' )*x. 230 : * 231 : IF( LSAME( UPLO, 'U' ) )THEN 232 : IF( INCX.EQ.1 )THEN 233 : DO 100, J = 1, N 234 : TEMP = X( J ) 235 : DO 90, I = 1, J - 1 236 : TEMP = TEMP - A( I, J )*X( I ) 237 : 90 CONTINUE 238 : IF( NOUNIT ) 239 : $ TEMP = TEMP/A( J, J ) 240 : X( J ) = TEMP 241 : 100 CONTINUE 242 : ELSE 243 : JX = KX 244 : DO 120, J = 1, N 245 : TEMP = X( JX ) 246 : IX = KX 247 : DO 110, I = 1, J - 1 248 : TEMP = TEMP - A( I, J )*X( IX ) 249 : IX = IX + INCX 250 : 110 CONTINUE 251 : IF( NOUNIT ) 252 : $ TEMP = TEMP/A( J, J ) 253 : X( JX ) = TEMP 254 : JX = JX + INCX 255 : 120 CONTINUE 256 : END IF 257 : ELSE 258 : IF( INCX.EQ.1 )THEN 259 : DO 140, J = N, 1, -1 260 : TEMP = X( J ) 261 : DO 130, I = N, J + 1, -1 262 : TEMP = TEMP - A( I, J )*X( I ) 263 : 130 CONTINUE 264 : IF( NOUNIT ) 265 : $ TEMP = TEMP/A( J, J ) 266 : X( J ) = TEMP 267 : 140 CONTINUE 268 : ELSE 269 : KX = KX + ( N - 1 )*INCX 270 : JX = KX 271 : DO 160, J = N, 1, -1 272 : TEMP = X( JX ) 273 : IX = KX 274 : DO 150, I = N, J + 1, -1 275 : TEMP = TEMP - A( I, J )*X( IX ) 276 : IX = IX - INCX 277 : 150 CONTINUE 278 : IF( NOUNIT ) 279 : $ TEMP = TEMP/A( J, J ) 280 : X( JX ) = TEMP 281 : JX = JX - INCX 282 : 160 CONTINUE 283 : END IF 284 : END IF 285 : END IF 286 : * 287 : C 288 : C Normal termination 289 : C 290 : RETURN 291 : C 292 : C Error handling 293 : C 294 : 9999 CONTINUE 295 : WRITE(IOIMP,*) 'Error number :',INFO 296 : WRITE(IOIMP,*) 'in subroutine dtrsv.' 296 : call erreur(21) 297 : RETURN 298 : * 299 : * End of DTRSV . 300 : * 301 : END 302 : 303 : 304 : 305 : 306 : 307 : 308 :
© Cast3M 2003 - All rights reserved.
Disclaimer