CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
dgesv.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dgesv (N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 
subroutine dgetf2 (M, N, A, LDA, IPIV, INFO)
 
subroutine dgetrf (M, N, A, LDA, IPIV, INFO)
 
subroutine dgetrs (TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
 
subroutine dlaswp (N, A, LDA, K1, K2, IPIV, INCX)
 
integer function ieeeck (ISPEC, ZERO, ONE)
 
integer function ilaenv (ISPEC, NAME, OPTS, N1, N2, N3, N4)
 
integer function idamax (n, dx, incx)
 
subroutine dger (M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
 
subroutine dswap (n, dx, incx, dy, incy)
 
subroutine dscal (n, da, dx, incx)
 
subroutine dtrsm (SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
 
subroutine dgemm (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
 

Function/Subroutine Documentation

◆ dgemm()

subroutine dgemm ( character*1  TRANSA,
character*1  TRANSB,
integer  M,
integer  N,
integer  K,
double precision  ALPHA,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision  BETA,
double precision, dimension( ldc, * )  C,
integer  LDC 
)
2170 * .. Scalar Arguments ..
2171  CHARACTER*1 transa, transb
2172  INTEGER m, n, k, lda, ldb, ldc
2173  DOUBLE PRECISION alpha, beta
2174 * .. Array Arguments ..
2175  DOUBLE PRECISION a( lda, * ), b( ldb, * ), c( ldc, * )
2176 * ..
2177 *
2178 * Purpose
2179 * =======
2180 *
2181 * DGEMM performs one of the matrix-matrix operations
2182 *
2183 * C := alpha*op( A )*op( B ) + beta*C,
2184 *
2185 * where op( X ) is one of
2186 *
2187 * op( X ) = X or op( X ) = X',
2188 *
2189 * alpha and beta are scalars, and A, B and C are matrices, with op( A )
2190 * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
2191 *
2192 * Parameters
2193 * ==========
2194 *
2195 * TRANSA - CHARACTER*1.
2196 * On entry, TRANSA specifies the form of op( A ) to be used in
2197 * the matrix multiplication as follows:
2198 *
2199 * TRANSA = 'N' or 'n', op( A ) = A.
2200 *
2201 * TRANSA = 'T' or 't', op( A ) = A'.
2202 *
2203 * TRANSA = 'C' or 'c', op( A ) = A'.
2204 *
2205 * Unchanged on exit.
2206 *
2207 * TRANSB - CHARACTER*1.
2208 * On entry, TRANSB specifies the form of op( B ) to be used in
2209 * the matrix multiplication as follows:
2210 *
2211 * TRANSB = 'N' or 'n', op( B ) = B.
2212 *
2213 * TRANSB = 'T' or 't', op( B ) = B'.
2214 *
2215 * TRANSB = 'C' or 'c', op( B ) = B'.
2216 *
2217 * Unchanged on exit.
2218 *
2219 * M - INTEGER.
2220 * On entry, M specifies the number of rows of the matrix
2221 * op( A ) and of the matrix C. M must be at least zero.
2222 * Unchanged on exit.
2223 *
2224 * N - INTEGER.
2225 * On entry, N specifies the number of columns of the matrix
2226 * op( B ) and the number of columns of the matrix C. N must be
2227 * at least zero.
2228 * Unchanged on exit.
2229 *
2230 * K - INTEGER.
2231 * On entry, K specifies the number of columns of the matrix
2232 * op( A ) and the number of rows of the matrix op( B ). K must
2233 * be at least zero.
2234 * Unchanged on exit.
2235 *
2236 * ALPHA - DOUBLE PRECISION.
2237 * On entry, ALPHA specifies the scalar alpha.
2238 * Unchanged on exit.
2239 *
2240 * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
2241 * k when TRANSA = 'N' or 'n', and is m otherwise.
2242 * Before entry with TRANSA = 'N' or 'n', the leading m by k
2243 * part of the array A must contain the matrix A, otherwise
2244 * the leading k by m part of the array A must contain the
2245 * matrix A.
2246 * Unchanged on exit.
2247 *
2248 * LDA - INTEGER.
2249 * On entry, LDA specifies the first dimension of A as declared
2250 * in the calling (sub) program. When TRANSA = 'N' or 'n' then
2251 * LDA must be at least max( 1, m ), otherwise LDA must be at
2252 * least max( 1, k ).
2253 * Unchanged on exit.
2254 *
2255 * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
2256 * n when TRANSB = 'N' or 'n', and is k otherwise.
2257 * Before entry with TRANSB = 'N' or 'n', the leading k by n
2258 * part of the array B must contain the matrix B, otherwise
2259 * the leading n by k part of the array B must contain the
2260 * matrix B.
2261 * Unchanged on exit.
2262 *
2263 * LDB - INTEGER.
2264 * On entry, LDB specifies the first dimension of B as declared
2265 * in the calling (sub) program. When TRANSB = 'N' or 'n' then
2266 * LDB must be at least max( 1, k ), otherwise LDB must be at
2267 * least max( 1, n ).
2268 * Unchanged on exit.
2269 *
2270 * BETA - DOUBLE PRECISION.
2271 * On entry, BETA specifies the scalar beta. When BETA is
2272 * supplied as zero then C need not be set on input.
2273 * Unchanged on exit.
2274 *
2275 * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
2276 * Before entry, the leading m by n part of the array C must
2277 * contain the matrix C, except when beta is zero, in which
2278 * case C need not be set on entry.
2279 * On exit, the array C is overwritten by the m by n matrix
2280 * ( alpha*op( A )*op( B ) + beta*C ).
2281 *
2282 * LDC - INTEGER.
2283 * On entry, LDC specifies the first dimension of C as declared
2284 * in the calling (sub) program. LDC must be at least
2285 * max( 1, m ).
2286 * Unchanged on exit.
2287 *
2288 *
2289 * Level 3 Blas routine.
2290 *
2291 * -- Written on 8-February-1989.
2292 * Jack Dongarra, Argonne National Laboratory.
2293 * Iain Duff, AERE Harwell.
2294 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
2295 * Sven Hammarling, Numerical Algorithms Group Ltd.
2296 *
2297 *
2298 * .. External Functions ..
2299  LOGICAL lsame
2300  EXTERNAL lsame
2301 * .. External Subroutines ..
2302  EXTERNAL xerbla
2303 * .. Intrinsic Functions ..
2304  INTRINSIC max
2305 * .. Local Scalars ..
2306  LOGICAL nota, notb
2307  INTEGER i, info, j, l, ncola, nrowa, nrowb
2308  DOUBLE PRECISION temp
2309 * .. Parameters ..
2310  DOUBLE PRECISION one , zero
2311  parameter( one = 1.0d+0, zero = 0.0d+0 )
2312 * ..
2313 * .. Executable Statements ..
2314 *
2315 * Set NOTA and NOTB as true if A and B respectively are not
2316 * transposed and set NROWA, NCOLA and NROWB as the number of rows
2317 * and columns of A and the number of rows of B respectively.
2318 *
2319  nota = lsame( transa, 'N' )
2320  notb = lsame( transb, 'N' )
2321  IF( nota )THEN
2322  nrowa = m
2323  ncola = k
2324  ELSE
2325  nrowa = k
2326  ncola = m
2327  END IF
2328  IF( notb )THEN
2329  nrowb = k
2330  ELSE
2331  nrowb = n
2332  END IF
2333 *
2334 * Test the input parameters.
2335 *
2336  info = 0
2337  IF( ( .NOT.nota ).AND.
2338  $ ( .NOT.lsame( transa, 'C' ) ).AND.
2339  $ ( .NOT.lsame( transa, 'T' ) ) )THEN
2340  info = 1
2341  ELSE IF( ( .NOT.notb ).AND.
2342  $ ( .NOT.lsame( transb, 'C' ) ).AND.
2343  $ ( .NOT.lsame( transb, 'T' ) ) )THEN
2344  info = 2
2345  ELSE IF( m .LT.0 )THEN
2346  info = 3
2347  ELSE IF( n .LT.0 )THEN
2348  info = 4
2349  ELSE IF( k .LT.0 )THEN
2350  info = 5
2351  ELSE IF( lda.LT.max( 1, nrowa ) )THEN
2352  info = 8
2353  ELSE IF( ldb.LT.max( 1, nrowb ) )THEN
2354  info = 10
2355  ELSE IF( ldc.LT.max( 1, m ) )THEN
2356  info = 13
2357  END IF
2358  IF( info.NE.0 )THEN
2359  CALL xerbla( 'DGEMM ', info )
2360  RETURN
2361  END IF
2362 *
2363 * Quick return if possible.
2364 *
2365  IF( ( m.EQ.0 ).OR.( n.EQ.0 ).OR.
2366  $ ( ( ( alpha.EQ.zero ).OR.( k.EQ.0 ) ).AND.( beta.EQ.one ) ) )
2367  $ RETURN
2368 *
2369 * And if alpha.eq.zero.
2370 *
2371  IF( alpha.EQ.zero )THEN
2372  IF( beta.EQ.zero )THEN
2373  DO 20, j = 1, n
2374  DO 10, i = 1, m
2375  c( i, j ) = zero
2376  10 CONTINUE
2377  20 CONTINUE
2378  ELSE
2379  DO 40, j = 1, n
2380  DO 30, i = 1, m
2381  c( i, j ) = beta*c( i, j )
2382  30 CONTINUE
2383  40 CONTINUE
2384  END IF
2385  RETURN
2386  END IF
2387 *
2388 * Start the operations.
2389 *
2390  IF( notb )THEN
2391  IF( nota )THEN
2392 *
2393 * Form C := alpha*A*B + beta*C.
2394 *
2395  DO 90, j = 1, n
2396  IF( beta.EQ.zero )THEN
2397  DO 50, i = 1, m
2398  c( i, j ) = zero
2399  50 CONTINUE
2400  ELSE IF( beta.NE.one )THEN
2401  DO 60, i = 1, m
2402  c( i, j ) = beta*c( i, j )
2403  60 CONTINUE
2404  END IF
2405  DO 80, l = 1, k
2406  IF( b( l, j ).NE.zero )THEN
2407  temp = alpha*b( l, j )
2408  DO 70, i = 1, m
2409  c( i, j ) = c( i, j ) + temp*a( i, l )
2410  70 CONTINUE
2411  END IF
2412  80 CONTINUE
2413  90 CONTINUE
2414  ELSE
2415 *
2416 * Form C := alpha*A'*B + beta*C
2417 *
2418  DO 120, j = 1, n
2419  DO 110, i = 1, m
2420  temp = zero
2421  DO 100, l = 1, k
2422  temp = temp + a( l, i )*b( l, j )
2423  100 CONTINUE
2424  IF( beta.EQ.zero )THEN
2425  c( i, j ) = alpha*temp
2426  ELSE
2427  c( i, j ) = alpha*temp + beta*c( i, j )
2428  END IF
2429  110 CONTINUE
2430  120 CONTINUE
2431  END IF
2432  ELSE
2433  IF( nota )THEN
2434 *
2435 * Form C := alpha*A*B' + beta*C
2436 *
2437  DO 170, j = 1, n
2438  IF( beta.EQ.zero )THEN
2439  DO 130, i = 1, m
2440  c( i, j ) = zero
2441  130 CONTINUE
2442  ELSE IF( beta.NE.one )THEN
2443  DO 140, i = 1, m
2444  c( i, j ) = beta*c( i, j )
2445  140 CONTINUE
2446  END IF
2447  DO 160, l = 1, k
2448  IF( b( j, l ).NE.zero )THEN
2449  temp = alpha*b( j, l )
2450  DO 150, i = 1, m
2451  c( i, j ) = c( i, j ) + temp*a( i, l )
2452  150 CONTINUE
2453  END IF
2454  160 CONTINUE
2455  170 CONTINUE
2456  ELSE
2457 *
2458 * Form C := alpha*A'*B' + beta*C
2459 *
2460  DO 200, j = 1, n
2461  DO 190, i = 1, m
2462  temp = zero
2463  DO 180, l = 1, k
2464  temp = temp + a( l, i )*b( j, l )
2465  180 CONTINUE
2466  IF( beta.EQ.zero )THEN
2467  c( i, j ) = alpha*temp
2468  ELSE
2469  c( i, j ) = alpha*temp + beta*c( i, j )
2470  END IF
2471  190 CONTINUE
2472  200 CONTINUE
2473  END IF
2474  END IF
2475 *
2476  RETURN
2477 *
2478 * End of DGEMM .
2479 *
max
#define max(a, b)
Definition: cascade.c:32

◆ dger()

subroutine dger ( integer  M,
integer  N,
double precision  ALPHA,
double precision, dimension( * )  X,
integer  INCX,
double precision, dimension( * )  Y,
integer  INCY,
double precision, dimension( lda, * )  A,
integer  LDA 
)
1488 * .. Scalar Arguments ..
1489  DOUBLE PRECISION alpha
1490  INTEGER incx, incy, lda, m, n
1491 * .. Array Arguments ..
1492  DOUBLE PRECISION a( lda, * ), x( * ), y( * )
1493 * ..
1494 *
1495 * Purpose
1496 * =======
1497 *
1498 * DGER performs the rank 1 operation
1499 *
1500 * A := alpha*x*y' + A,
1501 *
1502 * where alpha is a scalar, x is an m element vector, y is an n element
1503 * vector and A is an m by n matrix.
1504 *
1505 * Parameters
1506 * ==========
1507 *
1508 * M - INTEGER.
1509 * On entry, M specifies the number of rows of the matrix A.
1510 * M must be at least zero.
1511 * Unchanged on exit.
1512 *
1513 * N - INTEGER.
1514 * On entry, N specifies the number of columns of the matrix A.
1515 * N must be at least zero.
1516 * Unchanged on exit.
1517 *
1518 * ALPHA - DOUBLE PRECISION.
1519 * On entry, ALPHA specifies the scalar alpha.
1520 * Unchanged on exit.
1521 *
1522 * X - DOUBLE PRECISION array of dimension at least
1523 * ( 1 + ( m - 1 )*abs( INCX ) ).
1524 * Before entry, the incremented array X must contain the m
1525 * element vector x.
1526 * Unchanged on exit.
1527 *
1528 * INCX - INTEGER.
1529 * On entry, INCX specifies the increment for the elements of
1530 * X. INCX must not be zero.
1531 * Unchanged on exit.
1532 *
1533 * Y - DOUBLE PRECISION array of dimension at least
1534 * ( 1 + ( n - 1 )*abs( INCY ) ).
1535 * Before entry, the incremented array Y must contain the n
1536 * element vector y.
1537 * Unchanged on exit.
1538 *
1539 * INCY - INTEGER.
1540 * On entry, INCY specifies the increment for the elements of
1541 * Y. INCY must not be zero.
1542 * Unchanged on exit.
1543 *
1544 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
1545 * Before entry, the leading m by n part of the array A must
1546 * contain the matrix of coefficients. On exit, A is
1547 * overwritten by the updated matrix.
1548 *
1549 * LDA - INTEGER.
1550 * On entry, LDA specifies the first dimension of A as declared
1551 * in the calling (sub) program. LDA must be at least
1552 * max( 1, m ).
1553 * Unchanged on exit.
1554 *
1555 *
1556 * Level 2 Blas routine.
1557 *
1558 * -- Written on 22-October-1986.
1559 * Jack Dongarra, Argonne National Lab.
1560 * Jeremy Du Croz, Nag Central Office.
1561 * Sven Hammarling, Nag Central Office.
1562 * Richard Hanson, Sandia National Labs.
1563 *
1564 *
1565 * .. Parameters ..
1566  DOUBLE PRECISION zero
1567  parameter( zero = 0.0d+0 )
1568 * .. Local Scalars ..
1569  DOUBLE PRECISION temp
1570  INTEGER i, info, ix, j, jy, kx
1571 * .. External Subroutines ..
1572  EXTERNAL xerbla
1573 * .. Intrinsic Functions ..
1574  INTRINSIC max
1575 * ..
1576 * .. Executable Statements ..
1577 *
1578 * Test the input parameters.
1579 *
1580  info = 0
1581  IF ( m.LT.0 )THEN
1582  info = 1
1583  ELSE IF( n.LT.0 )THEN
1584  info = 2
1585  ELSE IF( incx.EQ.0 )THEN
1586  info = 5
1587  ELSE IF( incy.EQ.0 )THEN
1588  info = 7
1589  ELSE IF( lda.LT.max( 1, m ) )THEN
1590  info = 9
1591  END IF
1592  IF( info.NE.0 )THEN
1593  CALL xerbla( 'DGER ', info )
1594  RETURN
1595  END IF
1596 *
1597 * Quick return if possible.
1598 *
1599  IF( ( m.EQ.0 ).OR.( n.EQ.0 ).OR.( alpha.EQ.zero ) )
1600  $ RETURN
1601 *
1602 * Start the operations. In this version the elements of A are
1603 * accessed sequentially with one pass through A.
1604 *
1605  IF( incy.GT.0 )THEN
1606  jy = 1
1607  ELSE
1608  jy = 1 - ( n - 1 )*incy
1609  END IF
1610  IF( incx.EQ.1 )THEN
1611  DO 20, j = 1, n
1612  IF( y( jy ).NE.zero )THEN
1613  temp = alpha*y( jy )
1614  DO 10, i = 1, m
1615  a( i, j ) = a( i, j ) + x( i )*temp
1616  10 CONTINUE
1617  END IF
1618  jy = jy + incy
1619  20 CONTINUE
1620  ELSE
1621  IF( incx.GT.0 )THEN
1622  kx = 1
1623  ELSE
1624  kx = 1 - ( m - 1 )*incx
1625  END IF
1626  DO 40, j = 1, n
1627  IF( y( jy ).NE.zero )THEN
1628  temp = alpha*y( jy )
1629  ix = kx
1630  DO 30, i = 1, m
1631  a( i, j ) = a( i, j ) + x( ix )*temp
1632  ix = ix + incx
1633  30 CONTINUE
1634  END IF
1635  jy = jy + incy
1636  40 CONTINUE
1637  END IF
1638 *
1639  RETURN
1640 *
1641 * End of DGER .
1642 *
max
#define max(a, b)
Definition: cascade.c:32

◆ dgesv()

subroutine dgesv ( integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)
58 *
59 * -- LAPACK driver routine (version 3.0) --
60 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
61 * Courant Institute, Argonne National Lab, and Rice University
62 * March 31, 1993
63 *
64 * .. Scalar Arguments ..
65  INTEGER info, lda, ldb, n, nrhs
66 * ..
67 * .. Array Arguments ..
68  INTEGER ipiv( * )
69  DOUBLE PRECISION a( lda, * ), b( ldb, * )
70 * ..
71 *
72 * Purpose
73 * =======
74 *
75 * DGESV computes the solution to a real system of linear equations
76 * A * X = B,
77 * where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
78 *
79 * The LU decomposition with partial pivoting and row interchanges is
80 * used to factor A as
81 * A = P * L * U,
82 * where P is a permutation matrix, L is unit lower triangular, and U is
83 * upper triangular. The factored form of A is then used to solve the
84 * system of equations A * X = B.
85 *
86 * Arguments
87 * =========
88 *
89 * N (input) INTEGER
90 * The number of linear equations, i.e., the order of the
91 * matrix A. N >= 0.
92 *
93 * NRHS (input) INTEGER
94 * The number of right hand sides, i.e., the number of columns
95 * of the matrix B. NRHS >= 0.
96 *
97 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
98 * On entry, the N-by-N coefficient matrix A.
99 * On exit, the factors L and U from the factorization
100 * A = P*L*U; the unit diagonal elements of L are not stored.
101 *
102 * LDA (input) INTEGER
103 * The leading dimension of the array A. LDA >= max(1,N).
104 *
105 * IPIV (output) INTEGER array, dimension (N)
106 * The pivot indices that define the permutation matrix P;
107 * row i of the matrix was interchanged with row IPIV(i).
108 *
109 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
110 * On entry, the N-by-NRHS matrix of right hand side matrix B.
111 * On exit, if INFO = 0, the N-by-NRHS solution matrix X.
112 *
113 * LDB (input) INTEGER
114 * The leading dimension of the array B. LDB >= max(1,N).
115 *
116 * INFO (output) INTEGER
117 * = 0: successful exit
118 * < 0: if INFO = -i, the i-th argument had an illegal value
119 * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
120 * has been completed, but the factor U is exactly
121 * singular, so the solution could not be computed.
122 *
123 * =====================================================================
124 *
125 * .. External Subroutines ..
126  EXTERNAL dgetrf, dgetrs, xerbla
127 * ..
128 * .. Intrinsic Functions ..
129  INTRINSIC max
130 * ..
131 * .. Executable Statements ..
132 *
133 * Test the input parameters.
134 *
135  info = 0
136  IF( n.LT.0 ) THEN
137  info = -1
138  ELSE IF( nrhs.LT.0 ) THEN
139  info = -2
140  ELSE IF( lda.LT.max( 1, n ) ) THEN
141  info = -4
142  ELSE IF( ldb.LT.max( 1, n ) ) THEN
143  info = -7
144  END IF
145  IF( info.NE.0 ) THEN
146  CALL xerbla( 'DGESV ', -info )
147  RETURN
148  END IF
149 *
150 * Compute the LU factorization of A.
151 *
152  CALL dgetrf( n, n, a, lda, ipiv, info )
153  IF( info.EQ.0 ) THEN
154 *
155 * Solve the system A*X = B, overwriting B with X.
156 *
157  CALL dgetrs( 'No transpose', n, nrhs, a, lda, ipiv, b, ldb,
158  $ info )
159  END IF
160  RETURN
161 *
162 * End of DGESV
163 *
max
#define max(a, b)
Definition: cascade.c:32
dgetrs
subroutine dgetrs(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO)
Definition: dgesv.f:461
dgetrf
subroutine dgetrf(M, N, A, LDA, IPIV, INFO)
Definition: dgesv.f:301

◆ dgetf2()

subroutine dgetf2 ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)
166 *
167 * -- LAPACK routine (version 3.0) --
168 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
169 * Courant Institute, Argonne National Lab, and Rice University
170 * June 30, 1992
171 *
172 * .. Scalar Arguments ..
173  INTEGER info, lda, m, n
174 * ..
175 * .. Array Arguments ..
176  INTEGER ipiv( * )
177  DOUBLE PRECISION a( lda, * )
178 * ..
179 *
180 * Purpose
181 * =======
182 *
183 * DGETF2 computes an LU factorization of a general m-by-n matrix A
184 * using partial pivoting with row interchanges.
185 *
186 * The factorization has the form
187 * A = P * L * U
188 * where P is a permutation matrix, L is lower triangular with unit
189 * diagonal elements (lower trapezoidal if m > n), and U is upper
190 * triangular (upper trapezoidal if m < n).
191 *
192 * This is the right-looking Level 2 BLAS version of the algorithm.
193 *
194 * Arguments
195 * =========
196 *
197 * M (input) INTEGER
198 * The number of rows of the matrix A. M >= 0.
199 *
200 * N (input) INTEGER
201 * The number of columns of the matrix A. N >= 0.
202 *
203 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
204 * On entry, the m by n matrix to be factored.
205 * On exit, the factors L and U from the factorization
206 * A = P*L*U; the unit diagonal elements of L are not stored.
207 *
208 * LDA (input) INTEGER
209 * The leading dimension of the array A. LDA >= max(1,M).
210 *
211 * IPIV (output) INTEGER array, dimension (min(M,N))
212 * The pivot indices; for 1 <= i <= min(M,N), row i of the
213 * matrix was interchanged with row IPIV(i).
214 *
215 * INFO (output) INTEGER
216 * = 0: successful exit
217 * < 0: if INFO = -k, the k-th argument had an illegal value
218 * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
219 * has been completed, but the factor U is exactly
220 * singular, and division by zero will occur if it is used
221 * to solve a system of equations.
222 *
223 * =====================================================================
224 *
225 * .. Parameters ..
226  DOUBLE PRECISION one, zero
227  parameter( one = 1.0d+0, zero = 0.0d+0 )
228 * ..
229 * .. Local Scalars ..
230  INTEGER j, jp
231 * ..
232 * .. External Functions ..
233  INTEGER idamax
234  EXTERNAL idamax
235 * ..
236 * .. External Subroutines ..
237  EXTERNAL dger, dscal, dswap, xerbla
238 * ..
239 * .. Intrinsic Functions ..
240  INTRINSIC max, min
241 * ..
242 * .. Executable Statements ..
243 *
244 * Test the input parameters.
245 *
246  info = 0
247  IF( m.LT.0 ) THEN
248  info = -1
249  ELSE IF( n.LT.0 ) THEN
250  info = -2
251  ELSE IF( lda.LT.max( 1, m ) ) THEN
252  info = -4
253  END IF
254  IF( info.NE.0 ) THEN
255  CALL xerbla( 'DGETF2', -info )
256  RETURN
257  END IF
258 *
259 * Quick return if possible
260 *
261  IF( m.EQ.0 .OR. n.EQ.0 )
262  $ RETURN
263 *
264  DO 10 j = 1, min( m, n )
265 *
266 * Find pivot and test for singularity.
267 *
268  jp = j - 1 + idamax( m-j+1, a( j, j ), 1 )
269  ipiv( j ) = jp
270  IF( a( jp, j ).NE.zero ) THEN
271 *
272 * Apply the interchange to columns 1:N.
273 *
274  IF( jp.NE.j )
275  $ CALL dswap( n, a( j, 1 ), lda, a( jp, 1 ), lda )
276 *
277 * Compute elements J+1:M of J-th column.
278 *
279  IF( j.LT.m )
280  $ CALL dscal( m-j, one / a( j, j ), a( j+1, j ), 1 )
281 *
282  ELSE IF( info.EQ.0 ) THEN
283 *
284  info = j
285  END IF
286 *
287  IF( j.LT.min( m, n ) ) THEN
288 *
289 * Update trailing submatrix.
290 *
291  CALL dger( m-j, n-j, -one, a( j+1, j ), 1, a( j, j+1 ), lda,
292  $ a( j+1, j+1 ), lda )
293  END IF
294  10 CONTINUE
295  RETURN
296 *
297 * End of DGETF2
298 *
max
#define max(a, b)
Definition: cascade.c:32
min
#define min(a, b)
Definition: cascade.c:31
dscal
subroutine dscal(n, da, dx, incx)
Definition: dgesv.f:1746
idamax
integer function idamax(n, dx, incx)
Definition: dgesv.f:1448
dger
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
Definition: dgesv.f:1488
dswap
subroutine dswap(n, dx, incx, dy, incy)
Definition: dgesv.f:1689

◆ dgetrf()

subroutine dgetrf ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer  INFO 
)
301 *
302 * -- LAPACK routine (version 3.0) --
303 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
304 * Courant Institute, Argonne National Lab, and Rice University
305 * March 31, 1993
306 *
307 * .. Scalar Arguments ..
308  INTEGER info, lda, m, n
309 * ..
310 * .. Array Arguments ..
311  INTEGER ipiv( * )
312  DOUBLE PRECISION a( lda, * )
313 * ..
314 *
315 * Purpose
316 * =======
317 *
318 * DGETRF computes an LU factorization of a general M-by-N matrix A
319 * using partial pivoting with row interchanges.
320 *
321 * The factorization has the form
322 * A = P * L * U
323 * where P is a permutation matrix, L is lower triangular with unit
324 * diagonal elements (lower trapezoidal if m > n), and U is upper
325 * triangular (upper trapezoidal if m < n).
326 *
327 * This is the right-looking Level 3 BLAS version of the algorithm.
328 *
329 * Arguments
330 * =========
331 *
332 * M (input) INTEGER
333 * The number of rows of the matrix A. M >= 0.
334 *
335 * N (input) INTEGER
336 * The number of columns of the matrix A. N >= 0.
337 *
338 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
339 * On entry, the M-by-N matrix to be factored.
340 * On exit, the factors L and U from the factorization
341 * A = P*L*U; the unit diagonal elements of L are not stored.
342 *
343 * LDA (input) INTEGER
344 * The leading dimension of the array A. LDA >= max(1,M).
345 *
346 * IPIV (output) INTEGER array, dimension (min(M,N))
347 * The pivot indices; for 1 <= i <= min(M,N), row i of the
348 * matrix was interchanged with row IPIV(i).
349 *
350 * INFO (output) INTEGER
351 * = 0: successful exit
352 * < 0: if INFO = -i, the i-th argument had an illegal value
353 * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
354 * has been completed, but the factor U is exactly
355 * singular, and division by zero will occur if it is used
356 * to solve a system of equations.
357 *
358 * =====================================================================
359 *
360 * .. Parameters ..
361  DOUBLE PRECISION one
362  parameter( one = 1.0d+0 )
363 * ..
364 * .. Local Scalars ..
365  INTEGER i, iinfo, j, jb, nb
366 * ..
367 * .. External Subroutines ..
368  EXTERNAL dgemm, dgetf2, dlaswp, dtrsm, xerbla
369 * ..
370 * .. External Functions ..
371  INTEGER ilaenv
372  EXTERNAL ilaenv
373 * ..
374 * .. Intrinsic Functions ..
375  INTRINSIC max, min
376 * ..
377 * .. Executable Statements ..
378 *
379 * Test the input parameters.
380 *
381  info = 0
382  IF( m.LT.0 ) THEN
383  info = -1
384  ELSE IF( n.LT.0 ) THEN
385  info = -2
386  ELSE IF( lda.LT.max( 1, m ) ) THEN
387  info = -4
388  END IF
389  IF( info.NE.0 ) THEN
390  CALL xerbla( 'DGETRF', -info )
391  RETURN
392  END IF
393 *
394 * Quick return if possible
395 *
396  IF( m.EQ.0 .OR. n.EQ.0 )
397  $ RETURN
398 *
399 * Determine the block size for this environment.
400 *
401  nb = ilaenv( 1, 'DGETRF', ' ', m, n, -1, -1 )
402  IF( nb.LE.1 .OR. nb.GE.min( m, n ) ) THEN
403 *
404 * Use unblocked code.
405 *
406  CALL dgetf2( m, n, a, lda, ipiv, info )
407  ELSE
408 *
409 * Use blocked code.
410 *
411  DO 20 j = 1, min( m, n ), nb
412  jb = min( min( m, n )-j+1, nb )
413 *
414 * Factor diagonal and subdiagonal blocks and test for exact
415 * singularity.
416 *
417  CALL dgetf2( m-j+1, jb, a( j, j ), lda, ipiv( j ), iinfo )
418 *
419 * Adjust INFO and the pivot indices.
420 *
421  IF( info.EQ.0 .AND. iinfo.GT.0 )
422  $ info = iinfo + j - 1
423  DO 10 i = j, min( m, j+jb-1 )
424  ipiv( i ) = j - 1 + ipiv( i )
425  10 CONTINUE
426 *
427 * Apply interchanges to columns 1:J-1.
428 *
429  CALL dlaswp( j-1, a, lda, j, j+jb-1, ipiv, 1 )
430 *
431  IF( j+jb.LE.n ) THEN
432 *
433 * Apply interchanges to columns J+JB:N.
434 *
435  CALL dlaswp( n-j-jb+1, a( 1, j+jb ), lda, j, j+jb-1,
436  $ ipiv, 1 )
437 *
438 * Compute block row of U.
439 *
440  CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit', jb,
441  $ n-j-jb+1, one, a( j, j ), lda, a( j, j+jb ),
442  $ lda )
443  IF( j+jb.LE.m ) THEN
444 *
445 * Update trailing submatrix.
446 *
447  CALL dgemm( 'No transpose', 'No transpose', m-j-jb+1,
448  $ n-j-jb+1, jb, -one, a( j+jb, j ), lda,
449  $ a( j, j+jb ), lda, one, a( j+jb, j+jb ),
450  $ lda )
451  END IF
452  END IF
453  20 CONTINUE
454  END IF
455  RETURN
456 *
457 * End of DGETRF
458 *
max
#define max(a, b)
Definition: cascade.c:32
min
#define min(a, b)
Definition: cascade.c:31
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dgesv.f:2170
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: dgesv.f:880
dtrsm
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dgesv.f:1791
dlaswp
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
Definition: dgesv.f:611
dgetf2
subroutine dgetf2(M, N, A, LDA, IPIV, INFO)
Definition: dgesv.f:166

◆ dgetrs()

subroutine dgetrs ( character  TRANS,
integer  N,
integer  NRHS,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  INFO 
)
461 *
462 * -- LAPACK routine (version 3.0) --
463 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
464 * Courant Institute, Argonne National Lab, and Rice University
465 * March 31, 1993
466 *
467 * .. Scalar Arguments ..
468  CHARACTER trans
469  INTEGER info, lda, ldb, n, nrhs
470 * ..
471 * .. Array Arguments ..
472  INTEGER ipiv( * )
473  DOUBLE PRECISION a( lda, * ), b( ldb, * )
474 * ..
475 *
476 * Purpose
477 * =======
478 *
479 * DGETRS solves a system of linear equations
480 * A * X = B or A' * X = B
481 * with a general N-by-N matrix A using the LU factorization computed
482 * by DGETRF.
483 *
484 * Arguments
485 * =========
486 *
487 * TRANS (input) CHARACTER*1
488 * Specifies the form of the system of equations:
489 * = 'N': A * X = B (No transpose)
490 * = 'T': A'* X = B (Transpose)
491 * = 'C': A'* X = B (Conjugate transpose = Transpose)
492 *
493 * N (input) INTEGER
494 * The order of the matrix A. N >= 0.
495 *
496 * NRHS (input) INTEGER
497 * The number of right hand sides, i.e., the number of columns
498 * of the matrix B. NRHS >= 0.
499 *
500 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
501 * The factors L and U from the factorization A = P*L*U
502 * as computed by DGETRF.
503 *
504 * LDA (input) INTEGER
505 * The leading dimension of the array A. LDA >= max(1,N).
506 *
507 * IPIV (input) INTEGER array, dimension (N)
508 * The pivot indices from DGETRF; for 1<=i<=N, row i of the
509 * matrix was interchanged with row IPIV(i).
510 *
511 * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
512 * On entry, the right hand side matrix B.
513 * On exit, the solution matrix X.
514 *
515 * LDB (input) INTEGER
516 * The leading dimension of the array B. LDB >= max(1,N).
517 *
518 * INFO (output) INTEGER
519 * = 0: successful exit
520 * < 0: if INFO = -i, the i-th argument had an illegal value
521 *
522 * =====================================================================
523 *
524 * .. Parameters ..
525  DOUBLE PRECISION one
526  parameter( one = 1.0d+0 )
527 * ..
528 * .. Local Scalars ..
529  LOGICAL notran
530 * ..
531 * .. External Functions ..
532  LOGICAL lsame
533  EXTERNAL lsame
534 * ..
535 * .. External Subroutines ..
536  EXTERNAL dlaswp, dtrsm, xerbla
537 * ..
538 * .. Intrinsic Functions ..
539  INTRINSIC max
540 * ..
541 * .. Executable Statements ..
542 *
543 * Test the input parameters.
544 *
545  info = 0
546  notran = lsame( trans, 'N' )
547  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) .AND. .NOT.
548  $ lsame( trans, 'C' ) ) THEN
549  info = -1
550  ELSE IF( n.LT.0 ) THEN
551  info = -2
552  ELSE IF( nrhs.LT.0 ) THEN
553  info = -3
554  ELSE IF( lda.LT.max( 1, n ) ) THEN
555  info = -5
556  ELSE IF( ldb.LT.max( 1, n ) ) THEN
557  info = -8
558  END IF
559  IF( info.NE.0 ) THEN
560  CALL xerbla( 'DGETRS', -info )
561  RETURN
562  END IF
563 *
564 * Quick return if possible
565 *
566  IF( n.EQ.0 .OR. nrhs.EQ.0 )
567  $ RETURN
568 *
569  IF( notran ) THEN
570 *
571 * Solve A * X = B.
572 *
573 * Apply row interchanges to the right hand sides.
574 *
575  CALL dlaswp( nrhs, b, ldb, 1, n, ipiv, 1 )
576 *
577 * Solve L*X = B, overwriting B with X.
578 *
579  CALL dtrsm( 'Left', 'Lower', 'No transpose', 'Unit', n, nrhs,
580  $ one, a, lda, b, ldb )
581 *
582 * Solve U*X = B, overwriting B with X.
583 *
584  CALL dtrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', n,
585  $ nrhs, one, a, lda, b, ldb )
586  ELSE
587 *
588 * Solve A' * X = B.
589 *
590 * Solve U'*X = B, overwriting B with X.
591 *
592  CALL dtrsm( 'Left', 'Upper', 'Transpose', 'Non-unit', n, nrhs,
593  $ one, a, lda, b, ldb )
594 *
595 * Solve L'*X = B, overwriting B with X.
596 *
597  CALL dtrsm( 'Left', 'Lower', 'Transpose', 'Unit', n, nrhs, one,
598  $ a, lda, b, ldb )
599 *
600 * Apply row interchanges to the solution vectors.
601 *
602  CALL dlaswp( nrhs, b, ldb, 1, n, ipiv, -1 )
603  END IF
604 *
605  RETURN
606 *
607 * End of DGETRS
608 *
max
#define max(a, b)
Definition: cascade.c:32
dtrsm
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dgesv.f:1791
dlaswp
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
Definition: dgesv.f:611

