00001 /* sla_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 sla_lin_berr__(integer *n, integer *nz, integer *nrhs, 00017 real *res, real *ayb, real *berr) 00018 { 00019 /* System generated locals */ 00020 integer ayb_dim1, ayb_offset, res_dim1, res_offset, i__1, i__2; 00021 real r__1; 00022 00023 /* Local variables */ 00024 integer i__, j; 00025 real tmp, safe1; 00026 extern doublereal slamch_(char *); 00027 00028 00029 /* -- LAPACK routine (version 3.2.1) -- */ 00030 /* -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */ 00031 /* -- Jason Riedy of Univ. of California Berkeley. -- */ 00032 /* -- April 2009 -- */ 00033 00034 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ 00035 /* -- Univ. of California Berkeley and NAG Ltd. -- */ 00036 00037 /* .. */ 00038 /* .. Scalar Arguments .. */ 00039 /* .. */ 00040 /* .. Array Arguments .. */ 00041 /* .. */ 00042 00043 /* Purpose */ 00044 /* ======= */ 00045 00046 /* SLA_LIN_BERR computes componentwise relative backward error from */ 00047 /* the formula */ 00048 /* max(i) ( abs(R(i)) / ( abs(op(A_s))*abs(Y) + abs(B_s) )(i) ) */ 00049 /* where abs(Z) is the componentwise absolute value of the matrix */ 00050 /* or vector Z. */ 00051 00052 /* Arguments */ 00053 /* ========== */ 00054 00055 /* N (input) INTEGER */ 00056 /* The number of linear equations, i.e., the order of the */ 00057 /* matrix A. N >= 0. */ 00058 00059 /* NZ (input) INTEGER */ 00060 /* We add (NZ+1)*SLAMCH( 'Safe minimum' ) to R(i) in the numerator to */ 00061 /* guard against spuriously zero residuals. Default value is N. */ 00062 00063 /* NRHS (input) INTEGER */ 00064 /* The number of right hand sides, i.e., the number of columns */ 00065 /* of the matrices AYB, RES, and BERR. NRHS >= 0. */ 00066 00067 /* RES (input) REAL array, dimension (N,NRHS) */ 00068 /* The residual matrix, i.e., the matrix R in the relative backward */ 00069 /* error formula above. */ 00070 00071 /* AYB (input) REAL array, dimension (N, NRHS) */ 00072 /* The denominator in the relative backward error formula above, i.e., */ 00073 /* the matrix abs(op(A_s))*abs(Y) + abs(B_s). The matrices A, Y, and B */ 00074 /* are from iterative refinement (see sla_gerfsx_extended.f). */ 00075 00076 /* RES (output) REAL array, dimension (NRHS) */ 00077 /* The componentwise relative backward error from the formula above. */ 00078 00079 /* ===================================================================== */ 00080 00081 /* .. Local Scalars .. */ 00082 /* .. */ 00083 /* .. Intrinsic Functions .. */ 00084 /* .. */ 00085 /* .. External Functions .. */ 00086 /* .. */ 00087 /* .. Executable Statements .. */ 00088 00089 /* Adding SAFE1 to the numerator guards against spuriously zero */ 00090 /* residuals. A similar safeguard is in the SLA_yyAMV routine used */ 00091 /* to compute AYB. */ 00092 00093 /* Parameter adjustments */ 00094 --berr; 00095 ayb_dim1 = *n; 00096 ayb_offset = 1 + ayb_dim1; 00097 ayb -= ayb_offset; 00098 res_dim1 = *n; 00099 res_offset = 1 + res_dim1; 00100 res -= res_offset; 00101 00102 /* Function Body */ 00103 safe1 = slamch_("Safe minimum"); 00104 safe1 = (*nz + 1) * safe1; 00105 i__1 = *nrhs; 00106 for (j = 1; j <= i__1; ++j) { 00107 berr[j] = 0.f; 00108 i__2 = *n; 00109 for (i__ = 1; i__ <= i__2; ++i__) { 00110 if (ayb[i__ + j * ayb_dim1] != 0.f) { 00111 tmp = (safe1 + (r__1 = res[i__ + j * res_dim1], dabs(r__1))) / 00112 ayb[i__ + j * ayb_dim1]; 00113 /* Computing MAX */ 00114 r__1 = berr[j]; 00115 berr[j] = dmax(r__1,tmp); 00116 } 00117 00118 /* If AYB is exactly 0.0 (and if computed by SLA_yyAMV), then we know */ 00119 /* the true residual also must be exactly 0.0. */ 00120 00121 } 00122 } 00123 return 0; 00124 } /* sla_lin_berr__ */