00001 /* zla_lin_berr.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Subroutine */ int zla_lin_berr__(integer *n, integer *nz, integer *nrhs, 00017 doublecomplex *res, doublereal *ayb, doublereal *berr) 00018 { 00019 /* System generated locals */ 00020 integer ayb_dim1, ayb_offset, res_dim1, res_offset, i__1, i__2, i__3, 00021 i__4; 00022 doublereal d__1, d__2, d__3; 00023 doublecomplex z__1, z__2, z__3; 00024 00025 /* Builtin functions */ 00026 double d_imag(doublecomplex *); 00027 00028 /* Local variables */ 00029 integer i__, j; 00030 doublereal tmp, safe1; 00031 extern doublereal dlamch_(char *); 00032 00033 00034 /* -- LAPACK routine (version 3.2.1) -- */ 00035 /* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */ 00036 /* -- Jason Riedy of Univ. of California Berkeley. -- */ 00037 /* -- April 2009 -- */ 00038 00039 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ 00040 /* -- Univ. of California Berkeley and NAG Ltd. -- */ 00041 00042 /* .. */ 00043 /* .. Scalar Arguments .. */ 00044 /* .. */ 00045 /* .. Array Arguments .. */ 00046 /* .. */ 00047 00048 /* Purpose */ 00049 /* ======= */ 00050 00051 /* ZLA_LIN_BERR computes componentwise relative backward error from */ 00052 /* the formula */ 00053 /* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) */ 00054 /* where abs(Z) is the componentwise absolute value of the matrix */ 00055 /* or vector Z. */ 00056 00057 /* N (input) INTEGER */ 00058 /* The number of linear equations, i.e., the order of the */ 00059 /* matrix A. N >= 0. */ 00060 00061 /* NZ (input) INTEGER */ 00062 /* We add (NZ+1)*SLAMCH( 'Safe minimum' ) to R(i) in the numerator to */ 00063 /* guard against spuriously zero residuals. Default value is N. */ 00064 00065 /* NRHS (input) INTEGER */ 00066 /* The number of right hand sides, i.e., the number of columns */ 00067 /* of the matrices AYB, RES, and BERR. NRHS >= 0. */ 00068 00069 /* RES (input) DOUBLE PRECISION array, dimension (N,NRHS) */ 00070 /* The residual matrix, i.e., the matrix R in the relative backward */ 00071 /* error formula above. */ 00072 00073 /* AYB (input) DOUBLE PRECISION array, dimension (N, NRHS) */ 00074 /* The denominator in the relative backward error formula above, i.e., */ 00075 /* the matrix abs(op(A_s))*abs(Y) + abs(B_s). The matrices A, Y, and B */ 00076 /* are from iterative refinement (see zla_gerfsx_extended.f). */ 00077 00078 /* RES (output) COMPLEX*16 array, dimension (NRHS) */ 00079 /* The componentwise relative backward error from the formula above. */ 00080 00081 /* ===================================================================== */ 00082 00083 /* .. Local Scalars .. */ 00084 /* .. */ 00085 /* .. Intrinsic Functions .. */ 00086 /* .. */ 00087 /* .. External Functions .. */ 00088 /* .. */ 00089 /* .. Statement Functions .. */ 00090 /* .. */ 00091 /* .. Statement Function Definitions .. */ 00092 /* .. */ 00093 /* .. Executable Statements .. */ 00094 00095 /* Adding SAFE1 to the numerator guards against spuriously zero */ 00096 /* residuals. A similar safeguard is in the CLA_yyAMV routine used */ 00097 /* to compute AYB. */ 00098 00099 /* Parameter adjustments */ 00100 --berr; 00101 ayb_dim1 = *n; 00102 ayb_offset = 1 + ayb_dim1; 00103 ayb -= ayb_offset; 00104 res_dim1 = *n; 00105 res_offset = 1 + res_dim1; 00106 res -= res_offset; 00107 00108 /* Function Body */ 00109 safe1 = dlamch_("Safe minimum"); 00110 safe1 = (*nz + 1) * safe1; 00111 i__1 = *nrhs; 00112 for (j = 1; j <= i__1; ++j) { 00113 berr[j] = 0.; 00114 i__2 = *n; 00115 for (i__ = 1; i__ <= i__2; ++i__) { 00116 if (ayb[i__ + j * ayb_dim1] != 0.) { 00117 i__3 = i__ + j * res_dim1; 00118 d__3 = (d__1 = res[i__3].r, abs(d__1)) + (d__2 = d_imag(&res[ 00119 i__ + j * res_dim1]), abs(d__2)); 00120 z__3.r = d__3, z__3.i = 0.; 00121 z__2.r = safe1 + z__3.r, z__2.i = z__3.i; 00122 i__4 = i__ + j * ayb_dim1; 00123 z__1.r = z__2.r / ayb[i__4], z__1.i = z__2.i / ayb[i__4]; 00124 tmp = z__1.r; 00125 /* Computing MAX */ 00126 d__1 = berr[j]; 00127 berr[j] = max(d__1,tmp); 00128 } 00129 00130 /* If AYB is exactly 0.0 (and if computed by CLA_yyAMV), then we know */ 00131 /* the true residual also must be exactly 0.0. */ 00132 00133 } 00134 } 00135 return 0; 00136 } /* zla_lin_berr__ */