Go to the documentation of this file.00001
00002
00003
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
00054
00055 label50:
00056 if(SXI == ZERO)
00057 goto label200;
00058 if(fabs(SXI) > CUTLO)
00059 goto label85;
00060
00061
00062
00063 next=70;
00064 goto label105;
00065
00066
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
00078
00079
00080 label70:
00081 if(fabs(SXI) > CUTLO)
00082 goto label75;
00083
00084
00085
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
00099
00100 label75:
00101 sum= (sum*xmax)*xmax;
00102
00103
00104
00105 label85:
00106 hitest= CUTHI/(double)n;
00107
00108
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
00125
00126
00127 snrm= xmax*sqrt(sum);
00128 return(snrm);
00129 }