◆ dlaswp()

subroutine dlaswp ( integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer  K1,
integer  K2,
integer, dimension( * )  IPIV,
integer  INCX 
)
611 *
612 * -- LAPACK auxiliary routine (version 3.0) --
613 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
614 * Courant Institute, Argonne National Lab, and Rice University
615 * June 30, 1999
616 *
617 * .. Scalar Arguments ..
618  INTEGER incx, k1, k2, lda, n
619 * ..
620 * .. Array Arguments ..
621  INTEGER ipiv( * )
622  DOUBLE PRECISION a( lda, * )
623 * ..
624 *
625 * Purpose
626 * =======
627 *
628 * DLASWP performs a series of row interchanges on the matrix A.
629 * One row interchange is initiated for each of rows K1 through K2 of A.
630 *
631 * Arguments
632 * =========
633 *
634 * N (input) INTEGER
635 * The number of columns of the matrix A.
636 *
637 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
638 * On entry, the matrix of column dimension N to which the row
639 * interchanges will be applied.
640 * On exit, the permuted matrix.
641 *
642 * LDA (input) INTEGER
643 * The leading dimension of the array A.
644 *
645 * K1 (input) INTEGER
646 * The first element of IPIV for which a row interchange will
647 * be done.
648 *
649 * K2 (input) INTEGER
650 * The last element of IPIV for which a row interchange will
651 * be done.
652 *
653 * IPIV (input) INTEGER array, dimension (M*abs(INCX))
654 * The vector of pivot indices. Only the elements in positions
655 * K1 through K2 of IPIV are accessed.
656 * IPIV(K) = L implies rows K and L are to be interchanged.
657 *
658 * INCX (input) INTEGER
659 * The increment between successive values of IPIV. If IPIV
660 * is negative, the pivots are applied in reverse order.
661 *
662 * Further Details
663 * ===============
664 *
665 * Modified by
666 * R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA
667 *
668 * =====================================================================
669 *
670 * .. Local Scalars ..
671  INTEGER i, i1, i2, inc, ip, ix, ix0, j, k, n32
672  DOUBLE PRECISION temp
673 * ..
674 * .. Executable Statements ..
675 *
676 * Interchange row I with row IPIV(I) for each of rows K1 through K2.
677 *
678  IF( incx.GT.0 ) THEN
679  ix0 = k1
680  i1 = k1
681  i2 = k2
682  inc = 1
683  ELSE IF( incx.LT.0 ) THEN
684  ix0 = 1 + ( 1-k2 )*incx
685  i1 = k2
686  i2 = k1
687  inc = -1
688  ELSE
689  RETURN
690  END IF
691 *
692  n32 = ( n / 32 )*32
693  IF( n32.NE.0 ) THEN
694  DO 30 j = 1, n32, 32
695  ix = ix0
696  DO 20 i = i1, i2, inc
697  ip = ipiv( ix )
698  IF( ip.NE.i ) THEN
699  DO 10 k = j, j + 31
700  temp = a( i, k )
701  a( i, k ) = a( ip, k )
702  a( ip, k ) = temp
703  10 CONTINUE
704  END IF
705  ix = ix + incx
706  20 CONTINUE
707  30 CONTINUE
708  END IF
709  IF( n32.NE.n ) THEN
710  n32 = n32 + 1
711  ix = ix0
712  DO 50 i = i1, i2, inc
713  ip = ipiv( ix )
714  IF( ip.NE.i ) THEN
715  DO 40 k = n32, n
716  temp = a( i, k )
717  a( i, k ) = a( ip, k )
718  a( ip, k ) = temp
719  40 CONTINUE
720  END IF
721  ix = ix + incx
722  50 CONTINUE
723  END IF
724 *
725  RETURN
726 *
727 * End of DLASWP
728 *

