/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:09 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_ddasl6 s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include #include #include "p_ddasl6.h" #include "drddasl6.h" /* PARAMETER translations */ #define LIW (30 + NEQ) #define LRW (45 + (MAXORD + MAXCON + 4)*NEQ) #define MAXCON 0 #define MAXORD 5 #define N 40 #define NDIG 6 #define NEQ (SQ(N)) #define TOL (powi(10.e0,-NDIG)) /* end of PARAMETER translations */ /* COMMON translations */ struct t_counts { long int ke, kf, ks, kp; } counts; /* end of COMMON translations */ int main( ) { long int i, idid, imax, info[16], iwork[LIW], j, jmax; double atol[NEQ], maxu, rtol[NEQ], rwork[LRW], t, tout[11], u[N][N], uprime[N][N]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Atol = &atol[0] - 1; long *const Info = &info[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rtol = &rtol[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Tout = &tout[0] - 1; /* end of OFFSET VECTORS */ /*>> 2009-10-21 DRDDASL6 Krogh at do 120 upper limit km => km+1. *>> 2008-10-26 DRDDASL6 Krogh Moved Format statements up for C conv. *>> 2008-10-24 DRDDASL6 Krogh Removed an INCLUDE statement & cDEC$... *>> 2008-08-27 DRDDASL6 assign initial y' value, change user codes *>> 2001-10-11 DR DDASL6 R. J. Hanson Example 6 Code, with Download */ /* Driver program for the non-linear elliptic or Bratu problem. * The unknowns are the array entries of the solution stacked up * as one large vector. In the evaluation routine, the two-dimensional * aspect of the problem is used when evaluating residuals and * also linear system solves. The linear solves are done using * the GMRES algorithm. *--D replaces "?": ?DASLX, ?DASSF6, DR?DASL6 */ /*++S Default NDIG = 2 *++ Default NDIG = 6 *++ Substitute for NDIG below */ /* Note that the N here and in routine DDASSF6 must have the same value. */ /* The size of the work space LRW does not have the NEQ**2=N**4 term * because with INFO(5)=5, the Jacobian matrix data is stored in the * evaluation program ?DASF (= ?DASF6). */ /* Absolute and relative tolerances: */ for (i = 1; i <= NEQ; i++) { Atol[i] = TOL; Rtol[i] = TOL; } /* Setup default options: */ for (i = 1; i <= 16; i++) { Info[i] = 0; } /* Change some options from their defaults: * Partials and solutions are done in DDASSF6 below. */ Info[5] = 5; /* Restrict the BDF formulas to a lower order than 5. * This is not needed in this example. */ Info[9] = 4; /* Assign the starting values of U(*), UPRIME(*) in the * initial call to DDASSF6. */ /* Set output points for return of solution. */ Tout[1] = 0.e0; for (i = 2; i <= 11; i++) { Tout[i] = Tout[i - 1] + 1.e0; } for (i = 2; i <= 11; i++) { t = Tout[i - 1]; ddaslx( ddassf6, NEQ, &t, (double*)u, (double*)uprime, Tout[i], info, rtol, atol, &idid, rwork, LRW, iwork, LIW ); } /* Locate maximum value of solution. * It should be near the centroid of the unit square. */ maxu = u[0][0]; imax = 1; jmax = 1; for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { if (u[j - 1][i - 1] > maxu) { maxu = u[j - 1][i - 1]; imax = i; jmax = j; } } } /* Summarize results and amount of work requested. */ printf(" Example Results for a non-linear Elliptic PDE Problem \n\n Maximum value of solution with N by N system of ODEs\n"); printf(" u(Imax,Jmax), Imax, Jmax N, (N by N grid) T final\n"); printf(" %6.2f %4ld %4ld %4ld %6.2f\n\n", maxu, imax, jmax, N, Tout[11]); printf(" Number of requested evaluations\n Residuals Fresh Solve Steps Solutions GMRES products\n"); printf("%12ld %12ld %12ld %12ld\n", counts.ke, counts.kf, counts.ks, counts.kp); exit(0); } /* end of function */ /* PARAMETER translations */ #define KM 25 /* end of PARAMETER translations */ void /*FUNCTION*/ ddassf6( double *t, double u1[], double upr1[], double del1[], double *d, long *ldd, double *cj, long *ires, double rwork[], long iwork[]) { #define u(I_, J_) (*(u1+(I_)*(N)+(J_))) #define upr(I_, J_) (*(upr1+(I_)*(N)+(J_))) #define del(I_, J_) (*(del1+(I_)*(N)+(J_))) #define D(I_,J_) (*(d+(I_)*(*ldd)+(J_))) long int i, j, k, l; double cc, cs[KM], hh[KM + 1][KM-(0)+1], nrm, qq[KM][N][N], rnrm, rss, sn[KM], ss, tol, tt, ww[N][N], y[KM]; static double hsq3, hsqlam, hx, lambda, reltol, rr[N + 1-(0)+1][N + 1-(0)+1], scl[N][N], u0[N][N], xd[N][N]; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Cs = &cs[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Sn = &sn[0] - 1; double *const Y = &y[0] - 1; /* end of OFFSET VECTORS */ /* Routine for the non-linear elliptic Bratu problem. * This version computes initial data, * residuals, and solves the linear systems with GMRES. * A diagonal pre-conditioner is used, applied from both * sides implicitly. Other functions are left to the DDASLX. */ /* The problem size is N**2 equations. Up to KM iterations are used * in the GMRES algorithm. The results after KM steps are * returned to the integrator. (Have KM .le. N**2; usually KM * will be much smaller.) */ /* If N is changed in the calling program, have it the same value * here. */ /* This is the setup call. * It is used to compute quantities that need be set * once during the integration. */ if (*ires == 0) { counts.ke = 0; counts.kf = 0; counts.ks = 0; counts.kp = 0; hx = 1.e0/(1.e0 + N); hsq3 = 3.e0*SQ(hx); lambda = 6.e0; hsqlam = hsq3*lambda; /* Set initial values for U(*,*) and UPR(*,*): */ for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { u(j - 1, i - 1) = 0.e0; upr(j - 1, i - 1) = lambda; } } /* Set zero boundary values in array RR(,). * This makes it possible to use a natural double loop, and * avoid special code for the edges. */ for (j = 0; j <= (N + 1); j++) { rr[j][0] = 0.e0; rr[j][N + 1] = 0.e0; rr[N + 1][j] = 0.e0; rr[0][j] = 0.e0; } /* Relative tolerance for residual norm reduction in GMRES solver. * (For a fixed N, "optimize" RELT, ATOL(), RTOL(), and KM, * inside DDASSF6.) */ reltol = 9e-2; return; } /* The system residual value: */ if (*ires == 1) { /* Count number of residual evaluations: */ counts.ke += 1; /* Compute required exponential functions. Store the U(,) values * in a bordered array XD(,) so that the special edge condtions * are accounted for. The zero edge conditions in RR(,) are used. */ /* (Automatic differentiation can start here.) */ for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { del(j - 1, i - 1) = hsqlam*exp( u(j - 1, i - 1) ) - hsq3*upr(j - 1, i - 1) - 8.e0*u(j - 1, i - 1); rr[j][i] = u(j - 1, i - 1); } } /* This is the internal nine point formula: */ for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { del(j - 1, i - 1) += rr[j - 1][i - 1] + rr[j - 1][i] + rr[j - 1][i + 1] + rr[j][i - 1] + rr[j][i + 1] + rr[j + 1][i - 1] + rr[j + 1][i] + rr[j + 1][i + 1]; } } return; /* End of IF block for IRES .eq. 1. */ } /* The "Jacobian computation" is noted with this entry. * All that is needed in this case is to save the scalar * multiplying the partials with respect to y'. */ if (*ires == 2) { /* This passes the error code for singularity (IRES=0) to the integrator. * It is an issue as the user code solves the linear systems. * In this example we do not know about singularity until the solve step. * We precompute quantities that are expensive. */ *ires = 0; counts.kf += 1; for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { u0[j - 1][i - 1] = hsqlam*exp( u(j - 1, i - 1) ) - hsq3**cj - 8.e0; /* The pre-conditioner used is obtained from the matrix diagonal. * This yields a scaled problem, A = D * B * D. The system * is solved with respect to B, applying scaling as needed during * the GMRES algorithm. The diagonal matrix D is stored in SCL(,). */ scl[j - 1][i - 1] = 1.e0/sqrt( -u0[j - 1][i - 1] ); /* This array is used to hold the solution in GMRES. */ xd[j - 1][i - 1] = 0.e0; } } return; /* End of IF block for IRES .eq. 2. */ } /* Solve the corrector equation. */ if (*ires == 4) { /* Note number of corrector solve steps. */ counts.ks += 1; /* Use a rearranged version of the Saad and Schultz GMRES algorithm. */ rss = 0.e0; for (j = 1; j <= N; j++) { /* Compute the first residual, matrix * first approximation - RHS. * The solution is started at zero. * The difference in sign for this residual is accounted for in the * final approximate solution formula. */ for (i = 1; i <= N; i++) { /* Apply scaling to the first residual. */ rr[j][i] = -del(j - 1, i - 1)*scl[j - 1][i - 1]; /* Compute the sum of squares, for the eventual length * of the first residual. */ rss += SQ(rr[j][i]); } } nrm = sqrt( rss ); /* Set counter for number of GMRES basis vectors used. */ k = 0; for (i = 1; i <= (KM + 1); i++) { hh[i - 1][0] = 0.e0; } hh[0][0] = nrm; /* Set the relative tolerance for residual that terminates GMRES. * This is passed from the main program unit. */ tol = reltol*nrm; /* This is the primary iteration loop for the GMRES algorithm. * Total number of matrix-vector products required by GMRES: */ L_130: counts.kp += 1; /* See if residual is now small enough. */ if (fabs( hh[k][0] ) <= tol) goto L_280; /* Reciprocate factor for efficiency during scaling. */ rnrm = 1.e0/nrm; for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { tt = rnrm*rr[j][i]; qq[k][j - 1][i - 1] = tt; /* Account for pre-scaling before applying the product. */ rr[j][i] = tt*scl[j - 1][i - 1]; } } k += 1; /* Compute the product, matrix * column k of basis matrix. */ for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { /* Apply matrix product and post-scaling together. */ ww[j - 1][i - 1] = scl[j - 1][i - 1]*(u0[j - 1][i - 1]* rr[j][i] + rr[j - 1][i - 1] + rr[j][i - 1] + rr[j + 1][i - 1] + rr[j + 1][i] + rr[j - 1][i] + rr[j - 1][i + 1] + rr[j][i + 1] + rr[j + 1][i + 1]); } } /* Keep a working copy of product for orthogonalization updates. */ for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { rr[j][i] = ww[j - 1][i - 1]; } } /* Compute updates and the new column of Hessenberg matrix. */ for (l = 1; l <= k; l++) { tt = 0.e0; for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { tt += ww[j - 1][i - 1]*qq[l - 1][j - 1][i - 1]; } } hh[l - 1][k] = tt; for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { rr[j][i] += -tt*qq[l - 1][j - 1][i - 1]; } } } /* Apply previous rotations, that triangularized the Hessenberg matrix, * to the new column. */ for (i = 1; i <= (k - 1); i++) { cc = Cs[i]; ss = Sn[i]; tt = hh[i - 1][k]*cc + hh[i][k]*ss; hh[i][k] = -hh[i - 1][k]*ss + hh[i][k]*cc; hh[i - 1][k] = tt; } /* Compute the length of the updated residual vector. */ rss = 0.e0; for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { rss += SQ(rr[j][i]); } } /* Update the new entry in the developing Hessenberg matrix. */ nrm = sqrt( rss ); hh[k][k] = nrm; /* Compute the plane rotation that returns the Hessenberg * matrix to upper triangular form. */ tt = sqrt( rss + SQ(hh[k - 1][k]) ); cc = hh[k - 1][k]/tt; ss = hh[k][k]/tt; hh[k - 1][k] = tt; /* Apply rotation to the updated RHS. */ tt = hh[k - 1][0]*cc + hh[k][0]*ss; /* This is the new residual vector length. It is checked for * smallness at the next loop iteration. */ hh[k][0] = -hh[k - 1][0]*ss + hh[k][0]*cc; hh[k - 1][0] = tt; /* Save rotation parameters for later application to new * Hessenberg matrix columns. */ Cs[k] = cc; Sn[k] = ss; /* This check stops the iterations if there is no longer enough * array space. * This ends the primary iteration loop for the GMRES algorithm. */ if (k < KM) goto L_130; L_280: ; /* Solve the K by K upper triangular system for the multipliers * of the orthogonal basis of the approximate solution. */ for (j = k; j >= 1; j--) { tt = 0.e0; for (i = j + 1; i <= k; i++) { tt += Y[i]*hh[j - 1][i]; } Y[j] = (hh[j - 1][0] - tt)/hh[j - 1][j]; } /* Compute approximate solution. This is: initial approximation - * product of basis vectors multiplied by the y values. * This accounts for the change of sign on the first residual. */ for (l = 1; l <= k; l++) { for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { xd[j - 1][i - 1] += -qq[l - 1][j - 1][i - 1]*Y[l]; } } } /* This solution is therefore a starting approximation for the next one. */ for (j = 1; j <= N; j++) { for (i = 1; i <= N; i++) { /* Account for pre-scaling in solution returned. */ del(j - 1, i - 1) = xd[j - 1][i - 1]*scl[j - 1][i - 1]; xd[j - 1][i - 1] = 0.e0; } } /* End of IF block for IRES .eq. 4; end of GMRES algorithm application. */ } return; #undef D } /* end of function */