00001 /* dlarfp.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 dlarfp_(integer *n, doublereal *alpha, doublereal *x, 00017 integer *incx, doublereal *tau) 00018 { 00019 /* System generated locals */ 00020 integer i__1; 00021 doublereal d__1; 00022 00023 /* Builtin functions */ 00024 double d_sign(doublereal *, doublereal *); 00025 00026 /* Local variables */ 00027 integer j, knt; 00028 doublereal beta; 00029 extern doublereal dnrm2_(integer *, doublereal *, integer *); 00030 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 00031 integer *); 00032 doublereal xnorm; 00033 extern doublereal dlapy2_(doublereal *, doublereal *), dlamch_(char *); 00034 doublereal safmin, rsafmn; 00035 00036 00037 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00038 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00039 /* November 2006 */ 00040 00041 /* .. Scalar Arguments .. */ 00042 /* .. */ 00043 /* .. Array Arguments .. */ 00044 /* .. */ 00045 00046 /* Purpose */ 00047 /* ======= */ 00048 00049 /* DLARFP generates a real elementary reflector H of order n, such */ 00050 /* that */ 00051 00052 /* H * ( alpha ) = ( beta ), H' * H = I. */ 00053 /* ( x ) ( 0 ) */ 00054 00055 /* where alpha and beta are scalars, beta is non-negative, and x is */ 00056 /* an (n-1)-element real vector. H is represented in the form */ 00057 00058 /* H = I - tau * ( 1 ) * ( 1 v' ) , */ 00059 /* ( v ) */ 00060 00061 /* where tau is a real scalar and v is a real (n-1)-element */ 00062 /* vector. */ 00063 00064 /* If the elements of x are all zero, then tau = 0 and H is taken to be */ 00065 /* the unit matrix. */ 00066 00067 /* Otherwise 1 <= tau <= 2. */ 00068 00069 /* Arguments */ 00070 /* ========= */ 00071 00072 /* N (input) INTEGER */ 00073 /* The order of the elementary reflector. */ 00074 00075 /* ALPHA (input/output) DOUBLE PRECISION */ 00076 /* On entry, the value alpha. */ 00077 /* On exit, it is overwritten with the value beta. */ 00078 00079 /* X (input/output) DOUBLE PRECISION array, dimension */ 00080 /* (1+(N-2)*abs(INCX)) */ 00081 /* On entry, the vector x. */ 00082 /* On exit, it is overwritten with the vector v. */ 00083 00084 /* INCX (input) INTEGER */ 00085 /* The increment between elements of X. INCX > 0. */ 00086 00087 /* TAU (output) DOUBLE PRECISION */ 00088 /* The value tau. */ 00089 00090 /* ===================================================================== */ 00091 00092 /* .. Parameters .. */ 00093 /* .. */ 00094 /* .. Local Scalars .. */ 00095 /* .. */ 00096 /* .. External Functions .. */ 00097 /* .. */ 00098 /* .. Intrinsic Functions .. */ 00099 /* .. */ 00100 /* .. External Subroutines .. */ 00101 /* .. */ 00102 /* .. Executable Statements .. */ 00103 00104 /* Parameter adjustments */ 00105 --x; 00106 00107 /* Function Body */ 00108 if (*n <= 0) { 00109 *tau = 0.; 00110 return 0; 00111 } 00112 00113 i__1 = *n - 1; 00114 xnorm = dnrm2_(&i__1, &x[1], incx); 00115 00116 if (xnorm == 0.) { 00117 00118 /* H = [+/-1, 0; I], sign chosen so ALPHA >= 0 */ 00119 00120 if (*alpha >= 0.) { 00121 /* When TAU.eq.ZERO, the vector is special-cased to be */ 00122 /* all zeros in the application routines. We do not need */ 00123 /* to clear it. */ 00124 *tau = 0.; 00125 } else { 00126 /* However, the application routines rely on explicit */ 00127 /* zero checks when TAU.ne.ZERO, and we must clear X. */ 00128 *tau = 2.; 00129 i__1 = *n - 1; 00130 for (j = 1; j <= i__1; ++j) { 00131 x[(j - 1) * *incx + 1] = 0.; 00132 } 00133 *alpha = -(*alpha); 00134 } 00135 } else { 00136 00137 /* general case */ 00138 00139 d__1 = dlapy2_(alpha, &xnorm); 00140 beta = d_sign(&d__1, alpha); 00141 safmin = dlamch_("S") / dlamch_("E"); 00142 knt = 0; 00143 if (abs(beta) < safmin) { 00144 00145 /* XNORM, BETA may be inaccurate; scale X and recompute them */ 00146 00147 rsafmn = 1. / safmin; 00148 L10: 00149 ++knt; 00150 i__1 = *n - 1; 00151 dscal_(&i__1, &rsafmn, &x[1], incx); 00152 beta *= rsafmn; 00153 *alpha *= rsafmn; 00154 if (abs(beta) < safmin) { 00155 goto L10; 00156 } 00157 00158 /* New BETA is at most 1, at least SAFMIN */ 00159 00160 i__1 = *n - 1; 00161 xnorm = dnrm2_(&i__1, &x[1], incx); 00162 d__1 = dlapy2_(alpha, &xnorm); 00163 beta = d_sign(&d__1, alpha); 00164 } 00165 *alpha += beta; 00166 if (beta < 0.) { 00167 beta = -beta; 00168 *tau = -(*alpha) / beta; 00169 } else { 00170 *alpha = xnorm * (xnorm / *alpha); 00171 *tau = *alpha / beta; 00172 *alpha = -(*alpha); 00173 } 00174 i__1 = *n - 1; 00175 d__1 = 1. / *alpha; 00176 dscal_(&i__1, &d__1, &x[1], incx); 00177 00178 /* If BETA is subnormal, it may lose relative accuracy */ 00179 00180 i__1 = knt; 00181 for (j = 1; j <= i__1; ++j) { 00182 beta *= safmin; 00183 /* L20: */ 00184 } 00185 *alpha = beta; 00186 } 00187 00188 return 0; 00189 00190 /* End of DLARFP */ 00191 00192 } /* dlarfp_ */