◆ dscal()

subroutine dscal ( integer  n,
double precision  da,
double precision, dimension(*)  dx,
integer  incx 
)
1746 c
1747 c scales a vector by a constant.
1748 c uses unrolled loops for increment equal to one.
1749 c jack dongarra, linpack, 3/11/78.
1750 c modified 3/93 to return if incx .le. 0.
1751 c modified 12/3/93, array(1) declarations changed to array(*)
1752 c
1753  double precision da,dx(*)
1754  integer i,incx,m,mp1,n,nincx
1755 c
1756  if( n.le.0 .or. incx.le.0 )return
1757  if(incx.eq.1)go to 20
1758 c
1759 c code for increment not equal to 1
1760 c
1761  nincx = n*incx
1762  do 10 i = 1,nincx,incx
1763  dx(i) = da*dx(i)
1764  10 continue
1765  return
1766 c
1767 c code for increment equal to 1
1768 c
1769 c
1770 c clean-up loop
1771 c
1772  20 m = mod(n,5)
1773  if( m .eq. 0 ) go to 40
1774  do 30 i = 1,m
1775  dx(i) = da*dx(i)
1776  30 continue
1777  if( n .lt. 5 ) return
1778  40 mp1 = m + 1
1779  do 50 i = mp1,n,5
1780  dx(i) = da*dx(i)
1781  dx(i + 1) = da*dx(i + 1)
1782  dx(i + 2) = da*dx(i + 2)
1783  dx(i + 3) = da*dx(i + 3)
1784  dx(i + 4) = da*dx(i + 4)
1785  50 continue
1786  return

◆ dswap()

subroutine dswap ( integer  n,
double precision, dimension(*)  dx,
integer  incx,
double precision, dimension(*)  dy,
integer  incy 
)
1689 c
1690 c interchanges two vectors.
1691 c uses unrolled loops for increments equal one.
1692 c jack dongarra, linpack, 3/11/78.
1693 c modified 12/3/93, array(1) declarations changed to array(*)
1694 c
1695  double precision dx(*),dy(*),dtemp
1696  integer i,incx,incy,ix,iy,m,mp1,n
1697 c
1698  if(n.le.0)return
1699  if(incx.eq.1.and.incy.eq.1)go to 20
1700 c
1701 c code for unequal increments or equal increments not equal
1702 c to 1
1703 c
1704  ix = 1
1705  iy = 1
1706  if(incx.lt.0)ix = (-n+1)*incx + 1
1707  if(incy.lt.0)iy = (-n+1)*incy + 1
1708  do 10 i = 1,n
1709  dtemp = dx(ix)
1710  dx(ix) = dy(iy)
1711  dy(iy) = dtemp
1712  ix = ix + incx
1713  iy = iy + incy
1714  10 continue
1715  return
1716 c
1717 c code for both increments equal to 1
1718 c
1719 c
1720 c clean-up loop
1721 c
1722  20 m = mod(n,3)
1723  if( m .eq. 0 ) go to 40
1724  do 30 i = 1,m
1725  dtemp = dx(i)
1726  dx(i) = dy(i)
1727  dy(i) = dtemp
1728  30 continue
1729  if( n .lt. 3 ) return
1730  40 mp1 = m + 1
1731  do 50 i = mp1,n,3
1732  dtemp = dx(i)
1733  dx(i) = dy(i)
1734  dy(i) = dtemp
1735  dtemp = dx(i + 1)
1736  dx(i + 1) = dy(i + 1)
1737  dy(i + 1) = dtemp
1738  dtemp = dx(i + 2)
1739  dx(i + 2) = dy(i + 2)
1740  dy(i + 2) = dtemp
1741  50 continue
1742  return

