# include # include # include # include # include # include "roots_rc.h" /******************************************************************************/ void roots_rc ( int n, double x[], double fx[], double *ferr, double xnew[], double q[] ) /******************************************************************************/ /* Purpose: ROOTS_RC solves a system of nonlinear equations using reverse communication. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 January 2013 Author: Original FORTRAN77 version by Gaston Gonnet. C version by John Burkardt. Reference: Gaston Gonnet, On the Structure of Zero Finders, BIT Numerical Mathematics, Volume 17, Number 2, June 1977, pages 170-183. Parameters: Input, int N, the number of equations. Input, double X[N]. Before the first call, the user should set X to an initial guess or estimate for the root. Thereafter, the input value of X should be the output value of XNEW from the previous call. Input, double FX[N], the value of the function at XNEW. Output, double *FERR, the function error, that is, the sum of the absolute values of the most recently computed function vector. Output, double XNEW[N], a new point at which a function value is requested. Workspace, double Q[(2*N+2)*(N+2)]. Before the first call for a given problem, the user must set Q(2*N+1,1) to 0.0. */ { double damp; int i; int j; int jsma; int jsus; int lda; double sump; double t; lda = 2 * n + 2; *ferr = 0.0; for ( i = 0; i < n; i++ ) { *ferr = *ferr + fabs ( fx[i] ); } /* Initialization if Q(2*N+1,1) = 0.0. */ if ( q[2*n+1+0*lda] == 0.0 ) { for ( i = 1; i <= n; i++ ) { for ( j = 1; j <= n + 1; j++ ) { q[i-1+(j-1)*lda] = 0.0; q[i+(j-1)*lda] = 0.0; } q[i-1+(i-1)*lda] = 100.0; q[i+n-1+(i-1)*lda] = 1.0; } for ( j = 1; j <= n; j++ ) { q[2*n+(j-1)*lda] = r8_huge ( ); } for ( j = 1; j <= n; j++ ) { q[2*n+1+(j-1)*lda] = ( double ) ( n ); } for ( i = 1; i <= n; i++ ) { q[i+n-1+n*lda] = x[i-1]; } for ( i = 1; i <= n; i++ ) { q[i-1+n*lda] = fx[i-1]; } q[2*n+n*lda] = *ferr; q[2*n+1+n*lda] = 0.0; damp = 0.99; } else { jsus = 1; for ( i = 2; i <= n + 1; i++ ) { if ( ( double ) ( 2 * n ) <= q[2*n+1+(i-1)*lda] ) { q[2*n+(i-1)*lda] = r8_huge ( ); } if ( q[2*n+1+(jsus-1)*lda] < ( n + 3 ) / 2 ) { jsus = i; } if ( ( n + 3 ) / 2 <= q[2*n+1+(i-1)*lda] && q[2*n+(jsus-1)*lda] < q[2*n+(i-1)*lda] ) { jsus = i; } } for ( i = 1; i <= n; i++ ) { q[i+n-1+(jsus-1)*lda] = x[i-1]; q[i-1+(jsus-1)*lda] = fx[i-1]; } q[2*n+(jsus-1)*lda] = *ferr; q[2*n+1+(jsus-1)*lda] = 0; jsma = 1; damp = 0.0; for ( j = 1; j <= n + 1; j++ ) { if ( r8_huge ( ) / 10.0 < q[2*n+(j-1)*lda] ) { damp = 0.99; } if ( q[2*n+(j-1)*lda] < q[2*n+(jsma-1)*lda] ) { jsma = j; } } if ( jsma != n + 1 ) { for ( i = 1; i <= 2 * n + 2; i++ ) { t = q[i-1+(jsma-1)*lda]; q[i-1+(jsma-1)*lda] = q[i-1+n*lda]; q[i-1+n*lda] = t; } } } for ( i = 1; i <= n; i++ ) { q[i-1+(n+1)*lda] = q[i-1+n*lda]; } /* Call the linear equation solver, which should not destroy the matrix in Q(1:N,1:N), and should overwrite the solution into Q(1:N,N+2). */ r8mat_fs ( lda, n, q, q+(n+1)*lda ); sump = 0.0; for ( i = 1; i <= n; i++ ) { sump = sump + q[i-1+(n+1)*lda]; } if ( fabs ( 1.0 - sump ) <= 1.0E-10 ) { printf ( "\n" ); printf ( "ROOTS_RC - Fatal error!\n" ); printf ( " SUMP almost exactly 1.\n" ); printf ( " SUMP = %g\n", sump ); exit ( 1 ); } for ( i = 1; i <= n; i++ ) { xnew[i-1] = q[i+n-1+n*lda]; for ( j = 1; j <= n; j++ ) { xnew[i-1] = xnew[i-1] - q[i+n-1+(j-1)*lda] * q[j-1+(n+1)*lda]; } /* If system not complete, damp the solution. */ xnew[i-1] = xnew[i-1] / ( 1.0 - sump ) * ( 1.0 - damp ) + q[i+n-1+n*lda] * damp; } for ( j = 1; j <= n + 1; j++ ) { q[2*n+1+(j-1)*lda] = q[2*n+1+(j-1)*lda] + 1.0; } return; } /******************************************************************************/ double r8_huge ( ) /******************************************************************************/ /* Purpose: R8_HUGE returns a "huge" R8. Discussion: The value returned by this function is NOT required to be the maximum representable R8. This value varies from machine to machine, from compiler to compiler, and may cause problems when being printed. We simply want a "very large" but non-infinite number. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 October 2007 Author: John Burkardt Parameters: Output, double R8_HUGE, a "huge" R8 value. */ { double value; value = 1.0E+30; return value; } /******************************************************************************/ void r8mat_fs ( int lda, int n, double a[], double x[] ) /******************************************************************************/ /* Purpose: R8MAT_FS factors and solves a system with one right hand side. Discussion: This routine uses partial pivoting, but no pivot vector is required. This routine differs from R8MAT_FSS in two ways: * only one right hand side is allowed; * the input matrix A is not modified. An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 January 2013 Author: John Burkardt Parameters: Input, int LDA, the leading dimension of the matrix. Input, int N, the order of the matrix. N must be positive. Input, double A[N*N], the coefficient matrix of the linear system. The matrix is stored in an LDAxN array. Input/output, double X[N], on input, the right hand side of the linear system. On output, the solution of the linear system. */ { double *a2; int i; int ipiv; int j; int jcol; double piv; double t; a2 = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { a2[i+j*n] = a[i+j*lda]; } } for ( jcol = 1; jcol <= n; jcol++ ) { /* Find the maximum element in column I. */ piv = fabs ( a2[jcol-1+(jcol-1)*n] ); ipiv = jcol; for ( i = jcol+1; i <= n; i++ ) { if ( piv < fabs ( a2[i-1+(jcol-1)*n] ) ) { piv = fabs ( a2[i-1+(jcol-1)*n] ); ipiv = i; } } if ( piv == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_FS - Fatal error!\n" ); fprintf ( stderr, " Zero pivot on step %d\n", jcol ); exit ( 1 ); } /* Switch rows JCOL and IPIV, and X. */ if ( jcol != ipiv ) { for ( j = 1; j <= n; j++ ) { t = a2[jcol-1+(j-1)*n]; a2[jcol-1+(j-1)*n] = a2[ipiv-1+(j-1)*n]; a2[ipiv-1+(j-1)*n] = t; } t = x[jcol-1]; x[jcol-1] = x[ipiv-1]; x[ipiv-1] = t; } /* Scale the pivot row. */ t = a2[jcol-1+(jcol-1)*n]; a2[jcol-1+(jcol-1)*n] = 1.0; for ( j = jcol+1; j <= n; j++ ) { a2[jcol-1+(j-1)*n] = a2[jcol-1+(j-1)*n] / t; } x[jcol-1] = x[jcol-1] / t; /* Use the pivot row to eliminate lower entries in that column. */ for ( i = jcol+1; i <= n; i++ ) { if ( a2[i-1+(jcol-1)*n] != 0.0 ) { t = - a2[i-1+(jcol-1)*n]; a2[i-1+(jcol-1)*n] = 0.0; for ( j = jcol+1; j <= n; j++ ) { a2[i-1+(j-1)*n] = a2[i-1+(j-1)*n] + t * a2[jcol-1+(j-1)*n]; } x[i-1] = x[i-1] + t * x[jcol-1]; } } } /* Back solve. */ for ( jcol = n; 2 <= jcol; jcol-- ) { for ( i = 1; i < jcol; i++ ) { x[i-1] = x[i-1] - a2[i-1+(jcol-1)*n] * x[jcol-1]; } } free ( a2 ); return; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the GNU LGPL license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }