/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:13 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dwcom2 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_dwcom2.h" #include "drdwcom2.h" /* program DRDWCOM2 *>> 1996-05-28 DRDWCOM2 Krogh Removed implicit statement. *>> 1996-05-28 DRDWCOM2 Krogh Added external statement. *>> 1995-10-24 DRDWCOM2 Krogh Changed misleading message. *>> 1994-11-02 DRDWCOM2 Krogh Changes to use M77CON *>> 1992-05-15 DRDWCOM2 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 DWCOMP 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 DWCOMP 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 DWASIN & DWACOS to package. * CLL 8/31/87. Added DWRCHN and changed order of args in DWCHN. * CLL 9/1/87. Added DWSINH, DWCOSH, DWTANH. * 1992-05-15 CLL Removed "stop '... finished'" to simplify * comparison of output from different systems. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *--D 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; double 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 double a[NA]={-2.0e0,-1.0e0,0.75e0,2.5e0}; static double 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 */ double *const A = &a[0] - 1; LOGICAL32 *const Ac = &ac[0] - 1; LOGICAL32 *const As = &as[0] - 1; double *const Fix3 = &fix3[0] - 1; long *const Nn = &nn[0] - 1; double *const S = &s[0] - 1; double *const Save1 = &save1[0] - 1; double *const Save2 = &save2[0] - 1; double *const Temp = &temp[0] - 1; double *const Xexp = &xexp[0] - 1; double *const Xlog = &xlog[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ printf(" DRDWCOM2.. Demo driver for the whole DWCOMP package.\n" " Will print the numerical error in various calculations.\n"); pi = atan( 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; } dpascl( i - 1, s ); ii = ((i - 1)*i)/2 + 1; dvecp( 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 DWSET */ dwset( n, C9, -C9, &w[2][0] ); printf(" Error in SET =%11.3g\n", dmxdif( np1, &w[1][0], 1, &w[2][0], 1 )); /* Test DWSUM, DWSUM1, DWDIF, DWDIF1, DWPRO1 * */ dwsum1( n, FOUR, &x[0][0], &x[2][0] ); dwsum( n, &x[2][0], &x[1][0], &x[2][0] ); dwdif( n, &x[2][0], &x[0][0], &x[2][0] ); dwdif1( n, FOUR, &x[2][0], &x[2][0] ); dwpro1( n, -ONE, &x[2][0], &x[2][0] ); printf(" Error in -(4-(((4+x)+y))-x) - y =%11.3g\n", dmxdif( np1, &x[2][0], 1, &x[1][0], 1 )); /* Test DWQUO1, DWPRO1, and DWPWRI * */ dwquo1( n, TWO, &x[0][0], &x[2][0] ); dwpwri( n, -1, &x[2][0], &x[3][0] ); dwpro1( n, TWO, &x[3][0], &x[4][0] ); printf(" Error in 2*((2/x)**-1) - x =%11.3g\n", dmxdif( np1, &x[4][0], 1, &x[0][0], 1 )); /* Test DWPRO and DWQUO */ dwpro( n, &x[0][0], &x[1][0], &x[2][0] ); dwquo( n, &x[2][0], &x[1][0], &x[3][0] ); printf(" Error in (X1*X2)/X2 - X1 =%11.3g\n", dmxdif( 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"); } } dwquo( n, &x[2][0], &x[1][0], &x[2][0] ); printf(" Error in (X1*X2)/X2 - X1 =%11.3g\n", dmxdif( 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 DWLOG and DWEXP, uses DWCHN * */ dwlog( n, &w[0][0], &w[1][0] ); dwexp( n, &w[1][0], &w[2][0] ); printf(" Error in Exp(Log(2)) - 2 =%11.3g\n", dmxdif( np1, &w[2][0], 1, &w[0][0], 1 )); if (detail) { dvecp( &w[0][0], np1, "0W1 = T = 2." ); dvecp( &w[1][0], np1, "0W2 = LOG(W1)" ); dvecp( &w[2][0], np1, "0W3 = EXP(W2). SHOULD MATCH W1." ); dvecp( &w[3][0], np1, "0W4 = 0.+W1. SCALAR ADD." ); } /* Test DWCHN and DWRCHN */ dwset( n, HALF, ONE, &x[7][0] ); dwexp( n, &x[7][0], xexp ); dwset( n, HALF, ONE, &x[2][0] ); dwrchn( 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"); } } dwset( n, Xexp[1], ONE, &x[3][0] ); dwlog( n, &x[3][0], xlog ); dwset( n, x[3][0], ONE, &x[5][0] ); dwrchn( 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", dmxdif( np1, xlog, 1, &x[2][0], 1 )); printf(" Error in Exp(x) - Rev. of Log(x) =%11.3g\n", dmxdif( np1, xexp, 1, &x[5][0], 1 )); dcopy( np1, &x[0][0], 1, &x[6][0], 1 ); dcopy( np1, &x[6][0], 1, save1, 1 ); dwchn( n, xexp, &x[6][0] ); dcopy( np1, &x[6][0], 1, save2, 1 ); dwrchn( n, xexp, &x[6][0] ); printf(" Error in X7;CHN(XEXP,X7);RCHN(XEXP,X7) =%11.3g\n", dmxdif( 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 DWCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save2[i]); } printf("\n"); printf("call DWRCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } dcopy( np1, &x[0][0], 1, &x[6][0], 1 ); dcopy( np1, &x[6][0], 1, save1, 1 ); dwchn( n, xexp, &x[6][0] ); dcopy( np1, &x[6][0], 1, save2, 1 ); dwchn( n, xlog, &x[6][0] ); printf(" Error in X7;CHN(XEXP,X7);CHN(XLOG,X7) =%11.3g\n", dmxdif( 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 DWCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save2[i]); } printf("\n"); printf("call DWCHN ( N, X5, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", x[j - 1][i - 1]); } printf("\n"); } dcopy( np1, &x[0][0], 1, &x[6][0], 1 ); dcopy( np1, &x[6][0], 1, save1, 1 ); dwrchn( n, xexp, &x[6][0] ); dcopy( np1, &x[6][0], 1, save2, 1 ); dwrchn( n, xlog, &x[6][0] ); printf(" Error in X7;RCHN(XEXP,X7);RCHN(XLOG,X7)=%11.3g\n", dmxdif( 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 DWRCHN ( N, X2, X7). X7 =\n"); printf(" %3ld", j); for (i = 1; i <= NDIM; i++) { printf("%15.7g", Save2[i]); } printf("\n"); printf("call DWRCHN ( 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 DWPRO and DWSQRT * */ dcopy( np1, &w[0][0], 1, &w[3][0], 1 ); for (i = 1; i <= 3; i++) { dwpro( n, &w[3][0], &w[3][0], &w[3][0] ); if (i == 2) dcopy( np1, &w[3][0], 1, temp, 1 ); if (detail) dvecp( &w[3][0], np1, "0W4 = W4*W4" ); } dwsqrt( n, &w[3][0], &w[4][0] ); printf(" Error in Sqrt(x*x) - x =%11.3g\n", dmxdif( np1, temp, 1, &w[4][0], 1 )); if (detail) dvecp( &w[4][0], np1, "0W5=SQRT(W4)" ); /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Test of DWPWRI( ,0, , ) * */ dwpwri( n, 0, &x[0][0], &x[2][0] ); dwset( n, ONE, ZERO, &x[3][0] ); printf(" Error in x**0 - 1 =%11.3g\n", dmxdif( np1, &x[2][0], 1, &x[3][0], 1 )); /* Loop through tests of DWPWRI */ for (ipwr = 1; ipwr <= 3; ipwr++) { printf(" \n Test DWPWRI using I = %3ld\n", ipwr); dwpwri( n, ipwr, &x[0][0], &x[2][0] ); dcopy( np1, &x[0][0], 1, &x[3][0], 1 ); for (icount = 2; icount <= ipwr; icount++) { dwpro( n, &x[0][0], &x[3][0], &x[3][0] ); } printf(" Error in x**I - x * x * ... * x =%11.3g\n", dmxdif( np1, &x[2][0], 1, &x[3][0], 1 )); dwpwri( n, -ipwr, &x[0][0], &x[4][0] ); dwpro( n, &x[4][0], &x[3][0], &x[5][0] ); dwset( n, ONE, ZERO, &x[6][0] ); printf(" Error in x**(-I) * x**I - 1 =%11.3g\n", dmxdif( 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]); dcopy( np1, &w[0][0], 1, &w[4][0], 1 ); w[4][0] = A[ia]; /* Test DWSIN and DWASIN */ dwsin( n, &w[4][0], &w[5][0] ); if (As[ia]) { dwasin( n, &w[5][0], &w[12][0] ); printf(" Error in Asin(Sin(x)) - x =%11.3g\n", dmxdif( np1, &w[4][0], 1, &w[12][0], 1 )); } /* Test DWCOS and DWACOS */ dwcos( n, &w[4][0], &w[6][0] ); if (Ac[ia]) { dwacos( n, &w[6][0], &w[13][0] ); printf(" Error in Acos(Cos(x)) - x =%11.3g\n", dmxdif( np1, &w[4][0], 1, &w[13][0], 1 )); } /* Test DWTAN */ dwtan( n, &w[4][0], &w[7][0] ); dwquo( n, &w[5][0], &w[6][0], &w[8][0] ); printf(" Error in Tan(x) - Sin(x)/Cos(x) =%11.3g\n", dmxdif( np1, &w[7][0], 1, &w[8][0], 1 )); /* Test DWATN2, DWATAN */ dwatn2( n, &w[5][0], &w[6][0], &w[8][0] ); dwatan( n, &w[7][0], &w[9][0] ); dcopy( np1, &w[8][0], 1, temp, 1 ); Temp[1] += Fix3[ia]; printf(" Error in FIX + Atan2(y,x) - Atan(y/x) =%11.3g\n", dmxdif( np1, temp, 1, &w[9][0], 1 )); dwset( n, ONE, ZERO, &w[10][0] ); dwatn2( n, &w[7][0], &w[10][0], &w[11][0] ); printf(" Error in Atan2(y/x, 1) - Atan(y/x) =%11.3g\n", dmxdif( np1, &w[11][0], 1, &w[9][0], 1 )); if (detail) { dvecp( &w[4][0], np1, "0W5 = ANGLE(IA)" ); dvecp( &w[5][0], np1, "0W6 = SIN(W5)" ); dvecp( &w[12][0], np1, "0W13 = ASIN(W6). Compare with W5." ); dvecp( &w[6][0], np1, "0W7 = COS(W5)" ); dvecp( &w[13][0], np1, "0W13 = Acos(W7). Compare with W5." ); dvecp( &w[7][0], np1, "0W8 = W6/W7" ); dvecp( &w[8][0], np1, "0W9 = ATAN2(W6,W7). SHOULD MATCH W5" ); dvecp( &w[9][0], np1, "0W10 = ATAN(W8). SHOULD MATCH W5 OR DIFFER BY 3.14159265" ); dvecp( &w[11][0], np1, "0W12=ATAN2(W8,1.)" ); } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Test DWSINH, DWCOSH, TANH * First check identity: cosh**2 = sinh**2 + 1 * */ dwsinh( n, &x[0][0], &x[2][0] ); dwpro( n, &x[2][0], &x[2][0], &x[5][0] ); dwsum1( n, ONE, &x[5][0], &x[7][0] ); dwcosh( n, &x[0][0], &x[3][0] ); dwpro( n, &x[3][0], &x[3][0], &x[6][0] ); printf(" Error in cosh**2 - sinh**2 - 1 =%11.3g\n", dmxdif( np1, &x[6][0], 1, &x[7][0], 1 )); /* Check identity: Tanh = Sinh/Cosh * */ dwquo( n, &x[2][0], &x[3][0], &x[8][0] ); dwtanh( n, &x[0][0], &x[4][0] ); printf(" Error in tanh - sinh/cosh =%11.3g\n", dmxdif( np1, &x[4][0], 1, &x[8][0], 1 )); } exit(0); } /* end of function */ /* ================================================================== */ double /*FUNCTION*/ dmxdif( long n, double x[], long incx, double y[], long incy) { long int i, ix, iy; double dmxdif_v, temp; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const X = &x[0] - 1; double *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 = fmax( temp, fabs( X[ix] - Y[iy] ) ); } dmxdif_v = temp; return( dmxdif_v ); } /* end of function */