◆ dtrsm()

subroutine dtrsm ( character*1  SIDE,
character*1  UPLO,
character*1  TRANSA,
character*1  DIAG,
integer  M,
integer  N,
double precision  ALPHA,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB 
)
1791 * .. Scalar Arguments ..
1792  CHARACTER*1 side, uplo, transa, diag
1793  INTEGER m, n, lda, ldb
1794  DOUBLE PRECISION alpha
1795 * .. Array Arguments ..
1796  DOUBLE PRECISION a( lda, * ), b( ldb, * )
1797 * ..
1798 *
1799 * Purpose
1800 * =======
1801 *
1802 * DTRSM solves one of the matrix equations
1803 *
1804 * op( A )*X = alpha*B, or X*op( A ) = alpha*B,
1805 *
1806 * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
1807 * non-unit, upper or lower triangular matrix and op( A ) is one of
1808 *
1809 * op( A ) = A or op( A ) = A'.
1810 *
1811 * The matrix X is overwritten on B.
1812 *
1813 * Parameters
1814 * ==========
1815 *
1816 * SIDE - CHARACTER*1.
1817 * On entry, SIDE specifies whether op( A ) appears on the left
1818 * or right of X as follows:
1819 *
1820 * SIDE = 'L' or 'l' op( A )*X = alpha*B.
1821 *
1822 * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
1823 *
1824 * Unchanged on exit.
1825 *
1826 * UPLO - CHARACTER*1.
1827 * On entry, UPLO specifies whether the matrix A is an upper or
1828 * lower triangular matrix as follows:
1829 *
1830 * UPLO = 'U' or 'u' A is an upper triangular matrix.
1831 *
1832 * UPLO = 'L' or 'l' A is a lower triangular matrix.
1833 *
1834 * Unchanged on exit.
1835 *
1836 * TRANSA - CHARACTER*1.
1837 * On entry, TRANSA specifies the form of op( A ) to be used in
1838 * the matrix multiplication as follows:
1839 *
1840 * TRANSA = 'N' or 'n' op( A ) = A.
1841 *
1842 * TRANSA = 'T' or 't' op( A ) = A'.
1843 *
1844 * TRANSA = 'C' or 'c' op( A ) = A'.
1845 *
1846 * Unchanged on exit.
1847 *
1848 * DIAG - CHARACTER*1.
1849 * On entry, DIAG specifies whether or not A is unit triangular
1850 * as follows:
1851 *
1852 * DIAG = 'U' or 'u' A is assumed to be unit triangular.
1853 *
1854 * DIAG = 'N' or 'n' A is not assumed to be unit
1855 * triangular.
1856 *
1857 * Unchanged on exit.
1858 *
1859 * M - INTEGER.
1860 * On entry, M specifies the number of rows of B. M must be at
1861 * least zero.
1862 * Unchanged on exit.
1863 *
1864 * N - INTEGER.
1865 * On entry, N specifies the number of columns of B. N must be
1866 * at least zero.
1867 * Unchanged on exit.
1868 *
1869 * ALPHA - DOUBLE PRECISION.
1870 * On entry, ALPHA specifies the scalar alpha. When alpha is
1871 * zero then A is not referenced and B need not be set before
1872 * entry.
1873 * Unchanged on exit.
1874 *
1875 * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
1876 * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
1877 * Before entry with UPLO = 'U' or 'u', the leading k by k
1878 * upper triangular part of the array A must contain the upper
1879 * triangular matrix and the strictly lower triangular part of
1880 * A is not referenced.
1881 * Before entry with UPLO = 'L' or 'l', the leading k by k
1882 * lower triangular part of the array A must contain the lower
1883 * triangular matrix and the strictly upper triangular part of
1884 * A is not referenced.
1885 * Note that when DIAG = 'U' or 'u', the diagonal elements of
1886 * A are not referenced either, but are assumed to be unity.
1887 * Unchanged on exit.
1888 *
1889 * LDA - INTEGER.
1890 * On entry, LDA specifies the first dimension of A as declared
1891 * in the calling (sub) program. When SIDE = 'L' or 'l' then
1892 * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
1893 * then LDA must be at least max( 1, n ).
1894 * Unchanged on exit.
1895 *
1896 * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
1897 * Before entry, the leading m by n part of the array B must
1898 * contain the right-hand side matrix B, and on exit is
1899 * overwritten by the solution matrix X.
1900 *
1901 * LDB - INTEGER.
1902 * On entry, LDB specifies the first dimension of B as declared
1903 * in the calling (sub) program. LDB must be at least
1904 * max( 1, m ).
1905 * Unchanged on exit.
1906 *
1907 *
1908 * Level 3 Blas routine.
1909 *
1910 *
1911 * -- Written on 8-February-1989.
1912 * Jack Dongarra, Argonne National Laboratory.
1913 * Iain Duff, AERE Harwell.
1914 * Jeremy Du Croz, Numerical Algorithms Group Ltd.
1915 * Sven Hammarling, Numerical Algorithms Group Ltd.
1916 *
1917 *
1918 * .. External Functions ..
1919  LOGICAL lsame
1920  EXTERNAL lsame
1921 * .. External Subroutines ..
1922  EXTERNAL xerbla
1923 * .. Intrinsic Functions ..
1924  INTRINSIC max
1925 * .. Local Scalars ..
1926  LOGICAL lside, nounit, upper
1927  INTEGER i, info, j, k, nrowa
1928  DOUBLE PRECISION temp
1929 * .. Parameters ..
1930  DOUBLE PRECISION one , zero
1931  parameter( one = 1.0d+0, zero = 0.0d+0 )
1932 * ..
1933 * .. Executable Statements ..
1934 *
1935 * Test the input parameters.
1936 *
1937  lside = lsame( side , 'L' )
1938  IF( lside )THEN
1939  nrowa = m
1940  ELSE
1941  nrowa = n
1942  END IF
1943  nounit = lsame( diag , 'N' )
1944  upper = lsame( uplo , 'U' )
1945 *
1946  info = 0
1947  IF( ( .NOT.lside ).AND.
1948  $ ( .NOT.lsame( side , 'R' ) ) )THEN
1949  info = 1
1950  ELSE IF( ( .NOT.upper ).AND.
1951  $ ( .NOT.lsame( uplo , 'L' ) ) )THEN
1952  info = 2
1953  ELSE IF( ( .NOT.lsame( transa, 'N' ) ).AND.
1954  $ ( .NOT.lsame( transa, 'T' ) ).AND.
1955  $ ( .NOT.lsame( transa, 'C' ) ) )THEN
1956  info = 3
1957  ELSE IF( ( .NOT.lsame( diag , 'U' ) ).AND.
1958  $ ( .NOT.lsame( diag , 'N' ) ) )THEN
1959  info = 4
1960  ELSE IF( m .LT.0 )THEN
1961  info = 5
1962  ELSE IF( n .LT.0 )THEN
1963  info = 6
1964  ELSE IF( lda.LT.max( 1, nrowa ) )THEN
1965  info = 9
1966  ELSE IF( ldb.LT.max( 1, m ) )THEN
1967  info = 11
1968  END IF
1969  IF( info.NE.0 )THEN
1970  CALL xerbla( 'DTRSM ', info )
1971  RETURN
1972  END IF
1973 *
1974 * Quick return if possible.
1975 *
1976  IF( n.EQ.0 )
1977  $ RETURN
1978 *
1979 * And when alpha.eq.zero.
1980 *
1981  IF( alpha.EQ.zero )THEN
1982  DO 20, j = 1, n
1983  DO 10, i = 1, m
1984  b( i, j ) = zero
1985  10 CONTINUE
1986  20 CONTINUE
1987  RETURN
1988  END IF
1989 *
1990 * Start the operations.
1991 *
1992  IF( lside )THEN
1993  IF( lsame( transa, 'N' ) )THEN
1994 *
1995 * Form B := alpha*inv( A )*B.
1996 *
1997  IF( upper )THEN
1998  DO 60, j = 1, n
1999  IF( alpha.NE.one )THEN
2000  DO 30, i = 1, m
2001  b( i, j ) = alpha*b( i, j )
2002  30 CONTINUE
2003  END IF
2004  DO 50, k = m, 1, -1
2005  IF( b( k, j ).NE.zero )THEN
2006  IF( nounit )
2007  $ b( k, j ) = b( k, j )/a( k, k )
2008  DO 40, i = 1, k - 1
2009  b( i, j ) = b( i, j ) - b( k, j )*a( i, k )
2010  40 CONTINUE
2011  END IF
2012  50 CONTINUE
2013  60 CONTINUE
2014  ELSE
2015  DO 100, j = 1, n
2016  IF( alpha.NE.one )THEN
2017  DO 70, i = 1, m
2018  b( i, j ) = alpha*b( i, j )
2019  70 CONTINUE
2020  END IF
2021  DO 90 k = 1, m
2022  IF( b( k, j ).NE.zero )THEN
2023  IF( nounit )
2024  $ b( k, j ) = b( k, j )/a( k, k )
2025  DO 80, i = k + 1, m
2026  b( i, j ) = b( i, j ) - b( k, j )*a( i, k )
2027  80 CONTINUE
2028  END IF
2029  90 CONTINUE
2030  100 CONTINUE
2031  END IF
2032  ELSE
2033 *
2034 * Form B := alpha*inv( A' )*B.
2035 *
2036  IF( upper )THEN
2037  DO 130, j = 1, n
2038  DO 120, i = 1, m
2039  temp = alpha*b( i, j )
2040  DO 110, k = 1, i - 1
2041  temp = temp - a( k, i )*b( k, j )
2042  110 CONTINUE
2043  IF( nounit )
2044  $ temp = temp/a( i, i )
2045  b( i, j ) = temp
2046  120 CONTINUE
2047  130 CONTINUE
2048  ELSE
2049  DO 160, j = 1, n
2050  DO 150, i = m, 1, -1
2051  temp = alpha*b( i, j )
2052  DO 140, k = i + 1, m
2053  temp = temp - a( k, i )*b( k, j )
2054  140 CONTINUE
2055  IF( nounit )
2056  $ temp = temp/a( i, i )
2057  b( i, j ) = temp
2058  150 CONTINUE
2059  160 CONTINUE
2060  END IF
2061  END IF
2062  ELSE
2063  IF( lsame( transa, 'N' ) )THEN
2064 *
2065 * Form B := alpha*B*inv( A ).
2066 *
2067  IF( upper )THEN
2068  DO 210, j = 1, n
2069  IF( alpha.NE.one )THEN
2070  DO 170, i = 1, m
2071  b( i, j ) = alpha*b( i, j )
2072  170 CONTINUE
2073  END IF
2074  DO 190, k = 1, j - 1
2075  IF( a( k, j ).NE.zero )THEN
2076  DO 180, i = 1, m
2077  b( i, j ) = b( i, j ) - a( k, j )*b( i, k )
2078  180 CONTINUE
2079  END IF
2080  190 CONTINUE
2081  IF( nounit )THEN
2082  temp = one/a( j, j )
2083  DO 200, i = 1, m
2084  b( i, j ) = temp*b( i, j )
2085  200 CONTINUE
2086  END IF
2087  210 CONTINUE
2088  ELSE
2089  DO 260, j = n, 1, -1
2090  IF( alpha.NE.one )THEN
2091  DO 220, i = 1, m
2092  b( i, j ) = alpha*b( i, j )
2093  220 CONTINUE
2094  END IF
2095  DO 240, k = j + 1, n
2096  IF( a( k, j ).NE.zero )THEN
2097  DO 230, i = 1, m
2098  b( i, j ) = b( i, j ) - a( k, j )*b( i, k )
2099  230 CONTINUE
2100  END IF
2101  240 CONTINUE
2102  IF( nounit )THEN
2103  temp = one/a( j, j )
2104  DO 250, i = 1, m
2105  b( i, j ) = temp*b( i, j )
2106  250 CONTINUE
2107  END IF
2108  260 CONTINUE
2109  END IF
2110  ELSE
2111 *
2112 * Form B := alpha*B*inv( A' ).
2113 *
2114  IF( upper )THEN
2115  DO 310, k = n, 1, -1
2116  IF( nounit )THEN
2117  temp = one/a( k, k )
2118  DO 270, i = 1, m
2119  b( i, k ) = temp*b( i, k )
2120  270 CONTINUE
2121  END IF
2122  DO 290, j = 1, k - 1
2123  IF( a( j, k ).NE.zero )THEN
2124  temp = a( j, k )
2125  DO 280, i = 1, m
2126  b( i, j ) = b( i, j ) - temp*b( i, k )
2127  280 CONTINUE
2128  END IF
2129  290 CONTINUE
2130  IF( alpha.NE.one )THEN
2131  DO 300, i = 1, m
2132  b( i, k ) = alpha*b( i, k )
2133  300 CONTINUE
2134  END IF
2135  310 CONTINUE
2136  ELSE
2137  DO 360, k = 1, n
2138  IF( nounit )THEN
2139  temp = one/a( k, k )
2140  DO 320, i = 1, m
2141  b( i, k ) = temp*b( i, k )
2142  320 CONTINUE
2143  END IF
2144  DO 340, j = k + 1, n
2145  IF( a( j, k ).NE.zero )THEN
2146  temp = a( j, k )
2147  DO 330, i = 1, m
2148  b( i, j ) = b( i, j ) - temp*b( i, k )
2149  330 CONTINUE
2150  END IF
2151  340 CONTINUE
2152  IF( alpha.NE.one )THEN
2153  DO 350, i = 1, m
2154  b( i, k ) = alpha*b( i, k )
2155  350 CONTINUE
2156  END IF
2157  360 CONTINUE
2158  END IF
2159  END IF
2160  END IF
2161 *
2162  RETURN
2163 *
2164 * End of DTRSM .
2165 *
max
#define max(a, b)
Definition: cascade.c:32

◆ idamax()

integer function idamax ( integer  n,
double precision, dimension(*)  dx,
integer  incx 
)
1448 c
1449 c finds the index of element having max. absolute value.
1450 c jack dongarra, linpack, 3/11/78.
1451 c modified 3/93 to return if incx .le. 0.
1452 c modified 12/3/93, array(1) declarations changed to array(*)
1453 c
1454  double precision dx(*),dmax
1455  integer i,incx,ix,n
1456 c
1457  idamax = 0
1458  if( n.lt.1 .or. incx.le.0 ) return
1459  idamax = 1
1460  if(n.eq.1)return
1461  if(incx.eq.1)go to 20
1462 c
1463 c code for increment not equal to 1
1464 c
1465  ix = 1
1466  dmax = dabs(dx(1))
1467  ix = ix + incx
1468  do 10 i = 2,n
1469  if(dabs(dx(ix)).le.dmax) go to 5
1470  idamax = i
1471  dmax = dabs(dx(ix))
1472  5 ix = ix + incx
1473  10 continue
1474  return
1475 c
1476 c code for increment equal to 1
1477 c
1478  20 dmax = dabs(dx(1))
1479  do 30 i = 2,n
1480  if(dabs(dx(i)).le.dmax) go to 30
1481  idamax = i
1482  dmax = dabs(dx(i))
1483  30 continue
1484  return
idamax
integer function idamax(n, dx, incx)
Definition: dgesv.f:1448

◆ ieeeck()

integer function ieeeck ( integer  ISPEC,
real  ZERO,
real  ONE 
)
731 *
732 * -- LAPACK auxiliary routine (version 3.0) --
733 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
734 * Courant Institute, Argonne National Lab, and Rice University
735 * June 30, 1998
736 *
737 * .. Scalar Arguments ..
738  INTEGER ispec
739  REAL one, zero
740 * ..
741 *
742 * Purpose
743 * =======
744 *
745 * IEEECK is called from the ILAENV to verify that Infinity and
746 * possibly NaN arithmetic is safe (i.e. will not trap).
747 *
748 * Arguments
749 * =========
750 *
751 * ISPEC (input) INTEGER
752 * Specifies whether to test just for inifinity arithmetic
753 * or whether to test for infinity and NaN arithmetic.
754 * = 0: Verify infinity arithmetic only.
755 * = 1: Verify infinity and NaN arithmetic.
756 *
757 * ZERO (input) REAL
758 * Must contain the value 0.0
759 * This is passed to prevent the compiler from optimizing
760 * away this code.
761 *
762 * ONE (input) REAL
763 * Must contain the value 1.0
764 * This is passed to prevent the compiler from optimizing
765 * away this code.
766 *
767 * RETURN VALUE: INTEGER
768 * = 0: Arithmetic failed to produce the correct answers
769 * = 1: Arithmetic produced the correct answers
770 *
771 * .. Local Scalars ..
772  REAL nan1, nan2, nan3, nan4, nan5, nan6, neginf,
773  $ negzro, newzro, posinf
774 * ..
775 * .. Executable Statements ..
776  ieeeck = 1
777 *
778  posinf = one / zero
779  IF( posinf.LE.one ) THEN
780  ieeeck = 0
781  RETURN
782  END IF
783 *
784  neginf = -one / zero
785  IF( neginf.GE.zero ) THEN
786  ieeeck = 0
787  RETURN
788  END IF
789 *
790  negzro = one / ( neginf+one )
791  IF( negzro.NE.zero ) THEN
792  ieeeck = 0
793  RETURN
794  END IF
795 *
796  neginf = one / negzro
797  IF( neginf.GE.zero ) THEN
798  ieeeck = 0
799  RETURN
800  END IF
801 *
802  newzro = negzro + zero
803  IF( newzro.NE.zero ) THEN
804  ieeeck = 0
805  RETURN
806  END IF
807 *
808  posinf = one / newzro
809  IF( posinf.LE.one ) THEN
810  ieeeck = 0
811  RETURN
812  END IF
813 *
814  neginf = neginf*posinf
815  IF( neginf.GE.zero ) THEN
816  ieeeck = 0
817  RETURN
818  END IF
819 *
820  posinf = posinf*posinf
821  IF( posinf.LE.one ) THEN
822  ieeeck = 0
823  RETURN
824  END IF
825 *
826 *
827 *
828 *
829 * Return if we were only asked to check infinity arithmetic
830 *
831  IF( ispec.EQ.0 )
832  $ RETURN
833 *
834  nan1 = posinf + neginf
835 *
836  nan2 = posinf / neginf
837 *
838  nan3 = posinf / posinf
839 *
840  nan4 = posinf*zero
841 *
842  nan5 = neginf*negzro
843 *
844  nan6 = nan5*0.0
845 *
846  IF( nan1.EQ.nan1 ) THEN
847  ieeeck = 0
848  RETURN
849  END IF
850 *
851  IF( nan2.EQ.nan2 ) THEN
852  ieeeck = 0
853  RETURN
854  END IF
855 *
856  IF( nan3.EQ.nan3 ) THEN
857  ieeeck = 0
858  RETURN
859  END IF
860 *
861  IF( nan4.EQ.nan4 ) THEN
862  ieeeck = 0
863  RETURN
864  END IF
865 *
866  IF( nan5.EQ.nan5 ) THEN
867  ieeeck = 0
868  RETURN
869  END IF
870 *
871  IF( nan6.EQ.nan6 ) THEN
872  ieeeck = 0
873  RETURN
874  END IF
875 *
876  RETURN
ieeeck
integer function ieeeck(ISPEC, ZERO, ONE)
Definition: dgesv.f:731

◆ ilaenv()

integer function ilaenv ( integer  ISPEC,
character*( * )  NAME,
character*( * )  OPTS,
integer  N1,
integer  N2,
integer  N3,
integer  N4 
)
880 *
881 * -- LAPACK auxiliary routine (version 3.0) --
882 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
883 * Courant Institute, Argonne National Lab, and Rice University
884 * June 30, 1999
885 *
886 * .. Scalar Arguments ..
887  CHARACTER*( * ) name, opts
888  INTEGER ispec, n1, n2, n3, n4
889 * ..
890 *
891 * Purpose
892 * =======
893 *
894 * ILAENV is called from the LAPACK routines to choose problem-dependent
895 * parameters for the local environment. See ISPEC for a description of
896 * the parameters.
897 *
898 * This version provides a set of parameters which should give good,
899 * but not optimal, performance on many of the currently available
900 * computers. Users are encouraged to modify this subroutine to set
901 * the tuning parameters for their particular machine using the option
902 * and problem size information in the arguments.
903 *
904 * This routine will not function correctly if it is converted to all
905 * lower case. Converting it to all upper case is allowed.
906 *
907 * Arguments
908 * =========
909 *
910 * ISPEC (input) INTEGER
911 * Specifies the parameter to be returned as the value of
912 * ILAENV.
913 * = 1: the optimal blocksize; if this value is 1, an unblocked
914 * algorithm will give the best performance.
915 * = 2: the minimum block size for which the block routine
916 * should be used; if the usable block size is less than
917 * this value, an unblocked routine should be used.
918 * = 3: the crossover point (in a block routine, for N less
919 * than this value, an unblocked routine should be used)
920 * = 4: the number of shifts, used in the nonsymmetric
921 * eigenvalue routines
922 * = 5: the minimum column dimension for blocking to be used;
923 * rectangular blocks must have dimension at least k by m,
924 * where k is given by ILAENV(2,...) and m by ILAENV(5,...)
925 * = 6: the crossover point for the SVD (when reducing an m by n
926 * matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
927 * this value, a QR factorization is used first to reduce
928 * the matrix to a triangular form.)
929 * = 7: the number of processors
930 * = 8: the crossover point for the multishift QR and QZ methods
931 * for nonsymmetric eigenvalue problems.
932 * = 9: maximum size of the subproblems at the bottom of the
933 * computation tree in the divide-and-conquer algorithm
934 * (used by xGELSD and xGESDD)
935 * =10: ieee NaN arithmetic can be trusted not to trap
936 * =11: infinity arithmetic can be trusted not to trap
937 *
938 * NAME (input) CHARACTER*(*)
939 * The name of the calling subroutine, in either upper case or
940 * lower case.
941 *
942 * OPTS (input) CHARACTER*(*)
943 * The character options to the subroutine NAME, concatenated
944 * into a single character string. For example, UPLO = 'U',
945 * TRANS = 'T', and DIAG = 'N' for a triangular routine would
946 * be specified as OPTS = 'UTN'.
947 *
948 * N1 (input) INTEGER
949 * N2 (input) INTEGER
950 * N3 (input) INTEGER
951 * N4 (input) INTEGER
952 * Problem dimensions for the subroutine NAME; these may not all
953 * be required.
954 *
955 * (ILAENV) (output) INTEGER
956 * >= 0: the value of the parameter specified by ISPEC
957 * < 0: if ILAENV = -k, the k-th argument had an illegal value.
958 *
959 * Further Details
960 * ===============
961 *
962 * The following conventions have been used when calling ILAENV from the
963 * LAPACK routines:
964 * 1) OPTS is a concatenation of all of the character options to
965 * subroutine NAME, in the same order that they appear in the
966 * argument list for NAME, even if they are not used in determining
967 * the value of the parameter specified by ISPEC.
968 * 2) The problem dimensions N1, N2, N3, N4 are specified in the order
969 * that they appear in the argument list for NAME. N1 is used
970 * first, N2 second, and so on, and unused problem dimensions are
971 * passed a value of -1.
972 * 3) The parameter value returned by ILAENV is checked for validity in
973 * the calling subroutine. For example, ILAENV is used to retrieve
974 * the optimal blocksize for STRTRI as follows:
975 *
976 * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
977 * IF( NB.LE.1 ) NB = MAX( 1, N )
978 *
979 * =====================================================================
980 *
981 * .. Local Scalars ..
982  LOGICAL cname, sname
983  CHARACTER*1 c1
984  CHARACTER*2 c2, c4
985  CHARACTER*3 c3
986  CHARACTER*6 subnam
987  INTEGER i, ic, iz, nb, nbmin, nx
988 * ..
989 * .. Intrinsic Functions ..
990  INTRINSIC char, ichar, int, min, real
991 * ..
992 * .. External Functions ..
993  INTEGER ieeeck
994  EXTERNAL ieeeck
995 * ..
996 * .. Executable Statements ..
997 *
998 C GO TO ( 100, 100, 100, 400, 500, 600, 700, 800, 900, 1000,
999 C $ 1100 ) ISPEC
1000 C CHANGED COMPUTED GOTO: OBSOLETE
1001 C
1002  IF((ispec.EQ.1).OR.(ispec.EQ.2).OR.(ispec.EQ.3)) THEN
1003  GO TO 100
1004  ELSEIF(ispec.EQ.4) THEN
1005  GO TO 400
1006  ELSEIF(ispec.EQ.5) THEN
1007  GO TO 500
1008  ELSEIF(ispec.EQ.6) THEN
1009  GO TO 600
1010  ELSEIF(ispec.EQ.7) THEN
1011  GO TO 700
1012  ELSEIF(ispec.EQ.8) THEN
1013  GO TO 800
1014  ELSEIF(ispec.EQ.9) THEN
1015  GO TO 900
1016  ELSEIF(ispec.EQ.10) THEN
1017  GO TO 1000
1018  ELSEIF(ispec.EQ.11) THEN
1019  GO TO 1100
1020  ENDIF
1021 *
1022 * Invalid value for ISPEC
1023 *
1024  ilaenv = -1
1025  RETURN
1026 *
1027  100 CONTINUE
1028 *
1029 * Convert NAME to upper case if the first character is lower case.
1030 *
1031  ilaenv = 1
1032  subnam = name
1033  ic = ichar( subnam( 1:1 ) )
1034  iz = ichar( 'Z' )
1035  IF( iz.EQ.90 .OR. iz.EQ.122 ) THEN
1036 *
1037 * ASCII character set
1038 *
1039  IF( ic.GE.97 .AND. ic.LE.122 ) THEN
1040  subnam( 1:1 ) = char( ic-32 )
1041  DO 10 i = 2, 6
1042  ic = ichar( subnam( i:i ) )
1043  IF( ic.GE.97 .AND. ic.LE.122 )
1044  $ subnam( i:i ) = char( ic-32 )
1045  10 CONTINUE
1046  END IF
1047 *
1048  ELSE IF( iz.EQ.233 .OR. iz.EQ.169 ) THEN
1049 *
1050 * EBCDIC character set
1051 *
1052  IF( ( ic.GE.129 .AND. ic.LE.137 ) .OR.
1053  $ ( ic.GE.145 .AND. ic.LE.153 ) .OR.
1054  $ ( ic.GE.162 .AND. ic.LE.169 ) ) THEN
1055  subnam( 1:1 ) = char( ic+64 )
1056  DO 20 i = 2, 6
1057  ic = ichar( subnam( i:i ) )
1058  IF( ( ic.GE.129 .AND. ic.LE.137 ) .OR.
1059  $ ( ic.GE.145 .AND. ic.LE.153 ) .OR.
1060  $ ( ic.GE.162 .AND. ic.LE.169 ) )
1061  $ subnam( i:i ) = char( ic+64 )
1062  20 CONTINUE
1063  END IF
1064 *
1065  ELSE IF( iz.EQ.218 .OR. iz.EQ.250 ) THEN
1066 *
1067 * Prime machines: ASCII+128
1068 *
1069  IF( ic.GE.225 .AND. ic.LE.250 ) THEN
1070  subnam( 1:1 ) = char( ic-32 )
1071  DO 30 i = 2, 6
1072  ic = ichar( subnam( i:i ) )
1073  IF( ic.GE.225 .AND. ic.LE.250 )
1074  $ subnam( i:i ) = char( ic-32 )
1075  30 CONTINUE
1076  END IF
1077  END IF
1078 *
1079  c1 = subnam( 1:1 )
1080  sname = c1.EQ.'S' .OR. c1.EQ.'D'
1081  cname = c1.EQ.'C' .OR. c1.EQ.'Z'
1082  IF( .NOT.( cname .OR. sname ) )
1083  $ RETURN
1084  c2 = subnam( 2:3 )
1085  c3 = subnam( 4:6 )
1086  c4 = c3( 2:3 )
1087 *
1088  GO TO ( 110, 200, 300 ) ispec
1089 *
1090  110 CONTINUE
1091 *
1092 * ISPEC = 1: block size
1093 *
1094 * In these examples, separate code is provided for setting NB for
1095 * real and complex. We assume that NB will take the same value in
1096 * single or double precision.
1097 *
1098  nb = 1
1099 *
1100  IF( c2.EQ.'GE' ) THEN
1101  IF( c3.EQ.'TRF' ) THEN
1102  IF( sname ) THEN
1103  nb = 64
1104  ELSE
1105  nb = 64
1106  END IF
1107  ELSE IF( c3.EQ.'QRF' .OR. c3.EQ.'RQF' .OR. c3.EQ.'LQF' .OR.
1108  $ c3.EQ.'QLF' ) THEN
1109  IF( sname ) THEN
1110  nb = 32
1111  ELSE
1112  nb = 32
1113  END IF
1114  ELSE IF( c3.EQ.'HRD' ) THEN
1115  IF( sname ) THEN
1116  nb = 32
1117  ELSE
1118  nb = 32
1119  END IF
1120  ELSE IF( c3.EQ.'BRD' ) THEN
1121  IF( sname ) THEN
1122  nb = 32
1123  ELSE
1124  nb = 32
1125  END IF
1126  ELSE IF( c3.EQ.'TRI' ) THEN
1127  IF( sname ) THEN
1128  nb = 64
1129  ELSE
1130  nb = 64
1131  END IF
1132  END IF
1133  ELSE IF( c2.EQ.'PO' ) THEN
1134  IF( c3.EQ.'TRF' ) THEN
1135  IF( sname ) THEN
1136  nb = 64
1137  ELSE
1138  nb = 64
1139  END IF
1140  END IF
1141  ELSE IF( c2.EQ.'SY' ) THEN
1142  IF( c3.EQ.'TRF' ) THEN
1143  IF( sname ) THEN
1144  nb = 64
1145  ELSE
1146  nb = 64
1147  END IF
1148  ELSE IF( sname .AND. c3.EQ.'TRD' ) THEN
1149  nb = 32
1150  ELSE IF( sname .AND. c3.EQ.'GST' ) THEN
1151  nb = 64
1152  END IF
1153  ELSE IF( cname .AND. c2.EQ.'HE' ) THEN
1154  IF( c3.EQ.'TRF' ) THEN
1155  nb = 64
1156  ELSE IF( c3.EQ.'TRD' ) THEN
1157  nb = 32
1158  ELSE IF( c3.EQ.'GST' ) THEN
1159  nb = 64
1160  END IF
1161  ELSE IF( sname .AND. c2.EQ.'OR' ) THEN
1162  IF( c3( 1:1 ).EQ.'G' ) THEN
1163  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1164  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1165  $ c4.EQ.'BR' ) THEN
1166  nb = 32
1167  END IF
1168  ELSE IF( c3( 1:1 ).EQ.'M' ) THEN
1169  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1170  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1171  $ c4.EQ.'BR' ) THEN
1172  nb = 32
1173  END IF
1174  END IF
1175  ELSE IF( cname .AND. c2.EQ.'UN' ) THEN
1176  IF( c3( 1:1 ).EQ.'G' ) THEN
1177  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1178  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1179  $ c4.EQ.'BR' ) THEN
1180  nb = 32
1181  END IF
1182  ELSE IF( c3( 1:1 ).EQ.'M' ) THEN
1183  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1184  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1185  $ c4.EQ.'BR' ) THEN
1186  nb = 32
1187  END IF
1188  END IF
1189  ELSE IF( c2.EQ.'GB' ) THEN
1190  IF( c3.EQ.'TRF' ) THEN
1191  IF( sname ) THEN
1192  IF( n4.LE.64 ) THEN
1193  nb = 1
1194  ELSE
1195  nb = 32
1196  END IF
1197  ELSE
1198  IF( n4.LE.64 ) THEN
1199  nb = 1
1200  ELSE
1201  nb = 32
1202  END IF
1203  END IF
1204  END IF
1205  ELSE IF( c2.EQ.'PB' ) THEN
1206  IF( c3.EQ.'TRF' ) THEN
1207  IF( sname ) THEN
1208  IF( n2.LE.64 ) THEN
1209  nb = 1
1210  ELSE
1211  nb = 32
1212  END IF
1213  ELSE
1214  IF( n2.LE.64 ) THEN
1215  nb = 1
1216  ELSE
1217  nb = 32
1218  END IF
1219  END IF
1220  END IF
1221  ELSE IF( c2.EQ.'TR' ) THEN
1222  IF( c3.EQ.'TRI' ) THEN
1223  IF( sname ) THEN
1224  nb = 64
1225  ELSE
1226  nb = 64
1227  END IF
1228  END IF
1229  ELSE IF( c2.EQ.'LA' ) THEN
1230  IF( c3.EQ.'UUM' ) THEN
1231  IF( sname ) THEN
1232  nb = 64
1233  ELSE
1234  nb = 64
1235  END IF
1236  END IF
1237  ELSE IF( sname .AND. c2.EQ.'ST' ) THEN
1238  IF( c3.EQ.'EBZ' ) THEN
1239  nb = 1
1240  END IF
1241  END IF
1242  ilaenv = nb
1243  RETURN
1244 *
1245  200 CONTINUE
1246 *
1247 * ISPEC = 2: minimum block size
1248 *
1249  nbmin = 2
1250  IF( c2.EQ.'GE' ) THEN
1251  IF( c3.EQ.'QRF' .OR. c3.EQ.'RQF' .OR. c3.EQ.'LQF' .OR.
1252  $ c3.EQ.'QLF' ) THEN
1253  IF( sname ) THEN
1254  nbmin = 2
1255  ELSE
1256  nbmin = 2
1257  END IF
1258  ELSE IF( c3.EQ.'HRD' ) THEN
1259  IF( sname ) THEN
1260  nbmin = 2
1261  ELSE
1262  nbmin = 2
1263  END IF
1264  ELSE IF( c3.EQ.'BRD' ) THEN
1265  IF( sname ) THEN
1266  nbmin = 2
1267  ELSE
1268  nbmin = 2
1269  END IF
1270  ELSE IF( c3.EQ.'TRI' ) THEN
1271  IF( sname ) THEN
1272  nbmin = 2
1273  ELSE
1274  nbmin = 2
1275  END IF
1276  END IF
1277  ELSE IF( c2.EQ.'SY' ) THEN
1278  IF( c3.EQ.'TRF' ) THEN
1279  IF( sname ) THEN
1280  nbmin = 8
1281  ELSE
1282  nbmin = 8
1283  END IF
1284  ELSE IF( sname .AND. c3.EQ.'TRD' ) THEN
1285  nbmin = 2
1286  END IF
1287  ELSE IF( cname .AND. c2.EQ.'HE' ) THEN
1288  IF( c3.EQ.'TRD' ) THEN
1289  nbmin = 2
1290  END IF
1291  ELSE IF( sname .AND. c2.EQ.'OR' ) THEN
1292  IF( c3( 1:1 ).EQ.'G' ) THEN
1293  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1294  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1295  $ c4.EQ.'BR' ) THEN
1296  nbmin = 2
1297  END IF
1298  ELSE IF( c3( 1:1 ).EQ.'M' ) THEN
1299  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1300  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1301  $ c4.EQ.'BR' ) THEN
1302  nbmin = 2
1303  END IF
1304  END IF
1305  ELSE IF( cname .AND. c2.EQ.'UN' ) THEN
1306  IF( c3( 1:1 ).EQ.'G' ) THEN
1307  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1308  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1309  $ c4.EQ.'BR' ) THEN
1310  nbmin = 2
1311  END IF
1312  ELSE IF( c3( 1:1 ).EQ.'M' ) THEN
1313  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1314  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1315  $ c4.EQ.'BR' ) THEN
1316  nbmin = 2
1317  END IF
1318  END IF
1319  END IF
1320  ilaenv = nbmin
1321  RETURN
1322 *
1323  300 CONTINUE
1324 *
1325 * ISPEC = 3: crossover point
1326 *
1327  nx = 0
1328  IF( c2.EQ.'GE' ) THEN
1329  IF( c3.EQ.'QRF' .OR. c3.EQ.'RQF' .OR. c3.EQ.'LQF' .OR.
1330  $ c3.EQ.'QLF' ) THEN
1331  IF( sname ) THEN
1332  nx = 128
1333  ELSE
1334  nx = 128
1335  END IF
1336  ELSE IF( c3.EQ.'HRD' ) THEN
1337  IF( sname ) THEN
1338  nx = 128
1339  ELSE
1340  nx = 128
1341  END IF
1342  ELSE IF( c3.EQ.'BRD' ) THEN
1343  IF( sname ) THEN
1344  nx = 128
1345  ELSE
1346  nx = 128
1347  END IF
1348  END IF
1349  ELSE IF( c2.EQ.'SY' ) THEN
1350  IF( sname .AND. c3.EQ.'TRD' ) THEN
1351  nx = 32
1352  END IF
1353  ELSE IF( cname .AND. c2.EQ.'HE' ) THEN
1354  IF( c3.EQ.'TRD' ) THEN
1355  nx = 32
1356  END IF
1357  ELSE IF( sname .AND. c2.EQ.'OR' ) THEN
1358  IF( c3( 1:1 ).EQ.'G' ) THEN
1359  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1360  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1361  $ c4.EQ.'BR' ) THEN
1362  nx = 128
1363  END IF
1364  END IF
1365  ELSE IF( cname .AND. c2.EQ.'UN' ) THEN
1366  IF( c3( 1:1 ).EQ.'G' ) THEN
1367  IF( c4.EQ.'QR' .OR. c4.EQ.'RQ' .OR. c4.EQ.'LQ' .OR.
1368  $ c4.EQ.'QL' .OR. c4.EQ.'HR' .OR. c4.EQ.'TR' .OR.
1369  $ c4.EQ.'BR' ) THEN
1370  nx = 128
1371  END IF
1372  END IF
1373  END IF
1374  ilaenv = nx
1375  RETURN
1376 *
1377  400 CONTINUE
1378 *
1379 * ISPEC = 4: number of shifts (used by xHSEQR)
1380 *
1381  ilaenv = 6
1382  RETURN
1383 *
1384  500 CONTINUE
1385 *
1386 * ISPEC = 5: minimum column dimension (not used)
1387 *
1388  ilaenv = 2
1389  RETURN
1390 *
1391  600 CONTINUE
1392 *
1393 * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD)
1394 *
1395  ilaenv = int( REAL( MIN( N1, N2 ) )*1.6e0 )
1396  RETURN
1397 *
1398  700 CONTINUE
1399 *
1400 * ISPEC = 7: number of processors (not used)
1401 *
1402  ilaenv = 1
1403  RETURN
1404 *
1405  800 CONTINUE
1406 *
1407 * ISPEC = 8: crossover point for multishift (used by xHSEQR)
1408 *
1409  ilaenv = 50
1410  RETURN
1411 *
1412  900 CONTINUE
1413 *
1414 * ISPEC = 9: maximum size of the subproblems at the bottom of the
1415 * computation tree in the divide-and-conquer algorithm
1416 * (used by xGELSD and xGESDD)
1417 *
1418  ilaenv = 25
1419  RETURN
1420 *
1421  1000 CONTINUE
1422 *
1423 * ISPEC = 10: ieee NaN arithmetic can be trusted not to trap
1424 *
1425 c ILAENV = 0
1426  ilaenv = 1
1427  IF( ilaenv.EQ.1 ) THEN
1428  ilaenv = ieeeck( 0, 0.0, 1.0 )
1429  END IF
1430  RETURN
1431 *
1432  1100 CONTINUE
1433 *
1434 * ISPEC = 11: infinity arithmetic can be trusted not to trap
1435 *
1436 c ILAENV = 0
1437  ilaenv = 1
1438  IF( ilaenv.EQ.1 ) THEN
1439  ilaenv = ieeeck( 1, 0.0, 1.0 )
1440  END IF
1441  RETURN
1442 *
1443 * End of ILAENV
1444 *
min
#define min(a, b)
Definition: cascade.c:31
c1
static double * c1
Definition: mafillvcompmain.c:30
ieeeck
integer function ieeeck(ISPEC, ZERO, ONE)
Definition: dgesv.f:731
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: dgesv.f:880
Hosted by OpenAircraft.com, (Michigan UAV, LLC)