/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:19 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_swcom2 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_swcom2.h" #include "drswcom2.h" /* program DRSWCOM2 *>> 1996-05-28 DRSWCOM2 Krogh Removed implicit statement. *>> 1996-05-28 DRSWCOM2 Krogh Added external statement. *>> 1995-10-24 DRSWCOM2 Krogh Changed misleading message. *>> 1994-11-02 DRSWCOM2 Krogh Changes to use M77CON *>> 1992-05-15 DRSWCOM2 CLL Removed "stop '... finished'" *>> 1992-03-24 CLL Add parameter MXCASE. *>> 1992-03-12 CLL *>> 1987-10-28 Original time stamp * Demo driver for the SWCOMP package. This code was adapted from the * test driver to reduce the number of tests and the amount of output. * Here we have NCASES = 1 and NN(1) = 6, whereas the test driver had * NCASES = 7 and NN(1) = 0. * The package SWCOMP computes values and derivatives using the * approach described by Wengert. * C. L. Lawson, JPL, 1971. Revised for Fortran77, Jan 1987. * CLL 8/18/87. Added SWASIN & SWACOS to package. * CLL 8/31/87. Added SWRCHN and changed order of args in SWCHN. * CLL 9/1/87. Added SWSINH, SWCOSH, SWTANH. * 1992-05-15 CLL Removed "stop '... finished'" to simplify * comparison of output from different systems. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *--S replaces "?": DR?WCOM2, ?WCOMP *--& ?COPY, ?MXDIF, ?PASCL, ?VECP, ?WACOS, ?WASIN *--& ?WATAN, ?WATN2, ?WCHN, ?WCOS, ?WCOSH, ?WDIF, ?WDIF1, ?WEXP *--& ?WLOG, ?WPRO, ?WPRO1, ?WPWRI, ?WQUO, ?WQUO1, ?WRCHN, ?WSET *--& ?WSIN, ?WSINH, ?WSQRT, ?WSUM, ?WSUM1, ?WTAN, ?WTANH * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* PARAMETER translations */ #define C7 0.7e0 #define C8 0.8e0 #define C9 0.9e0 #define FOUR 4.0e0 #define HALF 0.50e0 #define MXCASE 7 #define NA 4 #define NCASES 1 #define NDIM (NMAX + 1) #define NMAX 6 #define NW 20 #define NX 10 #define ONE 1.0e0 #define TWO 2.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ int main( ) { long int i, ia, icase, icount, ii, ipwr, j, n, np1; float pi, s[295], save1[NDIM], save2[NDIM], temp[NDIM], w[NW][NDIM], x[NX][NDIM], xexp[NDIM], xlog[NDIM]; static long nn[MXCASE]={6,1,2,3,4,5,6}; static float a[NA]={-2.0e0,-1.0e0,0.75e0,2.5e0}; static float fix3[NA]={0.0e0,0.0e0,0.0e0,0.0e0}; static LOGICAL32 as[NA]={FALSE,TRUE,TRUE,FALSE}; static LOGICAL32 ac[NA]={FALSE,FALSE,TRUE,TRUE}; static LOGICAL32 detail = FALSE; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const A = &a[0] - 1; LOGICAL32 *const Ac = &ac[0] - 1; LOGICAL32 *const As = &as[0] - 1; float *const Fix3 = &fix3[0] - 1; long *const Nn = &nn[0] - 1; float *const S = &s[0] - 1; float *const Save1 = &save1[0] - 1; float *const Save2 = &save2[0] - 1; float *const Temp = &temp[0] - 1; float *const Xexp = &xexp[0] - 1; float *const Xlog = &xlog[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ printf(" DRSWCOM2.. Demo driver for the whole SWCOMP package.\n" " Will print the numerical error in various calculations.\n"); pi = atanf( ONE )*FOUR; Fix3[1] = pi; Fix3[4] = -pi; if (detail) { for (i = 1; i <= 6; i++) { for (j = 1; j <= 40; j++) { S[j] = ZERO; } spascl( i - 1, s ); ii = ((i - 1)*i)/2 + 1; svecp( s, ii, "0DPASCL" ); } } /* Set X(*,1) and X(*,2) */ x[0][0] = C8; x[1][0] = C7; for (i = 2; i <= NDIM; i++) { x[0][i - 1] = C8*x[0][i - 2]; x[1][i - 1] = -C7*x[1][i - 2]; } /* Loop through different values of N */ for (icase = 1; icase <= NCASES; icase++) { n = Nn[icase]; printf(" \n >>>>>> Tests with Nderiv = %3ld\n \n", n); np1 = n + 1; /* Set W(*,1) and W(*,2) */ w[0][0] = C9; w[1][0] = C9; w[0][1] = C9*C9; w[1][1] = -C9; for (i = 3; i <= NDIM; i++) { w[0][i - 1] = -C9*w[0][i - 2]; w[1][i - 1] = ZERO; } /* Test SWSET */ swset( n, C9, -C9, &w[2][0] ); printf(" Error in SET =%11.3g\n", smxdif( np1, &w[1][0], 1, &w[2][0], 1 )); /* Test SWSUM, SWSUM1, SWDIF, SWDIF1, SWPRO1 * */ swsum1( n, FOUR, &x[0][0], &x[2][0] ); swsum( n, &x[2][0], &x[1][0], &x[2][0] ); swdif( n, &x[2][0], &x[0][0], &x[2][0] ); swdif1( n, FOUR, &x[2][0], &x[2][0] ); swpro1( n, -ONE, &x[2][0], &x[2][0] ); printf(" Error in -(4-(((4+x)+y))-x) - y =%11.3g\n", smxdif( np1, &x[2][0], 1, &x[1][0], 1 )); /* Test SWQUO1, SWPRO1, and SWPWRI * */ swquo1( n, TWO, &x[0][0], &x[2][0] ); swpwri( n, -1, &x[2][0], &x[3][0] ); swpro1( n, TWO, &x[3][0], &x[4][0] ); printf(" Error in 2*((2/x)**-1) - x =%11.3g\n", smxdif( np1, &x[4][0], 1, &x[0][0], 1 )); /* Test SWPRO and SWQUO */ swpro( n, &x[0][0], &x[1][0], &x[2][0] ); swquo( n, &x[2][0], &x[1][0], &x[3][0] ); printf(" Error in (X1*X2)/X2 - X1 =%11.3g\n", smxdif( np1, &x[3][0], 1, &x[0][0], 1 )); if (detail) { printf("X1, X2, X3 = X1 * X2, X4 = X3/X2\n"); for (j = 1; j <= 4; j++) { printf(" %3ld", j); for (i = 1; i <= np1; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } } swquo( n, &x[2][0], &x[1][0], &x[2][0] ); printf(" Error in (X1*X2)/X2 - X1 =%11.3g\n", smxdif( np1, &x[2][0], 1, &x[0][0], 1 )); if (detail) { printf(" \n"); printf("X3 = X3/X2\n"); printf(" %3d", 3); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[2][i - 1]); } printf("\n"); } /* Test SWLOG and SWEXP, uses SWCHN * */ swlog( n, &w[0][0], &w[1][0] ); swexp( n, &w[1][0], &w[2][0] ); printf(" Error in Exp(Log(2)) - 2 =%11.3g\n", smxdif( np1, &w[2][0], 1, &w[0][0], 1 )); if (detail) { svecp( &w[0][0], np1, "0W1 = T = 2." ); svecp( &w[1][0], np1, "0W2 = LOG(W1)" ); svecp( &w[2][0], np1, "0W3 = EXP(W2). SHOULD MATCH W1." ); svecp( &w[3][0], np1, "0W4 = 0.+W1. SCALAR ADD." ); } /* Test SWCHN and SWRCHN */ swset( n, HALF, ONE, &x[7][0] ); swexp( n, &x[7][0], xexp ); swset( n, HALF, ONE, &x[2][0] ); swrchn( n, xexp, &x[2][0] ); if (detail) { printf(" \n"); printf("X1, X2 = EXP(X1), X3 = Reversion to Log\n"); for (j = 1; j <= 3; j++) { printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } } swset( n, Xexp[1], ONE, &x[3][0] ); swlog( n, &x[3][0], xlog ); swset( n, x[3][0], ONE, &x[5][0] ); swrchn( n, xlog, &x[5][0] ); if (detail) { printf(" \n"); printf("X4, X5 = LOG(X4), X6 = Reversion to Exp\n"); for (j = 4; j <= 6; j++) { printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } } printf(" Error in Log(x) - Rev. of Exp(x) =%11.3g\n", smxdif( np1, xlog, 1, &x[2][0], 1 )); printf(" Error in Exp(x) - Rev. of Log(x) =%11.3g\n", smxdif( np1, xexp, 1, &x[5][0], 1 )); scopy( np1, &x[0][0], 1, &x[6][0], 1 ); scopy( np1, &x[6][0], 1, save1, 1 ); swchn( n, xexp, &x[6][0] ); scopy( np1, &x[6][0], 1, save2, 1 ); swrchn( n, xexp, &x[6][0] ); printf(" Error in X7;CHN(XEXP,X7);RCHN(XEXP,X7) =%11.3g\n", smxdif( np1, save1, 1, &x[6][0], 1 )); if (detail) { printf(" \n"); printf("X7 =\n"); j = 7; printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save1[i]); } printf("\n"); printf("call SWCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save2[i]); } printf("\n"); printf("call SWRCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } scopy( np1, &x[0][0], 1, &x[6][0], 1 ); scopy( np1, &x[6][0], 1, save1, 1 ); swchn( n, xexp, &x[6][0] ); scopy( np1, &x[6][0], 1, save2, 1 ); swchn( n, xlog, &x[6][0] ); printf(" Error in X7;CHN(XEXP,X7);CHN(XLOG,X7) =%11.3g\n", smxdif( np1, save1, 1, &x[6][0], 1 )); if (detail) { printf(" \n"); printf("X7 =\n"); j = 7; printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save1[i]); } printf("\n"); printf("call SWCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save2[i]); } printf("\n"); printf("call SWCHN ( N, X5, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } scopy( np1, &x[0][0], 1, &x[6][0], 1 ); scopy( np1, &x[6][0], 1, save1, 1 ); swrchn( n, xexp, &x[6][0] ); scopy( np1, &x[6][0], 1, save2, 1 ); swrchn( n, xlog, &x[6][0] ); printf(" Error in X7;RCHN(XEXP,X7);RCHN(XLOG,X7)=%11.3g\n", smxdif( np1, save1, 1, &x[6][0], 1 )); if (detail) { printf(" \n"); printf("X7 =\n"); j = 7; printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save1[i]); } printf("\n"); printf("call SWRCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save2[i]); } printf("\n"); printf("call SWRCHN ( N, X5, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } /* Test SWPRO and SWSQRT * */ scopy( np1, &w[0][0], 1, &w[3][0], 1 ); for (i = 1; i <= 3; i++) { swpro( n, &w[3][0], &w[3][0], &w[3][0] ); if (i == 2) scopy( np1, &w[3][0], 1, temp, 1 ); if (detail) svecp( &w[3][0], np1, "0W4 = W4*W4" ); } swsqrt( n, &w[3][0], &w[4][0] ); printf(" Error in Sqrt(x*x) - x =%11.3g\n", smxdif( np1, temp, 1, &w[4][0], 1 )); if (detail) svecp( &w[4][0], np1, "0W5=SQRT(W4)" ); /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Test of SWPWRI( ,0, , ) * */ swpwri( n, 0, &x[0][0], &x[2][0] ); swset( n, ONE, ZERO, &x[3][0] ); printf(" Error in x**0 - 1 =%11.3g\n", smxdif( np1, &x[2][0], 1, &x[3][0], 1 )); /* Loop through tests of SWPWRI */ for (ipwr = 1; ipwr <= 3; ipwr++) { printf(" \n Test SWPWRI using I = %3ld\n", ipwr); swpwri( n, ipwr, &x[0][0], &x[2][0] ); scopy( np1, &x[0][0], 1, &x[3][0], 1 ); for (icount = 2; icount <= ipwr; icount++) { swpro( n, &x[0][0], &x[3][0], &x[3][0] ); } printf(" Error in x**I - x * x * ... * x =%11.3g\n", smxdif( np1, &x[2][0], 1, &x[3][0], 1 )); swpwri( n, -ipwr, &x[0][0], &x[4][0] ); swpro( n, &x[4][0], &x[3][0], &x[5][0] ); swset( n, ONE, ZERO, &x[6][0] ); printf(" Error in x**(-I) * x**I - 1 =%11.3g\n", smxdif( np1, &x[5][0], 1, &x[6][0], 1 )); } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Loop through quadrants to test trig functions. * */ for (ia = 1; ia <= NA; ia++) { printf(" \n Using angle argument,%11.4f\n", A[ia]); scopy( np1, &w[0][0], 1, &w[4][0], 1 ); w[4][0] = A[ia]; /* Test SWSIN and SWASIN */ swsin( n, &w[4][0], &w[5][0] ); if (As[ia]) { swasin( n, &w[5][0], &w[12][0] ); printf(" Error in Asin(Sin(x)) - x =%11.3g\n", smxdif( np1, &w[4][0], 1, &w[12][0], 1 )); } /* Test SWCOS and SWACOS */ swcos( n, &w[4][0], &w[6][0] ); if (Ac[ia]) { swacos( n, &w[6][0], &w[13][0] ); printf(" Error in Acos(Cos(x)) - x =%11.3g\n", smxdif( np1, &w[4][0], 1, &w[13][0], 1 )); } /* Test SWTAN */ swtan( n, &w[4][0], &w[7][0] ); swquo( n, &w[5][0], &w[6][0], &w[8][0] ); printf(" Error in Tan(x) - Sin(x)/Cos(x) =%11.3g\n", smxdif( np1, &w[7][0], 1, &w[8][0], 1 )); /* Test SWATN2, SWATAN */ swatn2( n, &w[5][0], &w[6][0], &w[8][0] ); swatan( n, &w[7][0], &w[9][0] ); scopy( np1, &w[8][0], 1, temp, 1 ); Temp[1] += Fix3[ia]; printf(" Error in FIX + Atan2(y,x) - Atan(y/x) =%11.3g\n", smxdif( np1, temp, 1, &w[9][0], 1 )); swset( n, ONE, ZERO, &w[10][0] ); swatn2( n, &w[7][0], &w[10][0], &w[11][0] ); printf(" Error in Atan2(y/x, 1) - Atan(y/x) =%11.3g\n", smxdif( np1, &w[11][0], 1, &w[9][0], 1 )); if (detail) { svecp( &w[4][0], np1, "0W5 = ANGLE(IA)" ); svecp( &w[5][0], np1, "0W6 = SIN(W5)" ); svecp( &w[12][0], np1, "0W13 = ASIN(W6). Compare with W5." ); svecp( &w[6][0], np1, "0W7 = COS(W5)" ); svecp( &w[13][0], np1, "0W13 = Acos(W7). Compare with W5." ); svecp( &w[7][0], np1, "0W8 = W6/W7" ); svecp( &w[8][0], np1, "0W9 = ATAN2(W6,W7). SHOULD MATCH W5" ); svecp( &w[9][0], np1, "0W10 = ATAN(W8). SHOULD MATCH W5 OR DIFFER BY 3.14159265" ); svecp( &w[11][0], np1, "0W12=ATAN2(W8,1.)" ); } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Test SWSINH, SWCOSH, TANH * First check identity: cosh**2 = sinh**2 + 1 * */ swsinh( n, &x[0][0], &x[2][0] ); swpro( n, &x[2][0], &x[2][0], &x[5][0] ); swsum1( n, ONE, &x[5][0], &x[7][0] ); swcosh( n, &x[0][0], &x[3][0] ); swpro( n, &x[3][0], &x[3][0], &x[6][0] ); printf(" Error in cosh**2 - sinh**2 - 1 =%11.3g\n", smxdif( np1, &x[6][0], 1, &x[7][0], 1 )); /* Check identity: Tanh = Sinh/Cosh * */ swquo( n, &x[2][0], &x[3][0], &x[8][0] ); swtanh( n, &x[0][0], &x[4][0] ); printf(" Error in tanh - sinh/cosh =%11.3g\n", smxdif( np1, &x[4][0], 1, &x[8][0], 1 )); } exit(0); } /* end of function */ /* ================================================================== */ float /*FUNCTION*/ smxdif( long n, float x[], long incx, float y[], long incy) { long int i, ix, iy; float smxdif_v, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const X = &x[0] - 1; float *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Compute max norm of difference between the N-vectors x and y. * The vectors are stored with a storage increment of INCX and INCY * between successive components. * C. L. Lawson, JPL, Sept 1987. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ temp = ZERO; ix = 1 - incx; iy = 1 - incy; for (i = 1; i <= n; i++) { ix += incx; iy += incy; temp = fmaxf( temp, fabsf( X[ix] - Y[iy] ) ); } smxdif_v = temp; return( smxdif_v ); } /* end of function */