/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:30:04 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dged.h" #include #include /* PARAMETER translations */ #define ONE 1.0e0 #define TEN 10.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dged( double *a, long lda, long n, long ipvt[], double det[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int _l0, _l1, i, k; double aterm, dt1, dt2; static double fbig, fsmall, step, tbig, tsmall; static LOGICAL32 first = TRUE; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Det = &det[0] - 1; long *const Ipvt = &ipvt[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. *>> 2001-05-25 DGED Krogh Minor change for making .f90 version. *>> 1996-03-30 DGED Krogh Added external statement. *>> 1994-10-20 DGED Krogh Changes to use M77CON *>> 1992-04-14 DGED CLL Using (*) in dimensions rather than (1). *>> 1990-01-19 CLL added test to avoid looping when some A(I,I) = 0. *>> 1987-08-20 DGED Lawson Initial code. *--D replaces "?": ?GED * * DGED computes the determinant of the matrix A * using the LU factorization of A given in the array A(). * ------------------------------------------------------------------ * Subroutine arguments * * A(,) [in] An array of size at least N x N. On entry must * contain the LU factors of an N x N matrix, A. It is * expected that this factorization will have been computed by * use of _GEFA, either directly or indirectly via use of * _GEFCO, _GEFS or _GEFSC. This subr will not alter the * contents of A(,) * * LDA [in] Leading dimensioning parameter for the array A(,). * * N [in] The order of the original matrix, A. * * IPVT() [in] An integer array of length at least N, containg a * record of the row interchanges made during factorization of * A. * * DET() [out] Array of length 2. On return will contain a * representation of the determinant of A of the form * DETERMINANT = DET(1) * 10.0**DET(2) * DET(2) will be integer valued: plus, minus or zero. * If the determinant is zero then DET(1) = 0.0 and DET(2) = 0. * If the determinant is nonzero the scaling chosen will not be * determined exclusively by the value of the determinant, but * will tend to choose DET(2) = 0.0 when the determinant is not * extremely large or small in magnitude. This is a change * from the normalization used in LINPACK. * ------------------------------------------------------------------ * Machine dependent values * * The following values are set depending on the machine underflow * and overflow limits. * FBIG = 10.0**STEP * FSMALL = 10.0**(-STEP) = 1.0/FBIG * TBIG = FBIG/10.0 = 10.0**(STEP-1) * TSMALL = FSMALL*10.0 = 10.0**(-(STEP-1)) = 1.0/TBIG * * where STEP is a positive integer value chosen so that FBIG**2 is * slightly smaller than the overflow limit and FBIG**(-2) is * slightly larger than the underflow limit. * If one wanted a machine-independent version of this code, one * could set STEP small enough so that FBIG**2 and FBIG**(-2) will * not be off-scale for any "reasonable" computer. For instance one * might set STEP = 10. However, setting STEP as large as possible * reduces the number of scaling multiplications and increases the * likelihood that the subroutine will be able to set DET(2) = zero. * These parameters added by C. L. Lawson, JPL, Aug 1987. * ------------------------------------------------------------------ * LINPACK. THIS VERSION DATED 08/14/78 . * CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. * Ref: LINPACK Users' Guide, by J. J. Dongarra, C. B. Moler, * J. R. Bunch, and G. W. Stewart, publ by Soc. for Indust. and Appl. * Math, Philadelphia, 1979. * Adapted from LINPACK for the JPL Math77 library by * C. L. Lawson, JPL, Aug 1987. * ------------------------------------------------------------------ * Subprograms referenced: D1MACH * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * Set machine-dependent values. */ if (first) { first = FALSE; k = fmin( log10( DBL_MAX/TEN ), -log10( DBL_MIN* TEN ) )/2; step = k - 1; fbig = pow(TEN,step); fsmall = ONE/fbig; tbig = fbig/TEN; tsmall = fsmall*TEN; } /* Initialize DT1 and DT2. */ dt1 = ONE; dt2 = ZERO; /* Begin main loop. */ for (i = 1; i <= n; i++) { if (Ipvt[i] != i) dt1 = -dt1; aterm = A(i - 1,i - 1); if (aterm == ZERO) { dt1 = ZERO; dt2 = ZERO; goto L_160; } /* Scale ATERM, if necessary, to bring its magnitude into * the range from TSMALL to TBIG. * */ if (fabs( aterm ) >= tbig) { L_10: aterm *= fsmall; dt2 += step; if (fabs( aterm ) >= tbig) goto L_10; } else if (fabs( aterm ) <= tsmall) { L_20: aterm *= fbig; dt2 -= step; if (fabs( aterm ) <= tsmall) goto L_20; } /* Scale DT1, if necessary, to bring its magnitude into * the range from TSMALL to TBIG. * */ if (fabs( dt1 ) >= tbig) { dt1 *= fsmall; dt2 += step; } else if (fabs( dt1 ) < tsmall) { dt1 *= fbig; dt2 -= step; } /* Update the determinant value. */ dt1 *= aterm; if (dt1 == ZERO) { dt2 = ZERO; goto L_160; } } L_160: ; Det[1] = dt1; Det[2] = dt2; return; #undef A } /* end of function */