/*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 <math.h>
#include "fcrt.h"
#include <stdio.h>
#include <stdlib.h>
#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 */