/*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_sucom2 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_sucom2.h" #include "drsucom2.h" /* program DRSUCOM2 *>> 1996-06-20 DRSUCOM2 Krogh Format changes for C conversion. *>> 1996-05-28 DRSUCOM2 Krogh Added external statement. *>> 1994-11-02 DRSUCOM2 Krogh Changes to use M77CON *>> 1994-08-12 DRSUCOM2 CLL New subroutine: SUSETN *>> 1992-05-15 CLL Removed "stop '... finished'" *>> 1992-04-21 CLL *>> 1992-03-16 CLL *>> 1992-03-12 CLL *>> 1987-10-30 Original time stamp * Demo driver for the SUCOMP package. This code was adapted from the * test driver to reduce the number of tests and the amount of output. * Here we have NCASES = 2 and MM(3:4) = (0,2), whereas the test driver * had NCASES = 12 and MM(3:4) = (0,0). * The SUCOMP package computes first and second * partial derivatives. * C.L.Lawson,JPL, 1969 Dec 4 * CLL, JPL, Jan 1987, Modified for Fortran 77 * CLL, JPL, Sept 1987, Added SUSINH, SUCOSH, SUTANH, SUTAN * 1992-05-15 CLL Removed "stop '... finished'" to simplify * comparison of output from different systems. * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - *--S replaces "?": DR?UCOM2, ?UCOMP, ?MXDIF, ?USETN *--& ?COPY , ?UACOS, ?UASIN, ?UATAN, ?UATN2, ?UCOS *--& ?UCOSH, ?UDIF , ?UDIF1, ?UEXP , ?ULOG , ?UPRO , ?UPRO1 *--& ?UPWRI, ?UQUO , ?UQUO1, ?USET , ?USIN , ?USINH, ?USQRT *--& ?USUM , ?USUM1, ?UTAN , ?UTANH * -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* PARAMETER translations */ #define A1 (-2.62e0) #define A2 1.57e0 #define C8 0.8e0 #define C9 0.9e0 #define IMAX 10 #define JMAX 72 #define NCASES 12 #define ONE 1.0e0 #define TEN 10.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ int main( ) { long int i, iang, iang4, icount, ipwr, ipwr7, j, key, kk, l, l2, loc, m1, m2, n; float ang, c, errmax, error, test[IMAX], u[JMAX][IMAX]; static long mm[24]={-1,-1,0,2,1,1,2,2,-1,-1,0,0,1,2,-1,-1,0,1, 2,2,-1,-1,0,2}; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Mm = &mm[0] - 1; float *const Test = &test[0] - 1; /* end of OFFSET VECTORS */ /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ printf(" DRSUCOM2.. Demo driver for the SUCOMP package.\n Computation of partial derivatives.\n Will print the numerical error in various calculations.\n These errors should be zero or small.\n"); c = TEN; n = 3; l = 1 + n + (n*(n + 1))/2; m1 = 0; m2 = 2; for (kk = 1; kk <= (2*NCASES - 1); kk += 2) { if (Mm[kk] < 0) { printf(" \n Test U-computation. Set U(,) = 0.\n"); /* Zero the U() array. */ for (i = 1; i <= IMAX; i++) { for (j = 1; j <= JMAX; j++) { u[j - 1][i - 1] = ZERO; } } /* Store values into U(,1) and U(,2) */ u[0][0] = C8; u[1][0] = -C9; for (i = 2; i <= l; i++) { u[0][i - 1] = -C8*u[0][i - 2]; u[1][i - 1] = -C9*u[1][i - 2]; } goto L_800; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Set M1 and M2 for new case. */ m1 = Mm[kk]; m2 = Mm[kk + 1]; susetn( n, m1, m2 ); printf(" \n N, M1, M2 =%10ld%10ld%10ld\n \n", n, m1, m2); if (m2 == 0) { l2 = 1; } else if (m2 == 1) { l2 = 1 + n; } else if (m2 == 2) { l2 = 1 + n + (n*(n + 1))/2; } errmax = ZERO; supro( &u[0][0], &u[1][0], &u[2][0] ); suquo( &u[2][0], &u[0][0], &u[3][0] ); error = smxdif( l2, &u[3][0], 1, &u[1][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in (u1*u2)/u1 - u2 =%11.3g\n", error); susqrt( &u[0][0], &u[4][0] ); supwri( 2, &u[4][0], &u[5][0] ); error = smxdif( l2, &u[5][0], 1, &u[0][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in sqrtf(u1)**2 - u1 =%11.3g\n", error); susin( &u[0][0], &u[6][0] ); sucos( &u[0][0], &u[7][0] ); suatn2( &u[6][0], &u[7][0], &u[8][0] ); error = smxdif( l2, &u[8][0], 1, &u[0][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in atan2f(sinf(u1),cosf(u1)) - u1 =%11.3g\n", error); suexp( &u[0][0], &u[9][0] ); sulog( &u[9][0], &u[10][0] ); error = smxdif( l2, &u[10][0], 1, &u[0][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in logf(expf(u1) - u1 =%11.3g\n", error); suquo( &u[6][0], &u[7][0], &u[11][0] ); suatan( &u[11][0], &u[12][0] ); error = smxdif( l2, &u[12][0], 1, &u[0][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in atanf(sinf(u1)/cosf(u1)) - u1 =%11.3g\n", error); sutan( &u[0][0], &u[43][0] ); error = smxdif( l2, &u[11][0], 1, &u[43][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in sinf(u1)/cosf(u1) - tanf(u1) =%11.3g\n", error); susinh( &u[0][0], &u[44][0] ); sucosh( &u[0][0], &u[45][0] ); sutanh( &u[0][0], &u[46][0] ); suquo( &u[44][0], &u[45][0], &u[47][0] ); error = smxdif( l2, &u[47][0], 1, &u[46][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in sinhf(u1)/coshf(u1) - tanhf(u1) =%11.3g\n", error); susum( &u[0][0], &u[1][0], &u[13][0] ); sudif( &u[13][0], &u[0][0], &u[14][0] ); error = smxdif( l2, &u[14][0], 1, &u[1][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in ((u1+u2)-u1) - u2 =%11.3g\n", error); susum1( c, &u[0][0], &u[15][0] ); sudif1( c, &u[15][0], &u[16][0] ); supro1( c, &u[16][0], &u[17][0] ); suquo1( c, &u[17][0], &u[18][0] ); suquo1( -ONE, &u[18][0], &u[19][0] ); error = smxdif( l2, &u[19][0], 1, &u[0][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in -1/(c/(c*(c-(c+u1)))) - u1 =%11.3g\n", error); suasin( &u[6][0], &u[20][0] ); error = smxdif( l2, &u[20][0], 1, &u[0][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in asinf(sinf(u1)) - u1 =%11.3g\n", error); suacos( &u[7][0], &u[21][0] ); error = smxdif( l2, &u[21][0], 1, &u[0][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in acosf(cosf(u1)) - u1 =%11.3g\n", error); /* Test of SUSET. We first use SUSUM1 to copy U(,1) * to U(,23+KEY) so it will have nonzero contents. * After calling SUSET it should be all zero except for 1 or 2 * elements. C will be stored to U(1,23+KEY). * When KEY = 0, everything else will be set to zero. * When KEY = 1, 2, or 3, and it will also set * U(KEY+1,23+KEY) = 1.0 * */ for (key = 0; key <= 3; key++) { susum1( ZERO, &u[0][0], &u[key + 22][0] ); suset( c, key, &u[key + 22][0] ); for (i = 1; i <= l; i++) { Test[i] = ZERO; } Test[1] = c; if (key > 0) Test[key + 1] = ONE; error = smxdif( l2, &u[key + 22][0], 1, test, 1 ); errmax = fmaxf( errmax, error ); printf(" Error in SUSET(C,KEY,...) =%11.3g\n", error); } /* Check quadrant resolution of SUATN2. * * Store values into U(,1). */ ang = A1 - A2; for (iang = 1; iang <= 4; iang++) { iang4 = 4*(iang - 1); ang += A2; scopy( l, &u[0][0], 1, &u[iang4 + 26][0], 1 ); u[iang4 + 26][0] = ang; susin( &u[iang4 + 26][0], &u[iang4 + 27][0] ); sucos( &u[iang4 + 26][0], &u[iang4 + 28][0] ); suatn2( &u[iang4 + 27][0], &u[iang4 + 28][0], &u[iang4 + 29][0] ); error = smxdif( l2, &u[iang4 + 29][0], 1, &u[iang4 + 26][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in atan2f(sinf(u1),cosf(u1)) - u1 =%11.3g\n", error); } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Test of SUPWRI( ,0, , ) * */ supwri( 0, &u[0][0], &u[48][0] ); suset( ONE, 0, &u[49][0] ); error = smxdif( l2, &u[48][0], 1, &u[49][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in x**0 - 1 =%11.3g\n", error); /* Loop through tests of SUPWRI */ for (ipwr = 1; ipwr <= 3; ipwr++) { ipwr7 = 7*(ipwr - 1); printf(" \n Test SUPWRI using I = %3ld\n", ipwr); supwri( ipwr, &u[0][0], &u[ipwr7 + 50][0] ); scopy( l2, &u[0][0], 1, &u[ipwr7 + 55][0], 1 ); loc = 0; for (icount = 2; icount <= ipwr; icount++) { supro( &u[0][0], &u[loc + ipwr7 + 55][0], &u[loc + ipwr7 + 56][0] ); loc += 1; } error = smxdif( l2, &u[ipwr7 + 50][0], 1, &u[loc + ipwr7 + 55][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in x**I - x * x * ... * x =%11.3g\n", error); supwri( -ipwr, &u[0][0], &u[ipwr7 + 52][0] ); supro( &u[ipwr7 + 52][0], &u[ipwr7 + 50][0], &u[ipwr7 + 53][0] ); suset( ONE, 0, &u[ipwr7 + 54][0] ); error = smxdif( l2, &u[ipwr7 + 53][0], 1, &u[ipwr7 + 54][0], 1 ); errmax = fmaxf( errmax, error ); printf(" Error in x**(-I) * x**I - 1 =%11.3g\n", error); } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ printf(" ***>>> ERRMAX =%11.3g\n", errmax); L_800: ; } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 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 */