/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:06 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dnrm2.h" #include #include double /*FUNCTION*/ dnrm2( long n, double x[], long incx) { long int _d_l, _d_m, _do0, _do1, _do2, _do3, _l0, i, id, iemn, iemx, nn; double dnrm2_v, sum, sumx, xa; static double asmall, dbig, dsmall, fbig, fsmall, tbig, tsmall; static double abig = 0.e0; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const X = &x[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1998-05-11 DNRM2 Krogh Minor changes for conversion to C. *>> 1996-08-29 DNRM2 Krogh Coded an entirely different algorithm *>> .... *>> 1985-08-02 DNRM2 Lawson Initial code. * * EUCLIDEAN NORM OF THE N-VECTOR STORED IN X() WITH STORAGE * INCREMENT INCX . * IF N .LE. 0 RETURN WITH RESULT = 0. * IF N .GE. 1 THEN INCX MUST BE .GE. 1 * * C.L.LAWSON, 1978 JAN 08 * * New algorithm avoids underflow as well as overflow, and avoids divides * as much as possible. F. Krogh, August 29, 1996. * * ************************* Variable Definitions *********************** * * Let b be the base for floating point, b**E be the smallest integer * power of b that overflows, b**e be the smallest integer power of b * that does not underflow, and d be the number of base b digits in a * a floating point number. The descriptions below use b, e, E, and d. * * ABIG = b ** ((E+5) / 4)). This is used as a factor to get the * final answer when combining SUMX from big XA, with SUM. * ASMALL = b ** ((e+d-E+5) / 4)). This is used as a factor to get the * final answer when combining SUMX from small XA, with SUM. * DBIG = b ** (2*((E+5)/4)). This is used as a factor to get the * final answer when only SUMX formed from big XA's is needed. * DSMALL = b ** (2*((e+d-E+5)/4)). This is used as a factor to get the * final answer when only SUMX formed from big XA's is needed. * FBIG = b ** (-2*((E+5)/4)). This is used as a multiplier when * accumulating SUMX for big XA. * FSMALL = b ** (-2*((e+d-E+5)/4)). This is used as a multiplier when * accumulating SUMX for small XA. * I Temporary index. * ID Number of base I1MACH(10) digits in Floating point number. * IEMN The minimum floating point exponent. * IEMX The maximum floating point exponent. * INCX Input, the increment (>0) between elements of X. * N Input, the number of elements in X. * NN N * INCX = last index processed. * SUM Place where sum of XA**2 is accumlated. * SUMX Place where scaled sum of XA**2 is accumlated. In first loop * this is for XA**2 that would likely underflow, and in second loop * this is for XA**2 that would likely overflow. * TBIG = b ** ((E - d - 1) / 2). If XA > TBIG, number is "big". Note * that for such XA's, (FBIG * XA) ** 2 should not underflow and the * accumlation overflows only if the final result would. * TSMALL = b ** ((e+1) / 2). If XA <= TSMALL, number is "small". Note * that for such XA's, (FSMALL * XA)**2 should not underflow and the * accumlation should not overlflow. * X Input, we are getting the L2 norm of X. * XA Contains base of floating point numbers when getting saved * parameters. Later contains abs(X(I)). * ------------------------------------------------------------------ *--D replaces "?": ?NRM2 * ------------------------------------------------------------------ */ /* Values of these saved parameters for D.P. IEEE arithmetic are: * ABIG= .2315841784746324E+78 DBIG= .5363123171977043E+155 * FBIG= .1864585182800050E-154 TBIG= .9989595361011182E+146 * ASMALL= .4887898181599363E-149 DSMALL= .2389154863368240E-298 * FSMALL= .4185580496821357E+299 TSMALL= .2983336292480080E-153 * Values of these saved parameters for S.P. IEEE arithmetic are: * ABIG= .8589935E+10 DBIG= .7378698E+20 * FBIG= .1355253E-19 TBIG= .2251800E+16 * ASMALL= .1387779E-16 DSMALL= .1925930E-33 * FSMALL= .5192297E+34 TSMALL= .2168404E-18 * */ /* ------------------------------------------------------------------ * * */ if (abig == 0.e0) { /*++ Code for (.N. == 'D') is active */ iemx = FLT_MAX_EXP; iemn = FLT_MIN_EXP; id = FLT_MANT_DIG; /*++ Code for (.N. == 'S') is inactive * IEMX = I1MACH(13) * IEMN = I1MACH(12) * ID = I1MACH(11) *++ END */ xa = (double)( FLT_RADIX ); abig = powi(xa,(iemx + 5)/4); dbig = SQ(abig); fbig = 1.e0/dbig; tbig = powi(xa,(iemx - id - 1)/2); asmall = powi(xa,(iemn + id - iemx + 5)/4); dsmall = SQ(asmall); fsmall = 1.e0/dsmall; tsmall = powi(xa,(iemn + 1)/2); } sum = 0.e0; if (n > 0) { nn = n*incx; sumx = 0.e0; /* Loop when no big number yet encountered. */ for (i = 1, _do0=DOCNT(i,nn,_do1 = incx); _do0 > 0; i += _do1, _do0--) { xa = fabs( X[i] ); if (xa < tsmall) { sumx += powi(fsmall*xa,2); } else { if (xa > tbig) goto L_200; sum += SQ(xa); } } if (sum != 0.e0) { if (sumx >= 1.e0) { if (sum < 1.e0) { sum = asmall*sqrt( fsmall*sum + dsmall*sumx ); goto L_400; } } sum = sqrt( sum ); } else { sum = dsmall*sqrt( sumx ); } goto L_400; L_200: sumx = 0.e0; /* Loop when we have at least one big number. */ for (i = i, _do2=DOCNT(i,nn,_do3 = incx); _do2 > 0; i += _do3, _do2--) { xa = fabs( X[i] ); if (xa > tsmall) { if (xa > tbig) { sumx += powi(fbig*xa,2); } else { sum += SQ(xa); } } } if ((sumx <= 1.e10) && (sum >= 1.e-10)) { sum = abig*sqrt( fbig*sum + dbig*sumx ); } else { sum = dbig*sqrt( sumx ); } } L_400: ; dnrm2_v = sum; return( dnrm2_v ); } /* end of function */