/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:15 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dgeco.h" #include /* PARAMETER translations */ #define ONE 1.0e0 #define ZERO 0.0e0 /* end of PARAMETER translations */ void /*FUNCTION*/ dgeco( double *a, long lda, long n, long ipvt[], double *rcond, double z[]) { #define A(I_,J_) (*(a+(I_)*(lda)+(J_))) long int info, j, k, kb, kp1, l; double anorm, ek, s, sm, t, wk, wkm, ynorm; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Ipvt = &ipvt[0] - 1; double *const Z = &z[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 DGECO Krogh Added external statement. *>> 1994-10-20 DGECO Krogh Changes to use M77CON *>> 1987-08-18 DGECO Lawson Initial code. *--D replaces "?": ?GECO, ?GEFA, ?AXPY,?DOT,?SCAL,?ASUM * * Given an Nth order matrix, A, this subroutine does the following * three steps: * 1. Compute a norm of A for internal use at Step 3. * 2. Call _GEFA to compute an LU factorization of A that overwrites * A in the array A(). * 3. Compute RCOND, an estimate of the reciprocal of the condition * number of A. * The total computing time for Steps 1 and 3 is relatively * economical, being about (9/N) times the time needed for Step 2. * ------------------------------------------------------------------ * Subroutine arguments * * A(,) [inout] An array of size at least N x N. On entry must * contain an N x N matrix, A, to be factored. On return will * contain the LU factors of A. * * LDA [in] Leading dimensioning parameter for the array A(,). * * N [in] The order of the matrix, A. * * IPVT() [out] An integer array of length at least N, into which * will be stored a record of the row interchanges made during * factorization of A. * * RCOND [out] An estimate of the reciprocal of the condition * number of the matrix, A. * Will satisfy ZERO .le. RCOND .le. ONE. * For the system A*X = B , relative perturmations * in A and B of size EPSILON may cause * relative perturbations in X of size EPSILON/RCOND . * If RCOND is so small that the expression * ONE + RCOND * truncated to working precision is not distinguishable from * ONE, then A is singular to working precision. * In particular, RCOND will be returned as ZERO if exact * singularity is detected or if the estimate underflows. * * Z() [out] An array of length at least N. Needed in this * subroutine as work space. On return will contain a vector * Z satisfying NORM(A*Z) = RCOND * NORM(A) * NORM(Z) . * If A is close to a singular matrix, then Z will be * an approximate null vector. * ------------------------------------------------------------------ * 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 directly: DGEFA, DAXPY,DDOT,DSCAL,DASUM * ------------------------------------------------------------------ */ /* ------------------------------------------------------------------ * COMPUTE 1-NORM OF A * */ anorm = ZERO; for (j = 1; j <= n; j++) { anorm = fmax( anorm, dasum( n, &A(j - 1,0), 1 ) ); } /* Replace A by its LU factorization * */ dgefa( a, lda, n, ipvt, &info ); /* RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . * ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . * TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE * CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE * TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID * OVERFLOW. * * SOLVE TRANS(U)*W = E * */ ek = ONE; for (j = 1; j <= n; j++) { Z[j] = ZERO; } for (k = 1; k <= n; k++) { if (Z[k] != ZERO) ek = sign( ek, -Z[k] ); if (fabs( ek - Z[k] ) > fabs( A(k - 1,k - 1) )) { s = fabs( A(k - 1,k - 1) )/fabs( ek - Z[k] ); dscal( n, s, z, 1 ); ek *= s; } wk = ek - Z[k]; wkm = -ek - Z[k]; s = fabs( wk ); sm = fabs( wkm ); if (A(k - 1,k - 1) != ZERO) { wk /= A(k - 1,k - 1); wkm /= A(k - 1,k - 1); } else { wk = ONE; wkm = ONE; } kp1 = k + 1; for (j = kp1; j <= n; j++) { sm += fabs( Z[j] + wkm*A(j - 1,k - 1) ); Z[j] += wk*A(j - 1,k - 1); s += fabs( Z[j] ); } if (s < sm) { t = wkm - wk; wk = wkm; for (j = kp1; j <= n; j++) { Z[j] += t*A(j - 1,k - 1); } } Z[k] = wk; } s = ONE/dasum( n, z, 1 ); dscal( n, s, z, 1 ); /* SOLVE TRANS(L)*Y = W * */ for (kb = 1; kb <= n; kb++) { k = n + 1 - kb; if (k < n) Z[k] += ddot( n - k, &A(k - 1,k), 1, &Z[k + 1], 1 ); if (fabs( Z[k] ) > ONE) { s = ONE/fabs( Z[k] ); dscal( n, s, z, 1 ); } l = Ipvt[k]; t = Z[l]; Z[l] = Z[k]; Z[k] = t; } s = ONE/dasum( n, z, 1 ); dscal( n, s, z, 1 ); ynorm = ONE; /* SOLVE L*V = Y * */ for (k = 1; k <= n; k++) { l = Ipvt[k]; t = Z[l]; Z[l] = Z[k]; Z[k] = t; if (k < n) daxpy( n - k, t, &A(k - 1,k), 1, &Z[k + 1], 1 ); if (fabs( Z[k] ) > ONE) { s = ONE/fabs( Z[k] ); dscal( n, s, z, 1 ); ynorm *= s; } } s = ONE/dasum( n, z, 1 ); dscal( n, s, z, 1 ); ynorm *= s; /* SOLVE U*Z = V * */ for (kb = 1; kb <= n; kb++) { k = n + 1 - kb; if (fabs( Z[k] ) > fabs( A(k - 1,k - 1) )) { s = fabs( A(k - 1,k - 1) )/fabs( Z[k] ); dscal( n, s, z, 1 ); ynorm *= s; } if (A(k - 1,k - 1) != ZERO) { Z[k] /= A(k - 1,k - 1); } else { Z[k] = ONE; } t = -Z[k]; daxpy( k - 1, t, &A(k - 1,0), 1, &Z[1], 1 ); } /* MAKE ZNORM = 1.0 */ s = ONE/dasum( n, z, 1 ); dscal( n, s, z, 1 ); ynorm *= s; if (anorm == ZERO) { *rcond = ZERO; } else { *rcond = ynorm/anorm; } return; #undef A } /* end of function */