/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:40 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "ddasls.h" #include #include void /*FUNCTION*/ ddasls( void (*ddasf)(), long neq, double *t, double y[], double yprime[], long info[], double ftol, double rnktol, double *c, long *ldc, long ltd, long *idid, double rwork[], long lrw, long iwork[], long liw) { #define C(I_,J_) (*(c+(I_)*(*ldc)+(J_))) LOGICAL32 banded, dense; long int _l0, i, i1, i2, iflag, ioff, irank, ires, irevc, iroff, irow, is, j, jcol, k, m, ml, mlp, mp, mu, mup, nrows; double bigval, cj, eps, one, regul, temp, zero; /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Info = &info[0] - 1; long *const Iwork = &iwork[0] - 1; double *const Rwork = &rwork[0] - 1; double *const Y = &y[0] - 1; double *const Yprime = &yprime[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 2008, Math a la Carte, Inc. * R. J. Hanson, Visual Numerics, Inc., 23 June 2008. * F. T. Krogh, Math a la Carte, Inc. */ /*>>2009/11/02 ddasls Krogh Avoided ref. to undefined IWORK locations. *>>2008/10/26 ddasls Krogh Changed 'No' no 'N' in call to dgemv. *>>2008/09/18 ddasls Krogh obtained from Hanson set up single prec. *--D replaces "?":?daslx, ?dasls, ?dasf, ?swap, ?scal, ?axpy, ?copy, *--& ?asum, ?gbfa, ?gbsl, ?gemv, ?gefa, ?gesl */ /* Provide starting values for y' == yprime(*) for a differential- * algebraic equation f(t,y,y') = 0. * This routine is normally to be used before a call to SDASLX * or DDASLX to get consistent initial conditions for y', * i.e. f(t,y,y') = 0 at the initial t. */ /* The assumptions are: * A) The system has index 0 or index 1 * B) Initial values for (t,y) are known but maybe not all * values of y'. * C) All values of y' have some meaningful initial estimates. * The partials f_t, f_y and f_y' are computable and continuous * at (t,y,y'). The rank of f_y' must be positive. * D) The initial values (t,y) are fixed. Any values for y' are * allowed to change to achieve f(t,y,y') = 0. */ /* This routine is organized so that a user can add additional cases in * the evaluation routine DDASF written for DDASLX that separately * computes f_t, f_y and f_y'. The evaluation of f(t,y,y') uses the * same flag as DDASLX. */ /* DDASF:EXTERNAL -- Depending on how INFO(5) is set, for forward or * reverse communication, this routine is used to compute four * terms f, f_t, f_y, and f_y' as requested. It has the form */ /* SUBROUTINE DDASF(T,Y,YPRIME,DELTA,C,LDC,CJ,IRES,RWORK,IWORK) */ /* The actions taken depend on the value of the IRES argument * to DDASF. The values of IRES that occur when using ddasls are: * = 0 Initialize subroutine parameters or any storage requirements. * If this does not apply, then ignore this value and return * without any further action. Initial conditions for t, y and y' * can be assigned to T, Y(:) and YPRIME(:). * = 1 Evaluate f(t,y,y') and return the result in DELTA(1:NEQ). * = 6 Evaluate the NEQ-vector f_t (evaluated at (t,y,y')) and return * the result in DELTA(1:NEQ) * = 7 Evaluate the NEQ by NEQ matrix f_y' and return the result in * C(1:NEQ,1:NEQ), for dense Jacobians. For banded Jacobians * return the result in C(1:ML+MU+1,1:NEQ) according to the * Linpack data structure for banded matrices. * = 8 Evaluate the NEQ by NEQ matrix f_y and return the result in * C(1:NEQ,1:NEQ). For banded Jacobians return the result in * C(1:ML+MU+1,1:NEQ).This matrix is required when the * system has index 1, i.e. when f_y' is rank-deficient. * * NEQ:IN This is the number of equations to be solved. NEQ .ge. 1. * * T:IN This is the current value of the independent variable. * Set this to the value of the initial point. * * Y(.ge. NEQ):IN This array contains the solution components at T. * Initialize to the values at the initial point. * * YPRIME(.ge. NEQ):INOUT This array contains the current derivatives * of the solution components at T. You must provide some intial * values, but they need not satisfy f(t,y,y')=0. * * INFO(.ge. 16):IN The only task of the code is to solve the system * f(t,y,y')=0 at the initial point. Values of INFO(:) are used * for the same purposes as will occur when integrating the system * using DDASLX. * When the Jacobians are banded matrices, |INFO(5)|=4. In that * case a regularization paramter is required. For this use * INFO(14) .eq. 0 to have the parameter with value macheps**(2./3.) * To use other values set INFO(14) > 0 and place the regularization * parameter in RWORK(INFO(14)). * FTOL:IN This is the tolerance required for the L_1 norm of * f(t,y,y'). Try macheps ** (2./3.) * size of values occuring in f. * RNKTOL:IN This is the relative column norm tolerance used for * determining rank deficiency of f_y'. It is also used for * determining the rank of certain rows of the combined matrix * [f_y' * f_y]. If this latter matrix is rank deficient the system * is defined to have an index .gt. 1. As such DDASLX does * not claim to solve it. Try macheps ** (2./3.). * * C(LDC, LTD):INOUT This is the working array where the partial * derivative matrices f_y' and f_y will be returned after * evaluation by ddasf or reverse communication. * * LDC:IN This is the leading dimension of the array C(:,:). It has * a value LDC .ge. NEQ when f_y' has full rank. This occurs with * an index 0 problem. When f_y' is rank deficient then * LDC .ge. 2*NEQ is required. When f_y' is rank deficient, * and the Jacobian matrices are banded, then * LDC .ge. 4*ML+2*MU+4 is required. This occurs when * the problem has index 1 or higher. * */ /* LDT:IN This is the second dimension of the array C(:,:). It has * a value LTD .ge. NEQ when f_y' has full rank. This occurs with * an index 0 problem for either dense or banded problesm. When * f_y' is rank deficient, and the Jacobian matrices are banded, * then LDT .ge. 2*NEQ is required. This occurs when the problem * has index 1 or higher. * * IDID:OUT This is an output flag that shows the status of solving * f(t,y,y')=0 for values of y'. * = 0 The system is index 0 or 1 and consistent values for y' * were computed * =-1 An error or exceptional condition was noted. Details are * flagged by various values in IWORK(1). * .ge.1 Compute requested value using reverse communication. * Re-enter the routine ddasls() for further computations. * IWORK(1) flag, giving more detail: * Currently IOFF=8, IROFF=2. * = 0 The system has index 0 and consistent values for y' * were computed * = 1 The system has index 1 and consistent values for y' * were computed * = 2 The system has index 0 but the system f(t,y,y') is * not consistent using the tolerance FTOL. It should * be consistent so FTOL may be too small for this problem. * Iterations are done while the L_1 norm of f is .lt. * 1/4 the L_1 norm of f from the previous iteration. * = 3 The system has index 1 but the equation f(t,y,y')=0 * is not consistent using the tolerance FTOL. * Iterations are done while the L_1 norm of f is .lt. * 1/4 the L_1 norm of f from the previous iteration. * = 4 The system appears to have an index .gt. 1. * = 5 The system has rank f_y' = 0. This is not a DAE. * = 6 The value NEQ .le. 0, i.e. no system * = 7 LDC is too small. Must be .ge. NEQ for index 0 problems, * dense Jacobians. * Must be .ge. (2*ML+MU+1), banded, index 0. * Here ML, MU are the lower and upper band widths. Provide the * lower ML and upper MU bandwidths by setting IWORK(1)=ML and * IWORK(2)=MU. * Must be .ge. 4*ML+2*MU+4 if problem is banded with index 1. * = 8 LDC is too small. Must be .ge. 2*NEQ for index 1 problems, * dense Jacobians. * Must be .ge. 2*(2*ML+1)+2*(MU+1), banded. * = 9 LRW is too small. Must be .ge. 2*NEQ+IROFF for dense * systems and .ge. 7*NEQ+IROFF for * banded systems. * = 10 LIW is too small. Must be .ge. 2*NEQ+IOFF * = 11, Must have |INFO(5)|=2 or 4. Routine supports dense * or banded matrices only. The user computes derivatives in * forward or reverse communication. * = 12, Must have FTOL .gt. 0. Now .le. 0. * = 13, Must have RNKTOL .ge. 0. Now .lt. 0. * = 14, Must have sub- and super- bandwidth parameters * ML, MU .ge. 0. Now one of them is .lt. 0. * = 15, Must have second dimension, LTD, of C(*,*) .ge. NEQ * for dense or banded systems. Value is now .lt. NEQ. * For banded systems of index 1 C(*,*) .ge. 2*NEQ. This * is checked when the index is .gt. 0, or f_y' is singular. * * RWORK(LRW .ge. 2*NEQ+IROFF):INOUT Working floating point * storage for dense systems, * LRW .ge. 7*NEQ+IROFF for banded systems. * IWORK(LIW .ge. 2*NEQ+IOFF):INOUT Working integer storage * User code can change values in these arrays past these minimum * limits for evaluation or other storage needs. */ /* INTERNAL VARIABLE DICTIONARY - * Floating Point Variables: * DELTA(1:NEQ) of DDASF is associated with RWORK(1+IROFF:NEQ+IROFF). * Used for evaluation of f(t,y,y') and f_t. * EPS - Machine precision for this accuracy. Obtained by a * call to D1MACH. * REGUL - Value of the regularization parameter used if the system * is banded and of index .gt. 1. Default is EPS**(2./3.) * but the user can provide an alternate value in * location RWORK(INFO(14)). * BIGVAL - Temporary value used to contain the largest magnitude * during full pivoting * CJ - A dummy floating point value, not used in this routine, but * it is an argument of DDASF. Set to 0 for all DDASF calls. * ONE, ZERO - Constants, 1 and 0 * TEMP - Temporary value used in swapping row values during * elimination * BIGVAL - Temporary value used for choosing pivots * CJ - Value used as temporary for forward communication with DDASF. * This is not used but the value is assigned 0. * ONE, ZERO - Local copies of 1.d0 and 0.d0. */ /* Integer Variables: * IRES - A flag for directing user action in routine ddasf * I - A dummy index for looping * IS - A dummy subscript for building interleaved matrix array * IROW - The row index of the largest magnitude in pivoting * IFLAG - Indicate a flat zero on the diagonal using * LINPACK LU factorization code, DGEFA. * IOFF - The offset value in IWORK(*) that is * used for storage of integer values and * communication flags. So IWORK(IOFF+1:IOFF+2*NEQ) * is used for pivoting information. * IWORK(1:IOFF) is used for storage and must remain */ /* unchanged during forward or reverse communication. * The present value is IOFF=8. * IROFF - The offset value in RWORK(*) that is used for * storage of floating point parameters. * So RWORK(IROFF+1:IROFF+2*NEQ) is used for storage * of intermediate values and requested values from * user code. The present value is IROFF=2. * * J - A dummy index for looping * JCOL - The column index of the largest magnitude in pivoting * I1,I2 - Temporary indices for placement of interleaved * matrix values * IREVC - Type of matrix, dense or banded, forward * or reverse communication mode * IRANK - Rank of f_y' for index 1 problems * ML, MU - For banded Jacobians, these are temps that hold * the number of lower and upper co-diagonals that * are non-zero. * M - A temporary integer holding subscripts for the * banded interleaved Jacobian. * MLP, MUP - For banded Jacobians, these are temps that hold * the number of lower and upper co-diagonals when * the interleaved and regularized matrix is constructed. * Have MLP=2*ML+1, MUP=max(2*MU,1). * MP - A temporary integer holding subscripts for the * banded interleaved Jacobian. * NROWS - A temporary that holds the number of rows to zero * in the C(:,:) array for computing f_y' and f_y, later. * Logical Variables: * DENSE, BANDED - Logical variables signifying that the * Jacobians are either dense or banded matrices. * Exactly one of these variables are .TRUE. * This information is taken from the input value * INFO(5). */ /* Variable typing, subroutine arguments */ /* Variable typing, local arguments or routines called */ ires = 0; one = 1.e0; zero = 0.e0; /* This defines the starting index in IWORK(*) where * row and column pivoting interchanges are stored. * Locations of IWORK(1:IOFF) are needed for possible reverse * communication storage and signalling. */ ioff = 8; /* This defines the starting index in RWORK(*) where * requested results for f and f_t will be placed. Locations * RWORK(1:IROFF) are needed for storage of routine values. */ iroff = 2; /* This routine only recognizes |IREVC| = 2 or 4. * These values correspond to dense or banded matrices. */ irevc = Info[5]; /* Check that these are the allowed values. */ i = labs( irevc ); dense = i == 2; banded = i == 4; if (!(dense || banded)) { Iwork[1] = 11; *idid = -1; goto L_370; } if (Info[1] == 1) { /* When reverse communication is used IRANK is saved. * With forward communication IRANK is not defined until * the LU factorization of f_y' is completed. */ irank = Iwork[6]; /* Extract band widths. If the Jacobian is dense these are not * used but it does not matter. */ ml = Iwork[7]; mu = Iwork[8]; } /* Give this parameter a defined value. It is not used * but the evaluator routine needs a slot for it. */ cj = zero; /* Re-enter to just past this subroutine return point when reverse * communication is used. */ switch (Info[1]) { case 1: goto L_40; case 2: goto L_70; case 3: goto L_170; case 4: goto L_190; } /* The first entry drops through here, INFO(1) == 0. * Check sizes of arrays and other parameters. Branch to * a block of code that does this and then branches back here. */ goto L_360; /* Set up counters and flags for reverse communication. * This is done just once because INFO(1) is .gt. 0 for * reverse communication entries. */ /* This is where the results for an evaluation of f or f_t * should be placed, RWORK(IWORK(4)+1:IWORK(4)+NEQ) */ L_10: Iwork[4] = iroff; Iwork[5] = 0; /* Make initialization call. This is done even when reverse * communication is used. */ ires = 0; (*ddasf)( t, y, yprime, &Rwork[iroff + 1], c, ldc, &cj, &ires, rwork, iwork ); /* Start of overall loop to compute Newton update for y' */ /* Make a call to ddasf to obtain the value of f(t,y,y'). The * result is returned in RWORK(IROFF+1:IROFF+NEQ). */ L_20: ires = 1; /* If IRES=1 or 6, clear out the spot where the results will go. * These are respectively the function and f_t. */ for (i = 1; i <= neq; i++) { Rwork[i + iroff] = zero; } if (irevc > 0) { (*ddasf)( t, y, yprime, &Rwork[iroff + 1], c, ldc, &cj, &ires, rwork, iwork ); } else { Info[1] = 1; goto L_380; } /* Branch back place if reverse communication in use. */ L_40: ; /* After one Newton update step, check the size of the * function f(t,y,y'). It should be zero but won't be * due to rounding and other approximations. Check its * value against the user tolerance, FTOL. */ if (Iwork[5] > 0) { /* Need to distinguish index 0 and 1 cases here. */ temp = dasum( neq, &Rwork[iroff + 1], 1 ); if (temp <= ftol) { *idid = 0; Iwork[1] = 0; /* Return L_1 norm of f here. */ Rwork[2] = temp; } else { /* After a Newton update step the defect or residual in the DAE * should be small. If it is not decreasing at an acceptable * rate (1/4) then flag a possible inconsistency. */ if (temp < 0.25*Rwork[2]) { /* Save the current norm of f for checking the next iteration. */ Rwork[2] = temp; goto L_50; } *idid = -1; Iwork[1] = 2; if (irank < neq) Iwork[1] = 3; /* Return L_1 norm of f here. */ Rwork[2] = temp; } if (*idid < 0) goto L_370; /* No errors here. Function is small, everything normal. * No error processing needed here. */ Info[1] = 0; return; } /* Branch point for looping while norm of f is decreasing * at an acceptable rate. */ L_50: ; /* Obtain the NEQ by NEQ partial matrix with respect to y' * Note that there is no check on the norm of f until at least one * update step for y'. This is to catch systems that are of * index .gt. 1, but the user has provided consistent initial data. */ ires = 7; /* Clear out the working matrix C(:,1:NEQ) before * computing the partials f_y'. */ if (dense) nrows = neq; if (banded) nrows = 2*ml + mu + 1; for (j = 1; j <= neq; j++) { for (i = 1; i <= nrows; i++) { C(j - 1,i - 1) = zero; } } if (irevc > 0) { (*ddasf)( t, y, yprime, &Rwork[iroff + 1], c, ldc, &cj, &ires, rwork, iwork ); } else { Info[1] = 2; goto L_380; } /* Branch back place if reverse communication in use. */ L_70: ; /* For the dense Jacobian case: * Compute LU factorization with full pivoting and explicit * interchanges. The use of full pivoting is to achieve the * best numerical stability and to put the small block, if * there is one, in the lower SE corner of the matrix. *---------------------------------------------------------------------- */ if (dense) { for (i = 1; i <= neq; i++) { /* Search in the remaining SE block for the largest magnitude. * Record the row and column where it is located. */ bigval = fabs( C(i - 1,i - 1) ); irow = i; jcol = i; for (j = i; j <= neq; j++) { for (k = i; k <= neq; k++) { if (fabs( C(j - 1,k - 1) ) <= bigval) goto L_80; bigval = fabs( C(j - 1,k - 1) ); irow = k; jcol = j; L_80: ; } } /* Record the row interchanges in IWORK(IOFF+1:IOFF+NEQ-1) * Record the column interchanges in IWORK(IOFF+NEQ+1:IOFF+2*NEQ) */ Iwork[i + ioff] = irow; Iwork[i + neq + ioff] = jcol; /* Interchange this column if necessary. Then make a rank * test and proceed on to the next phase if partials are * rank deficient. */ if (jcol > i) dswap( neq, &C(i - 1,0), 1, &C(jcol - 1,0), 1 ); if (irow > i) { temp = C(i - 1,i - 1); C(i - 1,i - 1) = C(i - 1,irow - 1); C(i - 1,irow - 1) = temp; } if (fabs( C(i - 1,i - 1) ) > rnktol*dasum( i, &C(i - 1,0), 1 )) { /* Compute multipliers and apply column operations to rest of matrix. */ temp = -one/C(i - 1,i - 1); dscal( neq - i, temp, &C(i - 1,i), 1 ); for (j = i + 1; j <= neq; j++) { temp = C(j - 1,irow - 1); /* Note that any interchange in a row is done just before * the column operation. */ if (irow > i) { C(j - 1,irow - 1) = C(j - 1,i - 1); C(j - 1,i - 1) = temp; } daxpy( neq - i, temp, &C(i - 1,i), 1, &C(j - 1,i), 1 ); } /* Apply to the right-hand side, f(t,y,y'). */ temp = Rwork[irow + iroff]; if (irow > i) { Rwork[irow + iroff] = Rwork[i + iroff]; Rwork[i + iroff] = temp; } daxpy( neq - i, temp, &C(i - 1,i), 1, &Rwork[i + 1 + iroff], 1 ); } else { /* Have rank deficiency. Proceed to index 1 phase. */ irank = i - 1; /* If IRANK == 0 the problem is an algebraic equation. * This exceptional case is flagged and a return is made. */ goto L_150; } } irank = neq; /* Have an index 0 problem. But it may not be consistent. * Back-solve for the update to y'. */ for (i = neq; i >= 1; i--) { temp = -Rwork[i + iroff]/C(i - 1,i - 1); daxpy( i - 1, temp, &C(i - 1,0), 1, &Rwork[iroff + 1], 1 ); Rwork[i + iroff] = -temp; } /* Rearrange solution corresponding to column pivoting. */ for (i = neq; i >= 1; i--) { j = Iwork[neq + i + ioff]; if (j > i) { temp = Rwork[j + iroff]; Rwork[j + iroff] = Rwork[i + iroff]; Rwork[i + iroff] = temp; } } /*---------------------------------------------------------------------- */ } else { /* Compute LU factorization of banded matrix. * The size of leading and trailing dimensions of * C(*,*) have been checked for adequate size. */ /* If there is room in the last rows of C(*,*) to save * f_y' then do so. This anticipates that the problem * may have index 1 or greater. Doing this copy saves * an evaluation of f_y' when the index is 1 or higher. */ if (*ldc > 4*ml + 2*mu + 3) { for (j = 1; j <= neq; j++) { dcopy( 2*ml + mu + 1, &C(j - 1,0), 1, &C(j - 1,ml*2 + mu + 1), 1 ); } } /* Save the current value of f(t,y,y') in * RWORK(2*NEQ+IROFF+1:3*NEQ+IROFF). This may be unnecessary if * the problem has index 0. If it has index 1 this saves an * evaluation of f(). The size of RWORK() has been checked. */ dcopy( neq, &Rwork[iroff + 1], 1, &Rwork[2*neq + iroff + 1], 1 ); dgbfa( c, *ldc, neq, ml, mu, &Iwork[ioff + 1], &iflag ); /* Check for a nearly rank deficient problem. * IF IFLAG .eq. 0 the code DGBFA did not encounter a zero diagonal. * A better check is to see if any diagonal terms are small * compared to a column L_1 norm. */ for (i = 1; i <= neq; i++) { if (fabs( C(i - 1,ml + mu) ) <= dasum( ml + mu + 1, &C(i - 1,0), 1 )*rnktol) { iflag = i - 1; /* Process the index 1 (or higher) case. This requires computing * df/dy and df/dt. */ goto L_150; } } irank = neq; /* There is no rank deficiency using the banded matrix * df/dy'. Compute the Newton step update and apply it */ dgbsl( c, *ldc, neq, ml, mu, &Iwork[ioff + 1], &Rwork[iroff + 1], 0 ); } /*---------------------------------------------------------------------- * Update the values of y' for index 0 problems. * Loop back to check size of f(t,y,y') again. */ daxpy( neq, -one, &Rwork[iroff + 1], 1, yprime, 1 ); /* This is an iteration counter. It is used * as a one-time check to allow one iteration before * checking the norm of f(t,y,y'). */ Iwork[5] += 1; /* Loop back to get another function evaluation * and check its norm against the user tolerance FTOL. */ goto L_20; /*----------------------------------------------------------------------- * This is the index 1 phase of solving for y'. */ L_150: ; if (dense) { /* Check if rank of f_y' .eq. 0. If so flag and return. * This gives the rank only when full pivoting was done. */ if (irank == 0) { *idid = -1; Iwork[1] = 5; goto L_370; } /* Copy the transformed value of f(t,y,y') and save it. */ dcopy( neq, &Rwork[iroff + 1], 1, &Rwork[neq + iroff + 1], 1 ); } /* The space in RWORK(:) is used for f_t. * Get the partial f_t. */ ires = 6; /* If IRES=1 or 6, clear out the spot where the results will go. * These are respectively the function and f_t. */ for (i = 1; i <= neq; i++) { Rwork[i + iroff] = zero; } if (irevc > 0) { (*ddasf)( t, y, yprime, &Rwork[iroff + 1], c, ldc, &cj, &ires, rwork, iwork ); } else { Info[1] = 3; goto L_380; } /* Branch back place if reverse communication in use. */ L_170: ; /* Get the partial f_y. This is placed in C(1:NEQ,1:NEQ) * so check that there is room. Require LDC .ge. 2*NEQ in the * dense Jacobian case and .gt. 4*ML+2*MU+3 in the banded case. */ if (dense) { if (*ldc < 2*neq) { *idid = -1; Iwork[1] = 8; Info[1] = 0; return; } } else { if (*ldc <= 4*ml + 2*mu + 3) { *idid = -1; Iwork[1] = 8; Info[1] = 0; return; } /* In the banded case the second dimension must * be at least 2*NEQ. */ if (ltd < 2*neq) { *idid = -1; Iwork[1] = 15; Info[1] = 0; return; } } ires = 8; /* Clear out rows of the working matrix before computing * the partials f_y. */ if (dense) nrows = neq; if (banded) nrows = 2*ml + mu + 1; for (j = 1; j <= neq; j++) { /* Copy the working store for the LU factorization into rows * NEQ+1:2*NEQ for dense problems. */ if (dense) { dcopy( neq, &C(j - 1,0), 1, &C(j - 1,neq), 1 ); } /* Zero the working rows in C(,) that will be used to * compute the matrix f_y. */ for (i = 1; i <= nrows; i++) { C(j - 1,i - 1) = zero; } } if (irevc > 0) { (*ddasf)( t, y, yprime, &Rwork[iroff + 1], c, ldc, &cj, &ires, rwork, iwork ); } else { Info[1] = 4; goto L_380; } /* Branch back place if reverse communication in use. */ L_190: ; /* Compute the matrix-vector operation: RWORK <- f_t + f_y*y' * Use a dense matrix-vector product operation. */ if (dense) { dgemv( 'N', neq, neq, one, c, *ldc, yprime, 1, one, &Rwork[iroff + 1], 1 ); } else if (banded) { for (j = 1; j <= neq; j++) { i1 = max( 1, j - mu ); i2 = min( neq, j + ml ); daxpy( i2 - i1 + 1, Yprime[j], &C(j - 1,i1 - j + ml + mu), 1, &Rwork[iroff + i1], 1 ); } } if (dense) { /* Apply the transformations and column pivoting to f_y. */ for (i = 1; i <= irank; i++) { /* Interchange columns as required. */ jcol = Iwork[i + neq + ioff]; if (jcol > i) dswap( neq, &C(i - 1,0), 1, &C(jcol - 1,0), 1 ); irow = Iwork[i + ioff]; for (j = 1; j <= neq; j++) { temp = C(j - 1,irow - 1); if (irow > i) { C(j - 1,irow - 1) = C(j - 1,i - 1); C(j - 1,i - 1) = temp; } if (i < neq) daxpy( neq - i, temp, &C(i - 1,neq + i), 1, &C(j - 1,i), 1 ); } /* Apply transformation to the right-hand side, f_t + f_y*y' */ temp = Rwork[irow + iroff]; if (irow > i) { Rwork[irow + iroff] = Rwork[i + iroff]; Rwork[i + iroff] = temp; } if (i < neq) daxpy( neq - i, temp, &C(i - 1,neq + i), 1, &Rwork[i + iroff + 1], 1 ); } /* Now use the IRANK by NEQ upper trapazoid of the transformed matrix * f_y' and the last NEQ-IRANK by NEQ block of the now transformed * matrix f_y. Eliminate using this special structure to complete the * factorization. The working space C(1:NEQ,1:NEQ) is used for * this last step. The resulting matrix should be of full rank but * may not be if the problem has index .gt. 1. */ for (i = 1; i <= neq; i++) { for (k = i + 1; k <= irank; k++) { /* Explicitly zero entries to achieve upper trapazoidal form. */ C(i - 1,k - 1) = zero; } /* Copy the upper trapazoidal matrix into working space. */ dcopy( min( i, irank ), &C(i - 1,neq), 1, &C(i - 1,0), 1 ); } /* Copy the first IRANK entries of the transformed f(t,y,y') * values into RWORK(:). */ dcopy( irank, &Rwork[neq + iroff + 1], 1, &Rwork[iroff + 1], 1 ); /* Use LINPACK code to compute LU factorization of the combined * transformed blocks. */ dgefa( c, *ldc, neq, &Iwork[ioff + 1], &iflag ); /* Check for a nearly rank deficient problem. * IF IFLAG .eq. 0 the code DGEFA did not encounter a zero diagonal. * A better check is to see if the diagonal terms are small * compared to a column L_1 norm. */ for (i = 1; i <= neq; i++) { if (fabs( C(i - 1,i - 1) ) <= dasum( i, &C(i - 1,0), 1 )* rnktol) { iflag = i - 1; *idid = -1; Iwork[1] = 4; goto L_370; /* Branch place for rank deficiency when it is not * expected for index 1 problems. This indicates * a value of the index greater than 1. */ } } /* Compute updates of y' using LINPACK code. */ dgesl( c, *ldc, neq, &Iwork[ioff + 1], &Rwork[iroff + 1], 0 ); /* Rearrange solution corresponding to column pivoting. * Note that only IRANK column permutations are recorded. */ for (i = irank; i >= 1; i--) { j = Iwork[neq + i + ioff]; if (j > i) { temp = Rwork[j + iroff]; Rwork[j + iroff] = Rwork[i + iroff]; Rwork[i + iroff] = temp; } } /* Update the values of y' for index 0 problems. * Loop back to check size of f(t,y,y') again. */ } else if (banded) { /* Compute the regularization parameter. If the user gives * it then for INFO(14) .gt. 0 it is passed in RWORK(INFO(14)). * If the user has not passed it then use macheps**(2./3.) for * the value. */ if (Info[14] > 0) { regul = Rwork[Info[14]]; } else { /* Get machine eps for this precision - */ eps = DBL_EPSILON; regul = pow(eps,2.e0/3.e0); } /* Construct the interleaved banded matrix from the banded * matrices f_y' and f_y. This requires use of the same * data structure that now contains [f_y * f_y'] to hold the * ultimate interleaved and regularized banded matrix. */ m = ml + mu + 1; mlp = 2*ml + 1; mup = max( 2*mu, 1 ); mp = mlp + mup + 1; for (j = neq; j >= 1; j--) { /* Copy the value J into columns 2*J-1 and 2*J of ultimate * matrix banded representation. */ /* Clear space in buffer area of RWORK(:). */ for (i = 1; i <= (2*neq); i++) { Rwork[3*neq + iroff + i] = zero; Rwork[5*neq + iroff + i] = zero; } /* Place columns 2*J-1 and 2*J of the interleaved matrix. * This involves values of column J from f_y, f_y' and the * regularization parameter, REGUL. Build the pair of columns * in the buffer area first. */ i1 = max( 1, j - mu ); i2 = min( neq, j + ml ); for (i = i1; i <= i2; i++) { is = i - j + m; /* Place values of f_y' into buffer space. */ Rwork[3*neq + iroff + 2*i - 1] = C(j - 1,is + ml*2 + mu); Rwork[5*neq + iroff + 2*i] = C(j - 1,is + ml*2 + mu); /* Place values of f_y into buffer space. */ Rwork[3*neq + iroff + 2*i] = C(j - 1,is - 1); } /* Place value of regularization parameter into buffer space. */ Rwork[5*neq + iroff + 2*j - 1] = regul; /* Clear target area for columns 2*J-1 and 2*J of interleaved * banded matrix. */ for (i = 1; i <= (2*mlp + mup + 1); i++) { C(j*2 - 2,i - 1) = zero; C(j*2 - 1,i - 1) = zero; } /* Copy from buffer area to target columns 2*J, then 2*J-1. */ i1 = max( 1, 2*j - mup ); i2 = min( 2*neq, 2*j + mlp ); for (i = i1; i <= i2; i++) { is = i - 2*j + mp; C(j*2 - 1,is - 1) = Rwork[5*neq + iroff + i]; } i1 = max( 1, 2*j - 1 - mup ); i2 = min( 2*neq, 2*j - 1 + mlp ); for (i = i1; i <= i2; i++) { is = i - (2*j - 1) + mp; C(j*2 - 2,is - 1) = Rwork[3*neq + iroff + i]; } } /* Interleave values of f and f_t + f_y * y' * These occupy array positions RWORK(NEQ+IROFF+1:3*NEQ+IROFF) * after the merge step. */ for (j = 1; j <= neq; j++) { /* Get f component into place */ Rwork[neq + iroff + 2*j - 1] = Rwork[2*neq + iroff + j]; /* Get f_t + f_y*y' component into alternate place */ Rwork[neq + iroff + 2*j] = Rwork[iroff + j]; } dgbfa( c, *ldc, 2*neq, mlp, mup, &Iwork[ioff + 1], &iflag ); /* Check for a zero diagonal. This should not happen for * index 1 problems. Declare the problem to have index .gt. 1. */ if (iflag > 0) { *idid = -1; Iwork[1] = 4; goto L_370; } dgbsl( c, *ldc, 2*neq, mlp, mup, &Iwork[ioff + 1], &Rwork[iroff + neq + 1], 0 ); for (j = 1; j <= neq; j++) { /* Get the changes for values of y', now interleaved */ Rwork[iroff + j] = Rwork[iroff + neq + 2*j - 1]; } } daxpy( neq, -one, &Rwork[iroff + 1], 1, yprime, 1 ); /* This is an iteration counter. It is used * as a one-time check to allow one iteration before * checking the norm of f(t,y,y'). */ Iwork[5] += 1; /* Loop back to get another function evaluation * and check its norm against the user tolerance FTOL. */ goto L_20; /* This block of code checks values of array sizes and other * parameters. If there are errors a return is made with * indicators flagged in IDID and IWORK(1). It is executed * just once per usage of ddasls. It is not re-entered when * reverse communication is used. */ L_360: ; if (banded) { /* For banded problems extract the lower and upper bandwidth. */ ml = Iwork[1]; mu = Iwork[2]; } Iwork[1] = 0; *idid = 0; irank = -1; if (neq <= 0) { Iwork[1] = 6; *idid = -1; goto L_370; } /* An additional requirement is that LDC .ge. 2*NEQ * for index 1 dense problems. This is checked as needed. */ if (dense) { if (*ldc < neq) { Iwork[1] = 7; *idid = -1; goto L_370; } } else { /* Test that bandwidth parameters are non-negative. */ if (ml < 0 || mu < 0) { Iwork[1] = 14; *idid = -1; goto L_370; } /* The row dimension of the working array must be at least * 2*ML+MU+1 for banded Jacobians. If this is an index 1 * problem the value must be larger, 2*(2*ML+1)*2*MU+1. * This is checked when an index 1 problem is encountered. */ if (*ldc <= 2*ml + mu) { Iwork[1] = 7; *idid = -1; goto L_370; } } /* In the dense or banded case the second dimension must * be at least NEQ. For banded index 1 problems the value must * be at least 2*NEQ. This is checked as needed. */ if (ltd < neq) { *idid = -1; Iwork[1] = 15; goto L_370; } /* Check sizes of additional work arrays. */ if (dense) { if (lrw < 2*neq + iroff) { Iwork[1] = 9; *idid = -1; goto L_370; } } else if (banded) { if (lrw < 7*neq + iroff) { Iwork[1] = 9; *idid = -1; goto L_370; } } if (liw < 2*neq + ioff) { Iwork[1] = 10; *idid = -1; goto L_370; } if (ftol <= zero) { Iwork[1] = 12; *idid = -1; goto L_370; } if (rnktol < zero) { Iwork[1] = 13; *idid = -1; goto L_370; } /* Set INFO(1)=0 so integration can start. First the problem * corresponding to IDID .lt. 0 must be fixed. */ L_370: ; if (*idid < 0) { /* Error messages can be put here. */ Info[1] = 0; return; } /* Save an initial dummy large value for the L_1 norm of f. * This is updated each iteration and is used to enforce * a minimal rate of decrease in the L_1 norm of f. */ Rwork[2] = DBL_MAX; goto L_10; /* This is the fixup and save place for reverse communication * function requests. All saved values are stored and the direction * flag, IRES, is stored for aiding the user in choosing * which evalution to make. No error processing here. */ L_380: ; Iwork[3] = ires; Iwork[6] = irank; Iwork[7] = ml; Iwork[8] = mu; *idid = 1; return; /* Last statement of routine ddasls */ #undef C } /* end of function */