/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:45 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */
#include
#include "fcrt.h"
#include "sdasls.h"
#include
#include
void /*FUNCTION*/ sdasls(
void (*sdasf)(),
long neq,
float *t,
float y[],
float yprime[],
long info[],
float ftol,
float rnktol,
float *c,
long *ldc,
long ltd,
long *idid,
float 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;
float 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;
float *const Rwork = &rwork[0] - 1;
float *const Y = &y[0] - 1;
float *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 sdasls Krogh Avoided ref. to undefined IWORK locations.
*>>2008/10/26 sdasls Krogh Changed 'No' no 'N' in call to sgemv.
*>>2008/09/18 sdasls Krogh obtained from Hanson set up single prec.
*--S 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 SDASLX 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 SDASF written for SDASLX that separately
* computes f_t, f_y and f_y'. The evaluation of f(t,y,y') uses the
* same flag as SDASLX. */
/* SDASF: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 SDASF(T,Y,YPRIME,DELTA,C,LDC,CJ,IRES,RWORK,IWORK) */
/* The actions taken depend on the value of the IRES argument
* to SDASF. The values of IRES that occur when using sdasls 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 SDASLX.
* 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 SDASLX 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 sdasf 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 sdasls() 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 SDASF 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 R1MACH.
* 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 SDASF. Set to 0 for all SDASF 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 SDASF.
* This is not used but the value is assigned 0.
* ONE, ZERO - Local copies of 1.e0 and 0.e0. */
/* Integer Variables:
* IRES - A flag for directing user action in routine sdasf
* 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, SGEFA.
* 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;
(*sdasf)( 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 sdasf 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)
{
(*sdasf)( 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 = sasum( 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)
{
(*sdasf)( 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 = fabsf( C(i - 1,i - 1) );
irow = i;
jcol = i;
for (j = i; j <= neq; j++)
{
for (k = i; k <= neq; k++)
{
if (fabsf( C(j - 1,k - 1) ) <= bigval)
goto L_80;
bigval = fabsf( 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)
sswap( 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 (fabsf( C(i - 1,i - 1) ) > rnktol*sasum( i, &C(i - 1,0),
1 ))
{
/* Compute multipliers and apply column operations to rest of matrix. */
temp = -one/C(i - 1,i - 1);
sscal( 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;
}
saxpy( 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;
}
saxpy( 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);
saxpy( 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++)
{
scopy( 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. */
scopy( neq, &Rwork[iroff + 1], 1, &Rwork[2*neq + iroff + 1],
1 );
sgbfa( c, *ldc, neq, ml, mu, &Iwork[ioff + 1], &iflag );
/* Check for a nearly rank deficient problem.
* IF IFLAG .eq. 0 the code SGBFA 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 (fabsf( C(i - 1,ml + mu) ) <= sasum( 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 */
sgbsl( 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. */
saxpy( 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. */
scopy( 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)
{
(*sdasf)( 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)
{
scopy( 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)
{
(*sdasf)( 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)
{
sgemv( '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 );
saxpy( 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)
sswap( 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)
saxpy( 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)
saxpy( 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. */
scopy( 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(:). */
scopy( irank, &Rwork[neq + iroff + 1], 1, &Rwork[iroff + 1],
1 );
/* Use LINPACK code to compute LU factorization of the combined
* transformed blocks. */
sgefa( c, *ldc, neq, &Iwork[ioff + 1], &iflag );
/* Check for a nearly rank deficient problem.
* IF IFLAG .eq. 0 the code SGEFA 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 (fabsf( C(i - 1,i - 1) ) <= sasum( 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. */
sgesl( 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 = FLT_EPSILON;
regul = powf(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];
}
sgbfa( 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;
}
sgbsl( 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];
}
}
saxpy( 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 sdasls. 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] = FLT_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 sdasls */
#undef C
} /* end of function */