1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
|
*> \brief \b DLA_GERFSX_EXTENDED improves the computed solution to a system of linear equations for general matrices by performing extra-precise iterative refinement and provides error bounds and backward error estimates for the solution.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLA_GERFSX_EXTENDED + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dla_gerfsx_extended.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dla_gerfsx_extended.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dla_gerfsx_extended.f">
*> [TXT]</a>
*> \endhtmlonly
*
* Definition:
* ===========
*
* SUBROUTINE DLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, NRHS, A,
* LDA, AF, LDAF, IPIV, COLEQU, C, B,
* LDB, Y, LDY, BERR_OUT, N_NORMS,
* ERRS_N, ERRS_C, RES, AYB, DY,
* Y_TAIL, RCOND, ITHRESH, RTHRESH,
* DZ_UB, IGNORE_CWISE, INFO )
*
* .. Scalar Arguments ..
* INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
* $ TRANS_TYPE, N_NORMS, ITHRESH
* LOGICAL COLEQU, IGNORE_CWISE
* DOUBLE PRECISION RTHRESH, DZ_UB
* ..
* .. Array Arguments ..
* INTEGER IPIV( * )
* DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
* $ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
* DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ),
* $ ERRS_N( NRHS, * ), ERRS_C( NRHS, * )
* ..
*
*
*> \par Purpose:
* =============
*>
*> \verbatim
*>
*>
*> DLA_GERFSX_EXTENDED improves the computed solution to a system of
*> linear equations by performing extra-precise iterative refinement
*> and provides error bounds and backward error estimates for the solution.
*> This subroutine is called by DGERFSX to perform iterative refinement.
*> In addition to normwise error bound, the code provides maximum
*> componentwise error bound if possible. See comments for ERRS_N
*> and ERRS_C for details of the error bounds. Note that this
*> subroutine is only responsible for setting the second fields of
*> ERRS_N and ERRS_C.
*> \endverbatim
*
* Arguments:
* ==========
*
*> \param[in] PREC_TYPE
*> \verbatim
*> PREC_TYPE is INTEGER
*> Specifies the intermediate precision to be used in refinement.
*> The value is defined by ILAPREC(P) where P is a CHARACTER and P
*> = 'S': Single
*> = 'D': Double
*> = 'I': Indigenous
*> = 'X' or 'E': Extra
*> \endverbatim
*>
*> \param[in] TRANS_TYPE
*> \verbatim
*> TRANS_TYPE is INTEGER
*> Specifies the transposition operation on A.
*> The value is defined by ILATRANS(T) where T is a CHARACTER and T
*> = 'N': No transpose
*> = 'T': Transpose
*> = 'C': Conjugate transpose
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*> N is INTEGER
*> The number of linear equations, i.e., the order of the
*> matrix A. N >= 0.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*> NRHS is INTEGER
*> The number of right-hand-sides, i.e., the number of columns of the
*> matrix B.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*> A is DOUBLE PRECISION array, dimension (LDA,N)
*> On entry, the N-by-N matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*> LDA is INTEGER
*> The leading dimension of the array A. LDA >= max(1,N).
*> \endverbatim
*>
*> \param[in] AF
*> \verbatim
*> AF is DOUBLE PRECISION array, dimension (LDAF,N)
*> The factors L and U from the factorization
*> A = P*L*U as computed by DGETRF.
*> \endverbatim
*>
*> \param[in] LDAF
*> \verbatim
*> LDAF is INTEGER
*> The leading dimension of the array AF. LDAF >= max(1,N).
*> \endverbatim
*>
*> \param[in] IPIV
*> \verbatim
*> IPIV is INTEGER array, dimension (N)
*> The pivot indices from the factorization A = P*L*U
*> as computed by DGETRF; row i of the matrix was interchanged
*> with row IPIV(i).
*> \endverbatim
*>
*> \param[in] COLEQU
*> \verbatim
*> COLEQU is LOGICAL
*> If .TRUE. then column equilibration was done to A before calling
*> this routine. This is needed to compute the solution and error
*> bounds correctly.
*> \endverbatim
*>
*> \param[in] C
*> \verbatim
*> C is DOUBLE PRECISION array, dimension (N)
*> The column scale factors for A. If COLEQU = .FALSE., C
*> is not accessed. If C is input, each element of C should be a power
*> of the radix to ensure a reliable solution and error estimates.
*> Scaling by powers of the radix does not cause rounding errors unless
*> the result underflows or overflows. Rounding errors during scaling
*> lead to refining with a matrix that is not equivalent to the
*> input matrix, producing error estimates that may not be
*> reliable.
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
*> The right-hand-side matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*> LDB is INTEGER
*> The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] Y
*> \verbatim
*> Y is DOUBLE PRECISION array, dimension (LDY,NRHS)
*> On entry, the solution matrix X, as computed by DGETRS.
*> On exit, the improved solution matrix Y.
*> \endverbatim
*>
*> \param[in] LDY
*> \verbatim
*> LDY is INTEGER
*> The leading dimension of the array Y. LDY >= max(1,N).
*> \endverbatim
*>
*> \param[out] BERR_OUT
*> \verbatim
*> BERR_OUT is DOUBLE PRECISION array, dimension (NRHS)
*> On exit, BERR_OUT(j) contains the componentwise relative backward
*> error for right-hand-side j from the formula
*> max(i) ( abs(RES(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
*> where abs(Z) is the componentwise absolute value of the matrix
*> or vector Z. This is computed by DLA_LIN_BERR.
*> \endverbatim
*>
*> \param[in] N_NORMS
*> \verbatim
*> N_NORMS is INTEGER
*> Determines which error bounds to return (see ERRS_N
*> and ERRS_C).
*> If N_NORMS >= 1 return normwise error bounds.
*> If N_NORMS >= 2 return componentwise error bounds.
*> \endverbatim
*>
*> \param[in,out] ERRS_N
*> \verbatim
*> ERRS_N is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
*> For each right-hand side, this array contains information about
*> various error bounds and condition numbers corresponding to the
*> normwise relative error, which is defined as follows:
*>
*> Normwise relative error in the ith solution vector:
*> max_j (abs(XTRUE(j,i) - X(j,i)))
*> ------------------------------
*> max_j abs(X(j,i))
*>
*> The array is indexed by the type of error information as described
*> below. There currently are up to three pieces of information
*> returned.
*>
*> The first index in ERRS_N(i,:) corresponds to the ith
*> right-hand side.
*>
*> The second index in ERRS_N(:,err) contains the following
*> three fields:
*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
*> reciprocal condition number is less than the threshold
*> sqrt(n) * slamch('Epsilon').
*>
*> err = 2 "Guaranteed" error bound: The estimated forward error,
*> almost certainly within a factor of 10 of the true error
*> so long as the next entry is greater than the threshold
*> sqrt(n) * slamch('Epsilon'). This error bound should only
*> be trusted if the previous boolean is true.
*>
*> err = 3 Reciprocal condition number: Estimated normwise
*> reciprocal condition number. Compared with the threshold
*> sqrt(n) * slamch('Epsilon') to determine if the error
*> estimate is "guaranteed". These reciprocal condition
*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
*> appropriately scaled matrix Z.
*> Let Z = S*A, where S scales each row by a power of the
*> radix so all absolute row sums of Z are approximately 1.
*>
*> This subroutine is only responsible for setting the second field
*> above.
*> See Lapack Working Note 165 for further details and extra
*> cautions.
*> \endverbatim
*>
*> \param[in,out] ERRS_C
*> \verbatim
*> ERRS_C is DOUBLE PRECISION array, dimension (NRHS, N_ERR_BNDS)
*> For each right-hand side, this array contains information about
*> various error bounds and condition numbers corresponding to the
*> componentwise relative error, which is defined as follows:
*>
*> Componentwise relative error in the ith solution vector:
*> abs(XTRUE(j,i) - X(j,i))
*> max_j ----------------------
*> abs(X(j,i))
*>
*> The array is indexed by the right-hand side i (on which the
*> componentwise relative error depends), and the type of error
*> information as described below. There currently are up to three
*> pieces of information returned for each right-hand side. If
*> componentwise accuracy is not requested (PARAMS(3) = 0.0), then
*> ERRS_C is not accessed. If N_ERR_BNDS < 3, then at most
*> the first (:,N_ERR_BNDS) entries are returned.
*>
*> The first index in ERRS_C(i,:) corresponds to the ith
*> right-hand side.
*>
*> The second index in ERRS_C(:,err) contains the following
*> three fields:
*> err = 1 "Trust/don't trust" boolean. Trust the answer if the
*> reciprocal condition number is less than the threshold
*> sqrt(n) * slamch('Epsilon').
*>
*> err = 2 "Guaranteed" error bound: The estimated forward error,
*> almost certainly within a factor of 10 of the true error
*> so long as the next entry is greater than the threshold
*> sqrt(n) * slamch('Epsilon'). This error bound should only
*> be trusted if the previous boolean is true.
*>
*> err = 3 Reciprocal condition number: Estimated componentwise
*> reciprocal condition number. Compared with the threshold
*> sqrt(n) * slamch('Epsilon') to determine if the error
*> estimate is "guaranteed". These reciprocal condition
*> numbers are 1 / (norm(Z^{-1},inf) * norm(Z,inf)) for some
*> appropriately scaled matrix Z.
*> Let Z = S*(A*diag(x)), where x is the solution for the
*> current right-hand side and S scales each row of
*> A*diag(x) by a power of the radix so all absolute row
*> sums of Z are approximately 1.
*>
*> This subroutine is only responsible for setting the second field
*> above.
*> See Lapack Working Note 165 for further details and extra
*> cautions.
*> \endverbatim
*>
*> \param[in] RES
*> \verbatim
*> RES is DOUBLE PRECISION array, dimension (N)
*> Workspace to hold the intermediate residual.
*> \endverbatim
*>
*> \param[in] AYB
*> \verbatim
*> AYB is DOUBLE PRECISION array, dimension (N)
*> Workspace. This can be the same workspace passed for Y_TAIL.
*> \endverbatim
*>
*> \param[in] DY
*> \verbatim
*> DY is DOUBLE PRECISION array, dimension (N)
*> Workspace to hold the intermediate solution.
*> \endverbatim
*>
*> \param[in] Y_TAIL
*> \verbatim
*> Y_TAIL is DOUBLE PRECISION array, dimension (N)
*> Workspace to hold the trailing bits of the intermediate solution.
*> \endverbatim
*>
*> \param[in] RCOND
*> \verbatim
*> RCOND is DOUBLE PRECISION
*> Reciprocal scaled condition number. This is an estimate of the
*> reciprocal Skeel condition number of the matrix A after
*> equilibration (if done). If this is less than the machine
*> precision (in particular, if it is zero), the matrix is singular
*> to working precision. Note that the error may still be small even
*> if this number is very small and the matrix appears ill-
*> conditioned.
*> \endverbatim
*>
*> \param[in] ITHRESH
*> \verbatim
*> ITHRESH is INTEGER
*> The maximum number of residual computations allowed for
*> refinement. The default is 10. For 'aggressive' set to 100 to
*> permit convergence using approximate factorizations or
*> factorizations other than LU. If the factorization uses a
*> technique other than Gaussian elimination, the guarantees in
*> ERRS_N and ERRS_C may no longer be trustworthy.
*> \endverbatim
*>
*> \param[in] RTHRESH
*> \verbatim
*> RTHRESH is DOUBLE PRECISION
*> Determines when to stop refinement if the error estimate stops
*> decreasing. Refinement will stop when the next solution no longer
*> satisfies norm(dx_{i+1}) < RTHRESH * norm(dx_i) where norm(Z) is
*> the infinity norm of Z. RTHRESH satisfies 0 < RTHRESH <= 1. The
*> default value is 0.5. For 'aggressive' set to 0.9 to permit
*> convergence on extremely ill-conditioned matrices. See LAWN 165
*> for more details.
*> \endverbatim
*>
*> \param[in] DZ_UB
*> \verbatim
*> DZ_UB is DOUBLE PRECISION
*> Determines when to start considering componentwise convergence.
*> Componentwise convergence is only considered after each component
*> of the solution Y is stable, which we define as the relative
*> change in each component being less than DZ_UB. The default value
*> is 0.25, requiring the first bit to be stable. See LAWN 165 for
*> more details.
*> \endverbatim
*>
*> \param[in] IGNORE_CWISE
*> \verbatim
*> IGNORE_CWISE is LOGICAL
*> If .TRUE. then ignore componentwise convergence. Default value
*> is .FALSE..
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*> INFO is INTEGER
*> = 0: Successful exit.
*> < 0: if INFO = -i, the ith argument to DGETRS had an illegal
*> value
*> \endverbatim
*
* Authors:
* ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup la_gerfsx_extended
*
* =====================================================================
SUBROUTINE DLA_GERFSX_EXTENDED( PREC_TYPE, TRANS_TYPE, N, NRHS,
$ A,
$ LDA, AF, LDAF, IPIV, COLEQU, C, B,
$ LDB, Y, LDY, BERR_OUT, N_NORMS,
$ ERRS_N, ERRS_C, RES, AYB, DY,
$ Y_TAIL, RCOND, ITHRESH, RTHRESH,
$ DZ_UB, IGNORE_CWISE, INFO )
*
* -- LAPACK computational routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, LDAF, LDB, LDY, N, NRHS, PREC_TYPE,
$ TRANS_TYPE, N_NORMS, ITHRESH
LOGICAL COLEQU, IGNORE_CWISE
DOUBLE PRECISION RTHRESH, DZ_UB
* ..
* .. Array Arguments ..
INTEGER IPIV( * )
DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
$ Y( LDY, * ), RES( * ), DY( * ), Y_TAIL( * )
DOUBLE PRECISION C( * ), AYB( * ), RCOND, BERR_OUT( * ),
$ ERRS_N( NRHS, * ), ERRS_C( NRHS, * )
* ..
*
* =====================================================================
*
* .. Local Scalars ..
CHARACTER TRANS
INTEGER CNT, I, J, X_STATE, Z_STATE, Y_PREC_STATE
DOUBLE PRECISION YK, DYK, YMIN, NORMY, NORMX, NORMDX, DXRAT,
$ DZRAT, PREVNORMDX, PREV_DZ_Z, DXRATMAX,
$ DZRATMAX, DX_X, DZ_Z, FINAL_DX_X, FINAL_DZ_Z,
$ EPS, HUGEVAL, INCR_THRESH
LOGICAL INCR_PREC
* ..
* .. Parameters ..
INTEGER UNSTABLE_STATE, WORKING_STATE, CONV_STATE,
$ NOPROG_STATE, BASE_RESIDUAL, EXTRA_RESIDUAL,
$ EXTRA_Y
PARAMETER ( UNSTABLE_STATE = 0, WORKING_STATE = 1,
$ CONV_STATE = 2, NOPROG_STATE = 3 )
PARAMETER ( BASE_RESIDUAL = 0, EXTRA_RESIDUAL = 1,
$ EXTRA_Y = 2 )
INTEGER FINAL_NRM_ERR_I, FINAL_CMP_ERR_I, BERR_I
INTEGER RCOND_I, NRM_RCOND_I, NRM_ERR_I, CMP_RCOND_I
INTEGER CMP_ERR_I, PIV_GROWTH_I
PARAMETER ( FINAL_NRM_ERR_I = 1, FINAL_CMP_ERR_I = 2,
$ BERR_I = 3 )
PARAMETER ( RCOND_I = 4, NRM_RCOND_I = 5, NRM_ERR_I = 6 )
PARAMETER ( CMP_RCOND_I = 7, CMP_ERR_I = 8,
$ PIV_GROWTH_I = 9 )
INTEGER LA_LINRX_ITREF_I, LA_LINRX_ITHRESH_I,
$ LA_LINRX_CWISE_I
PARAMETER ( LA_LINRX_ITREF_I = 1,
$ LA_LINRX_ITHRESH_I = 2 )
PARAMETER ( LA_LINRX_CWISE_I = 3 )
INTEGER LA_LINRX_TRUST_I, LA_LINRX_ERR_I,
$ LA_LINRX_RCOND_I
PARAMETER ( LA_LINRX_TRUST_I = 1, LA_LINRX_ERR_I = 2 )
PARAMETER ( LA_LINRX_RCOND_I = 3 )
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DCOPY, DGETRS, DGEMV,
$ BLAS_DGEMV_X,
$ BLAS_DGEMV2_X, DLA_GEAMV, DLA_WWADDW, DLAMCH,
$ CHLA_TRANSTYPE, DLA_LIN_BERR
DOUBLE PRECISION DLAMCH
CHARACTER CHLA_TRANSTYPE
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
* ..
* .. Executable Statements ..
*
IF ( INFO.NE.0 ) RETURN
TRANS = CHLA_TRANSTYPE(TRANS_TYPE)
EPS = DLAMCH( 'Epsilon' )
HUGEVAL = DLAMCH( 'Overflow' )
* Force HUGEVAL to Inf
HUGEVAL = HUGEVAL * HUGEVAL
* Using HUGEVAL may lead to spurious underflows.
INCR_THRESH = DBLE( N ) * EPS
*
DO J = 1, NRHS
Y_PREC_STATE = EXTRA_RESIDUAL
IF ( Y_PREC_STATE .EQ. EXTRA_Y ) THEN
DO I = 1, N
Y_TAIL( I ) = 0.0D+0
END DO
END IF
DXRAT = 0.0D+0
DXRATMAX = 0.0D+0
DZRAT = 0.0D+0
DZRATMAX = 0.0D+0
FINAL_DX_X = HUGEVAL
FINAL_DZ_Z = HUGEVAL
PREVNORMDX = HUGEVAL
PREV_DZ_Z = HUGEVAL
DZ_Z = HUGEVAL
DX_X = HUGEVAL
X_STATE = WORKING_STATE
Z_STATE = UNSTABLE_STATE
INCR_PREC = .FALSE.
DO CNT = 1, ITHRESH
*
* Compute residual RES = B_s - op(A_s) * Y,
* op(A) = A, A**T, or A**H depending on TRANS (and type).
*
CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
IF ( Y_PREC_STATE .EQ. BASE_RESIDUAL ) THEN
CALL DGEMV( TRANS, N, N, -1.0D+0, A, LDA, Y( 1, J ),
$ 1,
$ 1.0D+0, RES, 1 )
ELSE IF ( Y_PREC_STATE .EQ. EXTRA_RESIDUAL ) THEN
CALL BLAS_DGEMV_X( TRANS_TYPE, N, N, -1.0D+0, A, LDA,
$ Y( 1, J ), 1, 1.0D+0, RES, 1, PREC_TYPE )
ELSE
CALL BLAS_DGEMV2_X( TRANS_TYPE, N, N, -1.0D+0, A, LDA,
$ Y( 1, J ), Y_TAIL, 1, 1.0D+0, RES, 1, PREC_TYPE )
END IF
! XXX: RES is no longer needed.
CALL DCOPY( N, RES, 1, DY, 1 )
CALL DGETRS( TRANS, N, 1, AF, LDAF, IPIV, DY, N, INFO )
*
* Calculate relative changes DX_X, DZ_Z and ratios DXRAT, DZRAT.
*
NORMX = 0.0D+0
NORMY = 0.0D+0
NORMDX = 0.0D+0
DZ_Z = 0.0D+0
YMIN = HUGEVAL
*
DO I = 1, N
YK = ABS( Y( I, J ) )
DYK = ABS( DY( I ) )
IF ( YK .NE. 0.0D+0 ) THEN
DZ_Z = MAX( DZ_Z, DYK / YK )
ELSE IF ( DYK .NE. 0.0D+0 ) THEN
DZ_Z = HUGEVAL
END IF
YMIN = MIN( YMIN, YK )
NORMY = MAX( NORMY, YK )
IF ( COLEQU ) THEN
NORMX = MAX( NORMX, YK * C( I ) )
NORMDX = MAX( NORMDX, DYK * C( I ) )
ELSE
NORMX = NORMY
NORMDX = MAX( NORMDX, DYK )
END IF
END DO
IF ( NORMX .NE. 0.0D+0 ) THEN
DX_X = NORMDX / NORMX
ELSE IF ( NORMDX .EQ. 0.0D+0 ) THEN
DX_X = 0.0D+0
ELSE
DX_X = HUGEVAL
END IF
DXRAT = NORMDX / PREVNORMDX
DZRAT = DZ_Z / PREV_DZ_Z
*
* Check termination criteria
*
IF (.NOT.IGNORE_CWISE
$ .AND. YMIN*RCOND .LT. INCR_THRESH*NORMY
$ .AND. Y_PREC_STATE .LT. EXTRA_Y)
$ INCR_PREC = .TRUE.
IF ( X_STATE .EQ. NOPROG_STATE .AND. DXRAT .LE. RTHRESH )
$ X_STATE = WORKING_STATE
IF ( X_STATE .EQ. WORKING_STATE ) THEN
IF ( DX_X .LE. EPS ) THEN
X_STATE = CONV_STATE
ELSE IF ( DXRAT .GT. RTHRESH ) THEN
IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
INCR_PREC = .TRUE.
ELSE
X_STATE = NOPROG_STATE
END IF
ELSE
IF ( DXRAT .GT. DXRATMAX ) DXRATMAX = DXRAT
END IF
IF ( X_STATE .GT. WORKING_STATE ) FINAL_DX_X = DX_X
END IF
IF ( Z_STATE .EQ. UNSTABLE_STATE .AND. DZ_Z .LE. DZ_UB )
$ Z_STATE = WORKING_STATE
IF ( Z_STATE .EQ. NOPROG_STATE .AND. DZRAT .LE. RTHRESH )
$ Z_STATE = WORKING_STATE
IF ( Z_STATE .EQ. WORKING_STATE ) THEN
IF ( DZ_Z .LE. EPS ) THEN
Z_STATE = CONV_STATE
ELSE IF ( DZ_Z .GT. DZ_UB ) THEN
Z_STATE = UNSTABLE_STATE
DZRATMAX = 0.0D+0
FINAL_DZ_Z = HUGEVAL
ELSE IF ( DZRAT .GT. RTHRESH ) THEN
IF ( Y_PREC_STATE .NE. EXTRA_Y ) THEN
INCR_PREC = .TRUE.
ELSE
Z_STATE = NOPROG_STATE
END IF
ELSE
IF ( DZRAT .GT. DZRATMAX ) DZRATMAX = DZRAT
END IF
IF ( Z_STATE .GT. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
END IF
*
* Exit if both normwise and componentwise stopped working,
* but if componentwise is unstable, let it go at least two
* iterations.
*
IF ( X_STATE.NE.WORKING_STATE ) THEN
IF ( IGNORE_CWISE) GOTO 666
IF ( Z_STATE.EQ.NOPROG_STATE .OR. Z_STATE.EQ.CONV_STATE )
$ GOTO 666
IF ( Z_STATE.EQ.UNSTABLE_STATE .AND. CNT.GT.1 ) GOTO 666
END IF
IF ( INCR_PREC ) THEN
INCR_PREC = .FALSE.
Y_PREC_STATE = Y_PREC_STATE + 1
DO I = 1, N
Y_TAIL( I ) = 0.0D+0
END DO
END IF
PREVNORMDX = NORMDX
PREV_DZ_Z = DZ_Z
*
* Update solution.
*
IF ( Y_PREC_STATE .LT. EXTRA_Y ) THEN
CALL DAXPY( N, 1.0D+0, DY, 1, Y( 1, J ), 1 )
ELSE
CALL DLA_WWADDW( N, Y( 1, J ), Y_TAIL, DY )
END IF
END DO
* Target of "IF (Z_STOP .AND. X_STOP)". Sun's f77 won't EXIT.
666 CONTINUE
*
* Set final_* when cnt hits ithresh.
*
IF ( X_STATE .EQ. WORKING_STATE ) FINAL_DX_X = DX_X
IF ( Z_STATE .EQ. WORKING_STATE ) FINAL_DZ_Z = DZ_Z
*
* Compute error bounds
*
IF (N_NORMS .GE. 1) THEN
ERRS_N( J, LA_LINRX_ERR_I ) = FINAL_DX_X / (1 - DXRATMAX)
END IF
IF ( N_NORMS .GE. 2 ) THEN
ERRS_C( J, LA_LINRX_ERR_I ) = FINAL_DZ_Z / (1 - DZRATMAX)
END IF
*
* Compute componentwise relative backward error from formula
* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) )
* where abs(Z) is the componentwise absolute value of the matrix
* or vector Z.
*
* Compute residual RES = B_s - op(A_s) * Y,
* op(A) = A, A**T, or A**H depending on TRANS (and type).
*
CALL DCOPY( N, B( 1, J ), 1, RES, 1 )
CALL DGEMV( TRANS, N, N, -1.0D+0, A, LDA, Y(1,J), 1, 1.0D+0,
$ RES, 1 )
DO I = 1, N
AYB( I ) = ABS( B( I, J ) )
END DO
*
* Compute abs(op(A_s))*abs(Y) + abs(B_s).
*
CALL DLA_GEAMV ( TRANS_TYPE, N, N, 1.0D+0,
$ A, LDA, Y(1, J), 1, 1.0D+0, AYB, 1 )
CALL DLA_LIN_BERR ( N, N, 1, RES, AYB, BERR_OUT( J ) )
*
* End of loop for each RHS.
*
END DO
*
RETURN
*
* End of DLA_GERFSX_EXTENDED
*
END
|
