| /*************************************************************** |
| ***************************************************************** |
| ******************************************************************* |
| ***** ***** |
| ***** BLAS ***** |
| ***** Basic Linear Algebra Subroutines ***** |
| ***** Written in the C Programming Language. ***** |
| ***** ***** |
| ***** Functions include: ***** |
| ***** isamax, saxpy, saxpyx, scopy, sdot, snrm2, **** |
| ***** vexopy, vfill ***** |
| ***** ***** |
| ******************************************************************* |
| ***************************************************************** |
| ***************************************************************/ |
| #include <math.h> |
| #ifndef SMALLsp |
| #define SMALLsp 1.0e-45 /* Smallest (pos) binary float */ |
| #endif |
| #ifndef HUGEsp |
| #define HUGEsp 1.0e+38 /* Largest binary float */ |
| #endif |
| int isamax( n, sx, incx ) |
| float *sx; |
| int n, incx; |
| /* |
| PURPOSE |
| Finds the index of element having max. absolute value. |
| |
| INPUT |
| n Number of elements to check. |
| sx Vector to be checked. |
| incx Every incx-th element is checked. |
| |
| */ |
| { |
| float smax = 0.0e0; |
| int i, istmp = 0; |
| |
| #ifndef abs |
| #define abs(x) ((x)>0.0?(x):-(x)) |
| #endif |
| if( n <= 1 ) return( istmp ); |
| if( incx != 1 ) { |
| /* Code for increment not equal to 1. */ |
| if( incx < 0 ) sx = sx + ((-n+1)*incx + 1); |
| istmp = 0; |
| smax = abs( *sx ); |
| sx += incx; |
| for( i=1; i<n; i++, sx+=incx ) |
| if( abs( *sx ) > smax ) { |
| istmp = i; |
| smax = abs( *sx ); |
| } |
| return( istmp ); |
| } |
| /* Code for increment equal to 1. */ |
| istmp = 0; |
| smax = abs(*sx); |
| sx++; |
| for( i=1; i<n; i++, sx++ ) |
| if( abs( *sx ) > smax ) { |
| istmp = i; |
| smax = abs( *sx ); |
| } |
| return( istmp ); |
| } |
| |
| void saxpy( n, sa, sx, incx, sy, incy ) |
| float *sx, *sy, sa; |
| int n, incx, incy; |
| /* |
| PURPOSE |
| Vector times a scalar plus a vector. sy = sy + sa*sx. |
| |
| INPUT |
| n Number of elements to multiply. |
| sa Scalar to multiply by. |
| sx Pointer to float vector to scale. |
| incx Storage incrament for sx. |
| sy Pointer to float vector to add. |
| incy Storage incrament for sy. |
| |
| OUTPUT |
| sy sy = sy + sa*sx |
| */ |
| { |
| register int i; |
| |
| if( n<=0 || sa==0.0 ) return; |
| if( incx == incy ) { |
| if( incx == 1 ) { |
| /* Both increments = 1 */ |
| for( i=0; i<n; i++,sy++,sx++ ) |
| *sy += sa*(*sx); |
| return; |
| } |
| if( incx>0 ) { |
| /* Equal, positive, non-unit increments. */ |
| for( i=0; i<n; i++,sx+=incx,sy+=incx ) |
| *sy += sa*(*sx); |
| return; |
| } |
| } |
| /* Unequal or negative increments. */ |
| if( incx < 0 ) sx += ((-n+1)*incx + 1); |
| if( incy < 0 ) sy += ((-n+1)*incy + 1); |
| for( i=0; i<n; i++,sx+=incx,sy+=incy ) |
| *sy += sa*(*sx); |
| } |
| |
| void saxpyx( n, sa, sx, incx, sy, incy ) |
| float *sx, *sy, sa; |
| int n, incx, incy; |
| /* |
| PURPOSE |
| Vector times a scalar plus a vector. sx = sy + sa*sx. |
| |
| INPUT |
| n Number of elements to multiply. |
| sa Scalar to multiply by. |
| sx Pointer to float vector to scale. |
| incx Storage incrament for sx. |
| sy Pointer to float vector to add. |
| incy Storage incrament for sy. |
| |
| OUTPUT |
| sx sx = sy + sa*sx |
| */ |
| { |
| register i; |
| |
| if( n<=0 || sa==0.0 ) return; |
| if( incx == incy ) { |
| if( incx == 1 ) { |
| /* Both increments = 1 */ |
| for( i=0; i<n; i++, sx++, sy++ ) |
| *sx = *sy + sa*(*sx); |
| return; |
| } |
| if( incx>0 ) { |
| /* Equal, positive, non-unit increments. */ |
| for( i=0; i<n; i++, sx+=incx, sy+=incx) |
| *sx = *sy + sa*(*sx); |
| return; |
| } |
| } |
| /* Unequal or negative increments. */ |
| if( incx < 0 ) sx += ((-n+1)*incx + 1); |
| if( incy < 0 ) sy += ((-n+1)*incy + 1); |
| for( i=0; i<n; i++,sx+=incx,sy+=incy ) |
| *sx = *sy + sa*(*sx); |
| } |
| |
| void scopy( n, sx, incx, sy, incy ) |
| float *sx, *sy; |
| int n, incx, incy; |
| /* |
| PURPOSE |
| Copies vector sx into vector sy. |
| |
| INPUT |
| n Number of components to copy. |
| sx Source vector |
| incx Index increment for sx. |
| incy Index increment for sy. |
| |
| OUTPUT |
| sy Destination vector. |
| */ |
| { |
| register int i; |
| |
| if( n<1 ) return; |
| if( incx == incy ) { |
| if( incx == 1 ) { |
| /* Both increments = 1 */ |
| for( i=0; i<n; i++ ) |
| *(sy++) = *(sx++); |
| return; |
| } |
| if( incx > 0 ) { |
| /* Equal, positive, non-unit increments. */ |
| for( i=0; i<n; i++, sx+=incx, sy+=incx) |
| *sy = *sx; |
| return; |
| } |
| } |
| /* Non-equal or negative increments. */ |
| if( incx < 0 ) sx += ((-n+1)*incx + 1); |
| if( incy < 0 ) sy += ((-n+1)*incy + 1); |
| for( i=0; i<n; i++,sx+=incx,sy+=incy ) |
| (*sx) = (*sy); |
| return; |
| } |
| double sdot( n, sx, incx, sy, incy ) |
| float *sx, *sy; |
| int n, incx, incy; |
| /* |
| PURPOSE |
| Forms the dot product of a vector. |
| |
| INPUT |
| n Number of elements to sum. |
| sx Address of first element of x vector. |
| incx Incrament for the x vector. |
| sy Address of first element of y vector. |
| incy incrament for the y vector. |
| |
| OUPUT |
| sdot Dot product x and y. Double returned |
| due to `C' language features. |
| */ |
| { |
| register i; |
| float stemp = 0.0e0; |
| |
| if( n<1 ) return( stemp ); |
| if( incx == incy ) { |
| if( incx == 1 ) { |
| /* Both increments = 1 */ |
| for( i=0; i<n; i++, sx++, sy++ ) |
| stemp += (*sx)*(*sy); |
| return( stemp ); |
| } |
| if( incx>0 ) { |
| /* Equal, positive, non-unit increments. */ |
| for( i=0; i<n; i++, sx+=incx, sy+=incx) |
| stemp += (*sx)*(*sy); |
| return( stemp ); |
| } |
| } |
| /* Unequal or negative increments. */ |
| if( incx < 0 ) sx += ((-n+1)*incx + 1); |
| if( incy < 0 ) sy += ((-n+1)*incy + 1); |
| for( i=0; i<n; i++,sx+=incx,sy+=incy ) |
| stemp += (*sx)*(*sy); |
| return( stemp ); |
| } /* End of ---SDOT--- */ |
| |
| double snrm2( n, sx, incx ) |
| float *sx; |
| int n, incx; |
| /* |
| PURPOSE |
| Computes the Euclidean norm of sx while being |
| very careful of distructive underflow and overflow. |
| |
| INPUT |
| n Number of elements to use. |
| sx Address of first element of x vector. |
| incx Incrament for the x vector (>0). |
| |
| OUPUT |
| snrm2 Euclidean norm of sx. Returns double |
| due to `C' language features. |
| REMARKS |
| This algorithm proceeds in four steps. |
| 1) scan zero components. |
| 2) do phase 2 when component is near underflow. |
| */ |
| { |
| register int i; |
| int phase = 3; |
| double sum = 0.0e0, cutlo, cuthi, hitst, r1mach(); |
| float xmax; |
| |
| if( n<1 || incx<1 ) return( sum ); |
| |
| cutlo = sqrt( SMALLsp/r1mach() ); /* Calculate near underflow */ |
| cuthi = sqrt( HUGEsp ); /* Calculate near overflow */ |
| hitst = cuthi/(double) n; |
| i = 0; |
| |
| /* Zero Sum. */ |
| while( i < n && *sx == 0.0 ) { |
| i++; |
| sx += incx; |
| } |
| if( i>=n ) return( sum ); |
| |
| START: |
| if( abs( *sx ) > cutlo ) { |
| for( ; i<n; i++, sx+=incx ) { /* Loop over elements. */ |
| if( abs( *sx ) > hitst ) goto GOT_LARGE; |
| sum += (*sx) * (*sx); |
| } |
| sum = sqrt( sum ); |
| return( sum ); /* Sum completed normaly. */ |
| } |
| else { /* Small sum prepare for phase 2. */ |
| xmax = abs( *sx ); |
| sx += incx; |
| i++; |
| sum += 1.0; |
| for( ; i<n; i++, sx+=incx ) { |
| if( abs( *sx ) > cutlo ) { /* Got normal elem. Rescale and process. */ |
| sum = (sum*xmax)*xmax; |
| goto START; |
| } |
| if( abs( *sx ) > xmax ) { |
| sum = 1.0 + sum*(xmax /(*sx))*(xmax /(*sx)); |
| xmax = abs( *sx ); |
| continue; |
| } |
| sum += ((*sx)/xmax)*((*sx)/xmax); |
| } |
| return( (double)xmax*sqrt( sum ) ); |
| } /* End of small sum. */ |
| |
| GOT_LARGE: |
| sum = 1.0 + (sum/(*sx))/(*sx); /* Rescale and process. */ |
| xmax = abs( *sx ); |
| sx += incx; |
| i++; |
| for( ; i<n; i++, sx+=incx ) { |
| if( abs( *sx ) > xmax ) { |
| sum = 1.0 + sum*(xmax /(*sx))*(xmax /(*sx)); |
| xmax = abs( *sx ); |
| continue; |
| } |
| sum += ((*sx)/xmax)*((*sx)/xmax); |
| } |
| return( (double)xmax*sqrt( sum ) ); /* End of small sum. */ |
| } /* End of ---SDOT--- */ |
| |
| double r1mach() |
| /* --------------------------------------------------------------------- |
| This routine computes the unit roundoff for single precision |
| of the machine. This is defined as the smallest positive |
| machine number u such that 1.0 + u .ne. 1.0 |
| Returns a double due to `C' language features. |
| ----------------------------------------------------------------------*/ |
| { |
| float u = 1.0e0, comp; |
| |
| do { |
| u *= 0.5e0; |
| comp = 1.0e0 + u; |
| } |
| while( comp != 1.0e0 ); |
| return( (double)u*2.0e0 ); |
| } |
| /*-------------------- end of function r1mach ------------------------*/ |
| |
| int min0( n, a, b, c, d, e, f, g, h, i, j, k, l, m, o, p ) |
| /* |
| PURPOSE |
| Determine the minimum of the arguments a-p. |
| |
| INPUT |
| n Number of arguments to check 1 <= n <= 15. |
| a-p Integer arguments of which the minumum is desired. |
| |
| RETURNS |
| min0 Minimum of a thru p. |
| */ |
| int n, a, b, c, d, e, f, g, h, i, j, k, l, m, o, p; |
| { |
| int mt; |
| |
| if( n < 1 || n > 15 ) return( -1 ); |
| mt = a; |
| if( n == 1 ) return( mt ); |
| |
| if( mt > b ) mt = b; |
| if( n == 2 ) return( mt ); |
| |
| if( mt > c ) mt = c; |
| if( n == 3 ) return( mt ); |
| |
| if( mt > d ) mt = d; |
| if( n == 4 ) return( mt ); |
| |
| if( mt > e ) mt = e; |
| if( n == 5 ) return( mt ); |
| |
| if( mt > f ) mt = f; |
| if( n == 6 ) return( mt ); |
| |
| if( mt > g ) mt = g; |
| if( n == 7 ) return( mt ); |
| |
| if( mt > h ) mt = h; |
| if( n == 8 ) return( mt ); |
| |
| if( mt > i ) mt = i; |
| if( n == 9 ) return( mt ); |
| |
| if( mt > j ) mt = j; |
| if( n == 10 ) return( mt ); |
| |
| if( mt > k ) mt = k; |
| if( n == 11 ) return( mt ); |
| |
| if( mt > l ) mt = l; |
| if( n == 12 ) return( mt ); |
| |
| if( mt > m ) mt = m; |
| if( n == 13 ) return( mt ); |
| |
| if( mt > o ) mt = o; |
| if( n == 14 ) return( mt ); |
| |
| if( mt > p ) mt = p; |
| return( mt ); |
| } |
| int sscal( n, sa, sx, incx ) |
| int n, incx; |
| float sa, *sx; |
| /* |
| PURPOSE |
| Scales a vector by a constant. |
| |
| INPUT |
| n Number of components to scale. |
| sa Scale value. |
| sx Vector to scale. |
| incx Every incx-th element of sx will be scaled. |
| |
| OUTPUT |
| sx Scaled vector. |
| */ |
| { |
| int i; |
| |
| if( n < 1 ) return( 1 ); |
| |
| /* Code for increment not equal to 1.*/ |
| if( incx != 1 ) { |
| if( incx < 0 ) sx += (-n+1)*incx; |
| for( i=0; i<n; i++, sx+=incx ) |
| *sx *= sa; |
| return 0; |
| } |
| /* Code for unit increment. */ |
| for( i=0; i<n; i++, sx++ ) |
| *sx *= sa; |
| return 0; |
| } |
| void vexopy( n, v, x, y, itype ) |
| int n, itype; |
| float *v, *x, *y; |
| /* |
| Purpose: |
| To operate on the vectors x and y. |
| |
| Input: |
| n Number of elements to scale. |
| x First operand vector. |
| y Second operand vector. |
| itype Type of operation to perform: |
| itype = 1 => '+' |
| itype = 2 => '-' |
| |
| Output: |
| v Result vector of x op y. |
| */ |
| { |
| register int i; |
| |
| if( n<1 ) return; |
| |
| if( itype == 1 ) /* ADDITION. */ |
| for( i=0; i<n; i++, x++, y++, v++ ) |
| *v = *x + *y; |
| else /* SUBTRACTION. */ |
| for( i=0; i<n; i++, x++, y++, v++ ) |
| *v = *x - *y; |
| } |
| void vfill( n, v, val ) |
| int n; |
| float *v, val; |
| /* |
| Purpose |
| To fill the FLOAT vector v with the FLOAT value val. |
| Make sure that if you pass a value for val that you cast |
| it to float, viz., |
| vfill( 100, vector, (float)0.0 ); |
| */ |
| { |
| register int i; |
| |
| if( n<1 ) return; |
| for( i=0; i<n; i++, v++ ) |
| *v = val; |
| } |