You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

cgerfsx.f 28 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734
  1. *> \brief \b CGERFSX
  2. *
  3. * =========== DOCUMENTATION ===========
  4. *
  5. * Online html documentation available at
  6. * http://www.netlib.org/lapack/explore-html/
  7. *
  8. *> \htmlonly
  9. *> Download CGERFSX + dependencies
  10. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgerfsx.f">
  11. *> [TGZ]</a>
  12. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgerfsx.f">
  13. *> [ZIP]</a>
  14. *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgerfsx.f">
  15. *> [TXT]</a>
  16. *> \endhtmlonly
  17. *
  18. * Definition:
  19. * ===========
  20. *
  21. * SUBROUTINE CGERFSX( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
  22. * R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
  23. * ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
  24. * WORK, RWORK, INFO )
  25. *
  26. * .. Scalar Arguments ..
  27. * CHARACTER TRANS, EQUED
  28. * INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
  29. * $ N_ERR_BNDS
  30. * REAL RCOND
  31. * ..
  32. * .. Array Arguments ..
  33. * INTEGER IPIV( * )
  34. * COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  35. * $ X( LDX , * ), WORK( * )
  36. * REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
  37. * $ ERR_BNDS_NORM( NRHS, * ),
  38. * $ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
  39. * ..
  40. *
  41. *
  42. *> \par Purpose:
  43. * =============
  44. *>
  45. *> \verbatim
  46. *>
  47. *> CGERFSX improves the computed solution to a system of linear
  48. *> equations and provides error bounds and backward error estimates
  49. *> for the solution. In addition to normwise error bound, the code
  50. *> provides maximum componentwise error bound if possible. See
  51. *> comments for ERR_BNDS_NORM and ERR_BNDS_COMP for details of the
  52. *> error bounds.
  53. *>
  54. *> The original system of linear equations may have been equilibrated
  55. *> before calling this routine, as described by arguments EQUED, R
  56. *> and C below. In this case, the solution and error bounds returned
  57. *> are for the original unequilibrated system.
  58. *> \endverbatim
  59. *
  60. * Arguments:
  61. * ==========
  62. *
  63. *> \verbatim
  64. *> Some optional parameters are bundled in the PARAMS array. These
  65. *> settings determine how refinement is performed, but often the
  66. *> defaults are acceptable. If the defaults are acceptable, users
  67. *> can pass NPARAMS = 0 which prevents the source code from accessing
  68. *> the PARAMS argument.
  69. *> \endverbatim
  70. *>
  71. *> \param[in] TRANS
  72. *> \verbatim
  73. *> TRANS is CHARACTER*1
  74. *> Specifies the form of the system of equations:
  75. *> = 'N': A * X = B (No transpose)
  76. *> = 'T': A**T * X = B (Transpose)
  77. *> = 'C': A**H * X = B (Conjugate transpose)
  78. *> \endverbatim
  79. *>
  80. *> \param[in] EQUED
  81. *> \verbatim
  82. *> EQUED is CHARACTER*1
  83. *> Specifies the form of equilibration that was done to A
  84. *> before calling this routine. This is needed to compute
  85. *> the solution and error bounds correctly.
  86. *> = 'N': No equilibration
  87. *> = 'R': Row equilibration, i.e., A has been premultiplied by
  88. *> diag(R).
  89. *> = 'C': Column equilibration, i.e., A has been postmultiplied
  90. *> by diag(C).
  91. *> = 'B': Both row and column equilibration, i.e., A has been
  92. *> replaced by diag(R) * A * diag(C).
  93. *> The right hand side B has been changed accordingly.
  94. *> \endverbatim
  95. *>
  96. *> \param[in] N
  97. *> \verbatim
  98. *> N is INTEGER
  99. *> The order of the matrix A. N >= 0.
  100. *> \endverbatim
  101. *>
  102. *> \param[in] NRHS
  103. *> \verbatim
  104. *> NRHS is INTEGER
  105. *> The number of right hand sides, i.e., the number of columns
  106. *> of the matrices B and X. NRHS >= 0.
  107. *> \endverbatim
  108. *>
  109. *> \param[in] A
  110. *> \verbatim
  111. *> A is COMPLEX array, dimension (LDA,N)
  112. *> The original N-by-N matrix A.
  113. *> \endverbatim
  114. *>
  115. *> \param[in] LDA
  116. *> \verbatim
  117. *> LDA is INTEGER
  118. *> The leading dimension of the array A. LDA >= max(1,N).
  119. *> \endverbatim
  120. *>
  121. *> \param[in] AF
  122. *> \verbatim
  123. *> AF is COMPLEX array, dimension (LDAF,N)
  124. *> The factors L and U from the factorization A = P*L*U
  125. *> as computed by CGETRF.
  126. *> \endverbatim
  127. *>
  128. *> \param[in] LDAF
  129. *> \verbatim
  130. *> LDAF is INTEGER
  131. *> The leading dimension of the array AF. LDAF >= max(1,N).
  132. *> \endverbatim
  133. *>
  134. *> \param[in] IPIV
  135. *> \verbatim
  136. *> IPIV is INTEGER array, dimension (N)
  137. *> The pivot indices from CGETRF; for 1<=i<=N, row i of the
  138. *> matrix was interchanged with row IPIV(i).
  139. *> \endverbatim
  140. *>
  141. *> \param[in] R
  142. *> \verbatim
  143. *> R is REAL array, dimension (N)
  144. *> The row scale factors for A. If EQUED = 'R' or 'B', A is
  145. *> multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
  146. *> is not accessed.
  147. *> If R is accessed, each element of R should be a power of the radix
  148. *> to ensure a reliable solution and error estimates. Scaling by
  149. *> powers of the radix does not cause rounding errors unless the
  150. *> result underflows or overflows. Rounding errors during scaling
  151. *> lead to refining with a matrix that is not equivalent to the
  152. *> input matrix, producing error estimates that may not be
  153. *> reliable.
  154. *> \endverbatim
  155. *>
  156. *> \param[in] C
  157. *> \verbatim
  158. *> C is REAL array, dimension (N)
  159. *> The column scale factors for A. If EQUED = 'C' or 'B', A is
  160. *> multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
  161. *> is not accessed.
  162. *> If C is accessed, each element of C should be a power of the radix
  163. *> to ensure a reliable solution and error estimates. Scaling by
  164. *> powers of the radix does not cause rounding errors unless the
  165. *> result underflows or overflows. Rounding errors during scaling
  166. *> lead to refining with a matrix that is not equivalent to the
  167. *> input matrix, producing error estimates that may not be
  168. *> reliable.
  169. *> \endverbatim
  170. *>
  171. *> \param[in] B
  172. *> \verbatim
  173. *> B is COMPLEX array, dimension (LDB,NRHS)
  174. *> The right hand side matrix B.
  175. *> \endverbatim
  176. *>
  177. *> \param[in] LDB
  178. *> \verbatim
  179. *> LDB is INTEGER
  180. *> The leading dimension of the array B. LDB >= max(1,N).
  181. *> \endverbatim
  182. *>
  183. *> \param[in,out] X
  184. *> \verbatim
  185. *> X is COMPLEX array, dimension (LDX,NRHS)
  186. *> On entry, the solution matrix X, as computed by CGETRS.
  187. *> On exit, the improved solution matrix X.
  188. *> \endverbatim
  189. *>
  190. *> \param[in] LDX
  191. *> \verbatim
  192. *> LDX is INTEGER
  193. *> The leading dimension of the array X. LDX >= max(1,N).
  194. *> \endverbatim
  195. *>
  196. *> \param[out] RCOND
  197. *> \verbatim
  198. *> RCOND is REAL
  199. *> Reciprocal scaled condition number. This is an estimate of the
  200. *> reciprocal Skeel condition number of the matrix A after
  201. *> equilibration (if done). If this is less than the machine
  202. *> precision (in particular, if it is zero), the matrix is singular
  203. *> to working precision. Note that the error may still be small even
  204. *> if this number is very small and the matrix appears ill-
  205. *> conditioned.
  206. *> \endverbatim
  207. *>
  208. *> \param[out] BERR
  209. *> \verbatim
  210. *> BERR is REAL array, dimension (NRHS)
  211. *> Componentwise relative backward error. This is the
  212. *> componentwise relative backward error of each solution vector X(j)
  213. *> (i.e., the smallest relative change in any element of A or B that
  214. *> makes X(j) an exact solution).
  215. *> \endverbatim
  216. *>
  217. *> \param[in] N_ERR_BNDS
  218. *> \verbatim
  219. *> N_ERR_BNDS is INTEGER
  220. *> Number of error bounds to return for each right hand side
  221. *> and each type (normwise or componentwise). See ERR_BNDS_NORM and
  222. *> ERR_BNDS_COMP below.
  223. *> \endverbatim
  224. *>
  225. *> \param[out] ERR_BNDS_NORM
  226. *> \verbatim
  227. *> ERR_BNDS_NORM is REAL array, dimension (NRHS, N_ERR_BNDS)
  228. *> For each right-hand side, this array contains information about
  229. *> various error bounds and condition numbers corresponding to the
  230. *> normwise relative error, which is defined as follows:
  231. *>
  232. *> Normwise relative error in the ith solution vector:
  233. *> max_j (abs(XTRUE(j,i) - X(j,i)))
  234. *> ------------------------------
  235. *> max_j abs(X(j,i))
  236. *>
  237. *> The array is indexed by the type of error information as described
  238. *> below. There currently are up to three pieces of information
  239. *> returned.
  240. *>
  241. *> The first index in ERR_BNDS_NORM(i,:) corresponds to the ith
  242. *> right-hand side.
  243. *>
  244. *> The second index in ERR_BNDS_NORM(:,err) contains the following
  245. *> three fields:
  246. *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
  247. *> reciprocal condition number is less than the threshold
  248. *> sqrt(n) * slamch('Epsilon').
  249. *>
  250. *> err = 2 "Guaranteed" error bound: The estimated forward error,
  251. *> almost certainly within a factor of 10 of the true error
  252. *> so long as the next entry is greater than the threshold
  253. *> sqrt(n) * slamch('Epsilon'). This error bound should only
  254. *> be trusted if the previous boolean is true.
  255. *>
  256. *> err = 3 Reciprocal condition number: Estimated normwise
  257. *> reciprocal condition number. Compared with the threshold
  258. *> sqrt(n) * slamch('Epsilon') to determine if the error
  259. *> estimate is "guaranteed". These reciprocal condition
  260. *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
  261. *> appropriately scaled matrix Z.
  262. *> Let Z = S*A, where S scales each row by a power of the
  263. *> radix so all absolute row sums of Z are approximately 1.
  264. *>
  265. *> See Lapack Working Note 165 for further details and extra
  266. *> cautions.
  267. *> \endverbatim
  268. *>
  269. *> \param[out] ERR_BNDS_COMP
  270. *> \verbatim
  271. *> ERR_BNDS_COMP is REAL array, dimension (NRHS, N_ERR_BNDS)
  272. *> For each right-hand side, this array contains information about
  273. *> various error bounds and condition numbers corresponding to the
  274. *> componentwise relative error, which is defined as follows:
  275. *>
  276. *> Componentwise relative error in the ith solution vector:
  277. *> abs(XTRUE(j,i) - X(j,i))
  278. *> max_j ----------------------
  279. *> abs(X(j,i))
  280. *>
  281. *> The array is indexed by the right-hand side i (on which the
  282. *> componentwise relative error depends), and the type of error
  283. *> information as described below. There currently are up to three
  284. *> pieces of information returned for each right-hand side. If
  285. *> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
  286. *> ERR_BNDS_COMP is not accessed. If N_ERR_BNDS < 3, then at most
  287. *> the first (:,N_ERR_BNDS) entries are returned.
  288. *>
  289. *> The first index in ERR_BNDS_COMP(i,:) corresponds to the ith
  290. *> right-hand side.
  291. *>
  292. *> The second index in ERR_BNDS_COMP(:,err) contains the following
  293. *> three fields:
  294. *> err = 1 "Trust/don't trust" boolean. Trust the answer if the
  295. *> reciprocal condition number is less than the threshold
  296. *> sqrt(n) * slamch('Epsilon').
  297. *>
  298. *> err = 2 "Guaranteed" error bound: The estimated forward error,
  299. *> almost certainly within a factor of 10 of the true error
  300. *> so long as the next entry is greater than the threshold
  301. *> sqrt(n) * slamch('Epsilon'). This error bound should only
  302. *> be trusted if the previous boolean is true.
  303. *>
  304. *> err = 3 Reciprocal condition number: Estimated componentwise
  305. *> reciprocal condition number. Compared with the threshold
  306. *> sqrt(n) * slamch('Epsilon') to determine if the error
  307. *> estimate is "guaranteed". These reciprocal condition
  308. *> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
  309. *> appropriately scaled matrix Z.
  310. *> Let Z = S*(A*diag(x)), where x is the solution for the
  311. *> current right-hand side and S scales each row of
  312. *> A*diag(x) by a power of the radix so all absolute row
  313. *> sums of Z are approximately 1.
  314. *>
  315. *> See Lapack Working Note 165 for further details and extra
  316. *> cautions.
  317. *> \endverbatim
  318. *>
  319. *> \param[in] NPARAMS
  320. *> \verbatim
  321. *> NPARAMS is INTEGER
  322. *> Specifies the number of parameters set in PARAMS. If <= 0, the
  323. *> PARAMS array is never referenced and default values are used.
  324. *> \endverbatim
  325. *>
  326. *> \param[in,out] PARAMS
  327. *> \verbatim
  328. *> PARAMS is REAL array, dimension NPARAMS
  329. *> Specifies algorithm parameters. If an entry is < 0.0, then
  330. *> that entry will be filled with default value used for that
  331. *> parameter. Only positions up to NPARAMS are accessed; defaults
  332. *> are used for higher-numbered parameters.
  333. *>
  334. *> PARAMS(LA_LINRX_ITREF_I = 1) : Whether to perform iterative
  335. *> refinement or not.
  336. *> Default: 1.0
  337. *> = 0.0: No refinement is performed, and no error bounds are
  338. *> computed.
  339. *> = 1.0: Use the double-precision refinement algorithm,
  340. *> possibly with doubled-single computations if the
  341. *> compilation environment does not support DOUBLE
  342. *> PRECISION.
  343. *> (other values are reserved for future use)
  344. *>
  345. *> PARAMS(LA_LINRX_ITHRESH_I = 2) : Maximum number of residual
  346. *> computations allowed for refinement.
  347. *> Default: 10
  348. *> Aggressive: Set to 100 to permit convergence using approximate
  349. *> factorizations or factorizations other than LU. If
  350. *> the factorization uses a technique other than
  351. *> Gaussian elimination, the guarantees in
  352. *> err_bnds_norm and err_bnds_comp may no longer be
  353. *> trustworthy.
  354. *>
  355. *> PARAMS(LA_LINRX_CWISE_I = 3) : Flag determining if the code
  356. *> will attempt to find a solution with small componentwise
  357. *> relative error in the double-precision algorithm. Positive
  358. *> is true, 0.0 is false.
  359. *> Default: 1.0 (attempt componentwise convergence)
  360. *> \endverbatim
  361. *>
  362. *> \param[out] WORK
  363. *> \verbatim
  364. *> WORK is COMPLEX array, dimension (2*N)
  365. *> \endverbatim
  366. *>
  367. *> \param[out] RWORK
  368. *> \verbatim
  369. *> RWORK is REAL array, dimension (2*N)
  370. *> \endverbatim
  371. *>
  372. *> \param[out] INFO
  373. *> \verbatim
  374. *> INFO is INTEGER
  375. *> = 0: Successful exit. The solution to every right-hand side is
  376. *> guaranteed.
  377. *> < 0: If INFO = -i, the i-th argument had an illegal value
  378. *> > 0 and <= N: U(INFO,INFO) is exactly zero. The factorization
  379. *> has been completed, but the factor U is exactly singular, so
  380. *> the solution and error bounds could not be computed. RCOND = 0
  381. *> is returned.
  382. *> = N+J: The solution corresponding to the Jth right-hand side is
  383. *> not guaranteed. The solutions corresponding to other right-
  384. *> hand sides K with K > J may not be guaranteed as well, but
  385. *> only the first such right-hand side is reported. If a small
  386. *> componentwise error is not requested (PARAMS(3) = 0.0) then
  387. *> the Jth right-hand side is the first with a normwise error
  388. *> bound that is not guaranteed (the smallest J such
  389. *> that ERR_BNDS_NORM(J,1) = 0.0). By default (PARAMS(3) = 1.0)
  390. *> the Jth right-hand side is the first with either a normwise or
  391. *> componentwise error bound that is not guaranteed (the smallest
  392. *> J such that either ERR_BNDS_NORM(J,1) = 0.0 or
  393. *> ERR_BNDS_COMP(J,1) = 0.0). See the definition of
  394. *> ERR_BNDS_NORM(:,1) and ERR_BNDS_COMP(:,1). To get information
  395. *> about all of the right-hand sides check ERR_BNDS_NORM or
  396. *> ERR_BNDS_COMP.
  397. *> \endverbatim
  398. *
  399. * Authors:
  400. * ========
  401. *
  402. *> \author Univ. of Tennessee
  403. *> \author Univ. of California Berkeley
  404. *> \author Univ. of Colorado Denver
  405. *> \author NAG Ltd.
  406. *
  407. *> \date December 2016
  408. *
  409. *> \ingroup complexGEcomputational
  410. *
  411. * =====================================================================
  412. SUBROUTINE CGERFSX( TRANS, EQUED, N, NRHS, A, LDA, AF, LDAF, IPIV,
  413. $ R, C, B, LDB, X, LDX, RCOND, BERR, N_ERR_BNDS,
  414. $ ERR_BNDS_NORM, ERR_BNDS_COMP, NPARAMS, PARAMS,
  415. $ WORK, RWORK, INFO )
  416. *
  417. * -- LAPACK computational routine (version 3.7.0) --
  418. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  419. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  420. * December 2016
  421. *
  422. * .. Scalar Arguments ..
  423. CHARACTER TRANS, EQUED
  424. INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS, NPARAMS,
  425. $ N_ERR_BNDS
  426. REAL RCOND
  427. * ..
  428. * .. Array Arguments ..
  429. INTEGER IPIV( * )
  430. COMPLEX A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
  431. $ X( LDX , * ), WORK( * )
  432. REAL R( * ), C( * ), PARAMS( * ), BERR( * ),
  433. $ ERR_BNDS_NORM( NRHS, * ),
  434. $ ERR_BNDS_COMP( NRHS, * ), RWORK( * )
  435. * ..
  436. *
  437. * ==================================================================
  438. *
  439. * .. Parameters ..
  440. REAL ZERO, ONE
  441. PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  442. REAL ITREF_DEFAULT, ITHRESH_DEFAULT,
  443. $ COMPONENTWISE_DEFAULT
  444. REAL RTHRESH_DEFAULT, DZTHRESH_DEFAULT
  445. PARAMETER ( ITREF_DEFAULT = 1.0 )
  446. PARAMETER ( ITHRESH_DEFAULT = 10.0 )
  447. PARAMETER ( COMPONENTWISE_DEFAULT = 1.0 )
  448. PARAMETER ( RTHRESH_DEFAULT = 0.5 )
  449. PARAMETER ( DZTHRESH_DEFAULT = 0.25 )
  450. INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
  451. $ LA_LINRX_CWISE_I
  452. PARAMETER ( LA_LINRX_ITREF_I = 1,
  453. $ LA_LINRX_ITHRESH_I = 2 )
  454. PARAMETER ( LA_LINRX_CWISE_I = 3 )
  455. INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
  456. $ LA_LINRX_RCOND_I
  457. PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
  458. PARAMETER ( LA_LINRX_RCOND_I = 3 )
  459. * ..
  460. * .. Local Scalars ..
  461. CHARACTER(1) NORM
  462. LOGICAL ROWEQU, COLEQU, NOTRAN
  463. INTEGER J, TRANS_TYPE, PREC_TYPE, REF_TYPE
  464. INTEGER N_NORMS
  465. REAL ANORM, RCOND_TMP
  466. REAL ILLRCOND_THRESH, ERR_LBND, CWISE_WRONG
  467. LOGICAL IGNORE_CWISE
  468. INTEGER ITHRESH
  469. REAL RTHRESH, UNSTABLE_THRESH
  470. * ..
  471. * .. External Subroutines ..
  472. EXTERNAL XERBLA, CGECON, CLA_GERFSX_EXTENDED
  473. * ..
  474. * .. Intrinsic Functions ..
  475. INTRINSIC MAX, SQRT, TRANSFER
  476. * ..
  477. * .. External Functions ..
  478. EXTERNAL LSAME, ILATRANS, ILAPREC
  479. EXTERNAL SLAMCH, CLANGE, CLA_GERCOND_X, CLA_GERCOND_C
  480. REAL SLAMCH, CLANGE, CLA_GERCOND_X, CLA_GERCOND_C
  481. LOGICAL LSAME
  482. INTEGER ILATRANS, ILAPREC
  483. * ..
  484. * .. Executable Statements ..
  485. *
  486. * Check the input parameters.
  487. *
  488. INFO = 0
  489. TRANS_TYPE = ILATRANS( TRANS )
  490. REF_TYPE = INT( ITREF_DEFAULT )
  491. IF ( NPARAMS .GE. LA_LINRX_ITREF_I ) THEN
  492. IF ( PARAMS( LA_LINRX_ITREF_I ) .LT. 0.0 ) THEN
  493. PARAMS( LA_LINRX_ITREF_I ) = ITREF_DEFAULT
  494. ELSE
  495. REF_TYPE = PARAMS( LA_LINRX_ITREF_I )
  496. END IF
  497. END IF
  498. *
  499. * Set default parameters.
  500. *
  501. ILLRCOND_THRESH = REAL( N ) * SLAMCH( 'Epsilon' )
  502. ITHRESH = INT( ITHRESH_DEFAULT )
  503. RTHRESH = RTHRESH_DEFAULT
  504. UNSTABLE_THRESH = DZTHRESH_DEFAULT
  505. IGNORE_CWISE = COMPONENTWISE_DEFAULT .EQ. 0.0
  506. *
  507. IF ( NPARAMS.GE.LA_LINRX_ITHRESH_I ) THEN
  508. IF ( PARAMS( LA_LINRX_ITHRESH_I ).LT.0.0 ) THEN
  509. PARAMS(LA_LINRX_ITHRESH_I) = ITHRESH
  510. ELSE
  511. ITHRESH = INT( PARAMS( LA_LINRX_ITHRESH_I ) )
  512. END IF
  513. END IF
  514. IF ( NPARAMS.GE.LA_LINRX_CWISE_I ) THEN
  515. IF ( PARAMS( LA_LINRX_CWISE_I ).LT.0.0 ) THEN
  516. IF ( IGNORE_CWISE ) THEN
  517. PARAMS( LA_LINRX_CWISE_I ) = 0.0
  518. ELSE
  519. PARAMS( LA_LINRX_CWISE_I ) = 1.0
  520. END IF
  521. ELSE
  522. IGNORE_CWISE = PARAMS( LA_LINRX_CWISE_I ) .EQ. 0.0
  523. END IF
  524. END IF
  525. IF ( REF_TYPE .EQ. 0 .OR. N_ERR_BNDS .EQ. 0 ) THEN
  526. N_NORMS = 0
  527. ELSE IF ( IGNORE_CWISE ) THEN
  528. N_NORMS = 1
  529. ELSE
  530. N_NORMS = 2
  531. END IF
  532. *
  533. NOTRAN = LSAME( TRANS, 'N' )
  534. ROWEQU = LSAME( EQUED, 'R' ) .OR. LSAME( EQUED, 'B' )
  535. COLEQU = LSAME( EQUED, 'C' ) .OR. LSAME( EQUED, 'B' )
  536. *
  537. * Test input parameters.
  538. *
  539. IF( TRANS_TYPE.EQ.-1 ) THEN
  540. INFO = -1
  541. ELSE IF( .NOT.ROWEQU .AND. .NOT.COLEQU .AND.
  542. $ .NOT.LSAME( EQUED, 'N' ) ) THEN
  543. INFO = -2
  544. ELSE IF( N.LT.0 ) THEN
  545. INFO = -3
  546. ELSE IF( NRHS.LT.0 ) THEN
  547. INFO = -4
  548. ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  549. INFO = -6
  550. ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
  551. INFO = -8
  552. ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
  553. INFO = -13
  554. ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
  555. INFO = -15
  556. END IF
  557. IF( INFO.NE.0 ) THEN
  558. CALL XERBLA( 'CGERFSX', -INFO )
  559. RETURN
  560. END IF
  561. *
  562. * Quick return if possible.
  563. *
  564. IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
  565. RCOND = 1.0
  566. DO J = 1, NRHS
  567. BERR( J ) = 0.0
  568. IF ( N_ERR_BNDS .GE. 1 ) THEN
  569. ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0
  570. ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0
  571. END IF
  572. IF ( N_ERR_BNDS .GE. 2 ) THEN
  573. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 0.0
  574. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 0.0
  575. END IF
  576. IF ( N_ERR_BNDS .GE. 3 ) THEN
  577. ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 1.0
  578. ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 1.0
  579. END IF
  580. END DO
  581. RETURN
  582. END IF
  583. *
  584. * Default to failure.
  585. *
  586. RCOND = 0.0
  587. DO J = 1, NRHS
  588. BERR( J ) = 1.0
  589. IF ( N_ERR_BNDS .GE. 1 ) THEN
  590. ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0
  591. ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0
  592. END IF
  593. IF ( N_ERR_BNDS .GE. 2 ) THEN
  594. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0
  595. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0
  596. END IF
  597. IF ( N_ERR_BNDS .GE. 3 ) THEN
  598. ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = 0.0
  599. ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = 0.0
  600. END IF
  601. END DO
  602. *
  603. * Compute the norm of A and the reciprocal of the condition
  604. * number of A.
  605. *
  606. IF( NOTRAN ) THEN
  607. NORM = 'I'
  608. ELSE
  609. NORM = '1'
  610. END IF
  611. ANORM = CLANGE( NORM, N, N, A, LDA, RWORK )
  612. CALL CGECON( NORM, N, AF, LDAF, ANORM, RCOND, WORK, RWORK, INFO )
  613. *
  614. * Perform refinement on each right-hand side
  615. *
  616. IF ( REF_TYPE .NE. 0 ) THEN
  617. PREC_TYPE = ILAPREC( 'D' )
  618. IF ( NOTRAN ) THEN
  619. CALL CLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N,
  620. $ NRHS, A, LDA, AF, LDAF, IPIV, COLEQU, C, B,
  621. $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
  622. $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
  623. $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
  624. $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
  625. $ INFO )
  626. ELSE
  627. CALL CLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N,
  628. $ NRHS, A, LDA, AF, LDAF, IPIV, ROWEQU, R, B,
  629. $ LDB, X, LDX, BERR, N_NORMS, ERR_BNDS_NORM,
  630. $ ERR_BNDS_COMP, WORK, RWORK, WORK(N+1),
  631. $ TRANSFER (RWORK(1:2*N), (/ (ZERO, ZERO) /), N),
  632. $ RCOND, ITHRESH, RTHRESH, UNSTABLE_THRESH, IGNORE_CWISE,
  633. $ INFO )
  634. END IF
  635. END IF
  636. ERR_LBND = MAX( 10.0, SQRT( REAL( N ) ) ) * SLAMCH( 'Epsilon' )
  637. IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 1 ) THEN
  638. *
  639. * Compute scaled normwise condition number cond(A*C).
  640. *
  641. IF ( COLEQU .AND. NOTRAN ) THEN
  642. RCOND_TMP = CLA_GERCOND_C( TRANS, N, A, LDA, AF, LDAF, IPIV,
  643. $ C, .TRUE., INFO, WORK, RWORK )
  644. ELSE IF ( ROWEQU .AND. .NOT. NOTRAN ) THEN
  645. RCOND_TMP = CLA_GERCOND_C( TRANS, N, A, LDA, AF, LDAF, IPIV,
  646. $ R, .TRUE., INFO, WORK, RWORK )
  647. ELSE
  648. RCOND_TMP = CLA_GERCOND_C( TRANS, N, A, LDA, AF, LDAF, IPIV,
  649. $ C, .FALSE., INFO, WORK, RWORK )
  650. END IF
  651. DO J = 1, NRHS
  652. *
  653. * Cap the error at 1.0.
  654. *
  655. IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
  656. $ .AND. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .GT. 1.0 )
  657. $ ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0
  658. *
  659. * Threshold the error (see LAWN).
  660. *
  661. IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
  662. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = 1.0
  663. ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 0.0
  664. IF ( INFO .LE. N ) INFO = N + J
  665. ELSE IF (ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) .LT. ERR_LBND)
  666. $ THEN
  667. ERR_BNDS_NORM( J, LA_LINRX_ERR_I ) = ERR_LBND
  668. ERR_BNDS_NORM( J, LA_LINRX_TRUST_I ) = 1.0
  669. END IF
  670. *
  671. * Save the condition number.
  672. *
  673. IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
  674. ERR_BNDS_NORM( J, LA_LINRX_RCOND_I ) = RCOND_TMP
  675. END IF
  676. END DO
  677. END IF
  678. IF ( N_ERR_BNDS .GE. 1 .AND. N_NORMS .GE. 2 ) THEN
  679. *
  680. * Compute componentwise condition number cond(A*diag(Y(:,J))) for
  681. * each right-hand side using the current solution as an estimate of
  682. * the true solution. If the componentwise error estimate is too
  683. * large, then the solution is a lousy estimate of truth and the
  684. * estimated RCOND may be too optimistic. To avoid misleading users,
  685. * the inverse condition number is set to 0.0 when the estimated
  686. * cwise error is at least CWISE_WRONG.
  687. *
  688. CWISE_WRONG = SQRT( SLAMCH( 'Epsilon' ) )
  689. DO J = 1, NRHS
  690. IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .LT. CWISE_WRONG )
  691. $ THEN
  692. RCOND_TMP = CLA_GERCOND_X( TRANS, N, A, LDA, AF, LDAF,
  693. $ IPIV, X(1,J), INFO, WORK, RWORK )
  694. ELSE
  695. RCOND_TMP = 0.0
  696. END IF
  697. *
  698. * Cap the error at 1.0.
  699. *
  700. IF ( N_ERR_BNDS .GE. LA_LINRX_ERR_I
  701. $ .AND. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) .GT. 1.0 )
  702. $ ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0
  703. *
  704. * Threshold the error (see LAWN).
  705. *
  706. IF ( RCOND_TMP .LT. ILLRCOND_THRESH ) THEN
  707. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = 1.0
  708. ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 0.0
  709. IF ( PARAMS( LA_LINRX_CWISE_I ) .EQ. 1.0
  710. $ .AND. INFO.LT.N + J ) INFO = N + J
  711. ELSE IF ( ERR_BNDS_COMP( J, LA_LINRX_ERR_I )
  712. $ .LT. ERR_LBND ) THEN
  713. ERR_BNDS_COMP( J, LA_LINRX_ERR_I ) = ERR_LBND
  714. ERR_BNDS_COMP( J, LA_LINRX_TRUST_I ) = 1.0
  715. END IF
  716. *
  717. * Save the condition number.
  718. *
  719. IF ( N_ERR_BNDS .GE. LA_LINRX_RCOND_I ) THEN
  720. ERR_BNDS_COMP( J, LA_LINRX_RCOND_I ) = RCOND_TMP
  721. END IF
  722. END DO
  723. END IF
  724. *
  725. RETURN
  726. *
  727. * End of CGERFSX
  728. *
  729. END