/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:44 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "durev.h" #include /* COMMON translations */ struct t_ucom1 { long int n, m1, m2; } ucom1; /* end of COMMON translations */ void /*FUNCTION*/ durev( double *xt, double *tx, long ndim, double *rcond, long iwork[], double *work) { #define XT(I_,J_) (*(xt+(I_)*(ndim)+(J_))) #define TX(I_,J_) (*(tx+(I_)*(ndim)+(J_))) #define WORK(I_,J_,K_) (*(work+(I_)*(ucom1.n)*(ucom1.n)+(J_)*(ucom1.n)+\ (K_))) long int i, it, itxx, ix, ixtt, j, k; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iwork = &iwork[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-03-30 DUREV Krogh Added external statement. *>> 1994-10-20 DUREV Krogh Changes to use M77CON *>> 1994-08-04 DUREV CLL Changed name of common from UCOM to UCOM1. *>> 1992-02-17 CLL *>> 1990-12-12 CLL DUREV Initial code. * This subr does what is sometimes called series reversion. * Regarding N variables x sub j as functions of * N variables t sub i, and given a set of values of the t's and * the values of the x's and the 1st and 2nd partial derivatives of * the x's with respect to the t's evaluated at this set of t values, * this subr computes values at this point of the 1st and 2nd partial * derivatives of the t's with respect to the x's. * It is required that the matrix of 1st partials of the x's with * respect to the t's must be nonsingular. * ------------------------------------------------------------------ * Subroutine Arguments * * XT(,) [in float] Array containing values of N variables x sub 1 * through x sub N, along with their partial derivatives of orders * 1 and 2 with respect to N variables t sub 1 through t sub N, all * evaluated at the set of t values given in TX(). * Data for x sub i is in XT(1:(N+2)*(N+1)/2),i). * * TX(,) [inout float] On entry TX(1,i), i = 1,..., N, must contain * values of the N variables t sub 1 thru t sub N. * Other values in the array TX() on entry are irrelevant. * This subr will compute 1st and 2nd partials of the t's with * respect to the x's and store these results in * TX(2:(N+2)*(N+1)/2),1:N). * * NDIM [in, integer] Leading dimension for the arrays XT() and TX(). * Require NDIM .ge. (N+2)*(N+1)/2. * * RCOND [out float] Estimate of the reciprocal condition number of * the matrix of 1st partials of x's with respect to the t's. * RCOND will satisfy 0.0 .le. RCOND .le. 1.0. * Values near 1.0 indicate a well conditioned matrix, small values * indicate poor conditioning, zero indicates a singular matrix. In * this latter case the subr will return without computing partials * of the t's with respect to the x's. * * IWORK() [scratch, integer] Integer work space for this subroutine. * * WORK(N,N,3) [scratch, float] Floating-point work space for this * subroutine. We use this as 3 NxN arrays. In the calls to DGECO * and DGEI the last argument needs only N scratch locations. * ------------------------------------------------------------------ * N, M1, M2 in common /UCOM1/ * * N [in] Number of components in the (conceptual) x and t vectors. * M1, M2 [in] We assume 0 .le. M1 .le. M2 .le. 2. In other subrs of * the [D/S]UCOMP package, M1 and M2 select computation of partial * derivatives of orders M1 through M2, assuming all needed partial * derivatives of orders less than M1 are available. * This subr differs from others in the [D/S]UCOMP package in that * derivs of order 0, i.e., the values of t sub i for i = 1, ..., N, * associated with x sub j for j = 1, ..., N, must always be input to * this subr. Thus this subr treats M1 = 0 like M1 = 1. * If M2 .eq. 0 this subr does nothing. * If M2 .eq. 1 this subr computes first partials of t w.r.t. x. * If M1 .le. 1 and M2 .eq. 2 this subr computes first and second * partials of t w.r.t. x. * If M1 .eq. 2 and M2 .eq. 2 this subr assumes first partials of t * w.r.t. x are available and computes second partials of t w.r.t. x. * ------------------------------------------------------------------ *--D replaces "?": ?UREV, ?DOT, ?GECO, ?GEI * Also uses ERMSG * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ */ if (ucom1.m2 == 0) return; if (ucom1.m1 <= 1) { /* Copy the matrix of 1st partials of x w.r.t. t to WORK2(*,*,1) * */ for (ix = 1; ix <= ucom1.n; ix++) { for (it = 1; it <= ucom1.n; it++) { WORK(0,it - 1,ix - 1) = XT(ix - 1,it); } } /* Compute an LU factorization and RCOND for this Jacobian matrix. * */ dgeco( &WORK(0,0,0), ucom1.n, ucom1.n, iwork, rcond, &WORK(2,0,0) ); if (*rcond == 0.0e0) { ermsg( "DUREV", 6, 0, "Singular Jacobian matrix.", '.' ); return; } /* Compute the inverse Jacobian matrix in WORK(*,*,1) and copy * its elements into TX(). These elements are the 1st partials of * the t's w.r.t. the x's. * */ dgei( &WORK(0,0,0), ucom1.n, ucom1.n, iwork, &WORK(2,0,0) ); for (it = 1; it <= ucom1.n; it++) { for (ix = 1; ix <= ucom1.n; ix++) { TX(it - 1,ix) = WORK(0,ix - 1,it - 1); } } } else { /* Copy 1st partials of t's w.r.t. x's from TX() to WORK(*,*,1). * */ for (it = 1; it <= ucom1.n; it++) { for (ix = 1; ix <= ucom1.n; ix++) { WORK(0,ix - 1,it - 1) = TX(it - 1,ix); } } } /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ if (ucom1.m2 < 2) return; /* -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Main loop on K to compute 2nd partials * of t sub K w.r.t. the x's. * */ for (k = 1; k <= ucom1.n; k++) { /* Compute weighted sum of the Hessian matrices of the x's * using the 1st partial of t sub K w.r.t. x sub i as the * weight for the Hessian matrix of x sub i. * Store the resulting symmetric matrix in WORK(*,*,2). * */ ixtt = 1 + ucom1.n; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { ixtt += 1; WORK(1,j - 1,i - 1) = ddot( ucom1.n, &WORK(0,0,k - 1), ucom1.n, &XT(0,ixtt - 1), ndim ); WORK(1,i - 1,j - 1) = WORK(1,j - 1,i - 1); } } /* Multiply the matrix in WORK(*,*,2) by WORK(*,*,1) from the * right and by the transpose of WORK(*,*,1) from the left and * invert the sign. Put the result of the first multiplication * into WORK(*,*,3) and the lower triangle of the (symmetric) * final result into TX(*,K). * */ for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= ucom1.n; j++) { WORK(2,j - 1,i - 1) = ddot( ucom1.n, &WORK(1,0,i - 1), ucom1.n, &WORK(0,j - 1,0), 1 ); } } itxx = 1 + ucom1.n; for (i = 1; i <= ucom1.n; i++) { for (j = 1; j <= i; j++) { itxx += 1; TX(k - 1,itxx - 1) = -ddot( ucom1.n, &WORK(0,i - 1,0), 1, &WORK(2,j - 1,0), 1 ); } } } return; #undef WORK #undef XT #undef TX } /* end of function */