snrm2.c
Go to the documentation of this file.
00001 /* Euclidean norm of the n-vector stored in sx()
00002    with storage -- Simplified Version
00003    Ver.1.0, May,26,1988.                         */
00004 
00005 #include "arith.h"
00006 #define CUTLO 4.441e-16
00007 #define CUTHI 1.304e19
00008 #define ZERO 0.0e0
00009 #define ONE 1.0e0
00010 #define SXI sx[isx+i][jsx]
00011 
00012 REAL snrm2(n,sx,isx,jsx,incx)
00013 int n,isx,jsx,incx;
00014 MATRIX sx;
00015 {
00016     int i,j,next,nn;
00017     REAL hitest,sum,xmax,snrm;
00018 
00019     if(n <=0)
00020         return(0.0);
00021 
00022     next=30;
00023     sum= ZERO;
00024     nn= n*incx;
00025 
00026     i=0;
00027 
00028   label20:
00029     switch(next){
00030       case 30:
00031         goto label30;
00032         break;
00033 
00034       case 50:
00035         goto label50;
00036         break;
00037 
00038       case 70:
00039         goto label70;
00040         break;
00041 
00042       case 110:
00043         goto label110;
00044         break;
00045     }
00046 
00047   label30:
00048     if(fabs(SXI) > CUTLO)
00049         goto label85;
00050     next=50;
00051     xmax=ZERO;
00052 
00053     /* phase 1. sum is zero */
00054 
00055   label50:
00056     if(SXI == ZERO)
00057         goto label200;
00058     if(fabs(SXI) > CUTLO)
00059         goto label85;
00060 
00061     /* prepare for pahse 2. */
00062 
00063     next=70;
00064     goto label105;
00065 
00066     /* prepare for phase 4. */
00067 
00068   label100:
00069     i=j;
00070     next=110;
00071     sum=(sum/SXI)/SXI;
00072 
00073   label105:
00074     xmax=fabs(SXI);
00075     goto label115;
00076 
00077     /* phase 2. sum is small.
00078                 scale to avoid destructive underflow. */
00079 
00080   label70:
00081     if(fabs(SXI) > CUTLO)
00082         goto label75;
00083 
00084     /* common code for phases 2 and 4.
00085        in phase 4 sum is large.  scale to avoid overflow. */
00086 
00087   label110:
00088     if(fabs(SXI) <= xmax)
00089         goto label115;
00090     sum= ONE + sum*(xmax/SXI)*(xmax/SXI);
00091     xmax=fabs(SXI);
00092     goto label200;
00093 
00094   label115:
00095     sum= sum+(SXI/xmax)*(SXI/xmax);
00096     goto label200;
00097 
00098     /* prepare for phase 3. */
00099 
00100   label75:
00101     sum= (sum*xmax)*xmax;
00102 
00103     /* for real or D.P. set hitest= CUTHI/n */
00104 
00105   label85:
00106     hitest= CUTHI/(double)n;
00107 
00108     /* phase 3. sum is mid-range. no scaling. */
00109 
00110     for(j=i; j<nn; j += incx){
00111         if(fabs(SXI) >= hitest)
00112             goto label100;
00113         sum= sum+ sx[isx+j][jsx]*sx[isx+j][jsx];
00114     }
00115 
00116     snrm= sqrt(sum);
00117     return(snrm);
00118 
00119   label200:
00120     i= i+incx;
00121     if(i < nn)
00122         goto label20;
00123 
00124     /* end of main loop
00125        computer square root and adjust for scaling. */
00126 
00127     snrm= xmax*sqrt(sum);
00128     return(snrm);
00129 }


euslisp
Author(s): Toshihiro Matsui
autogenerated on Thu Mar 9 2017 04:57:50