SUBROUTINE AUGINP(OLDFIL,NEWFIL,NEEDSI,ERRCD,ERRMSG)
C
C FUNCTION:
CF
CF
C USAGE:
CU
CU
C INPUTS:
CI
CI
C OUTPUTS:
CO
CO
C ALGORITHM:
CA
CA
C MACHINE DEPENDENCIES:
CM
CM
C HISTORY:
CH
CH written by:
CH date:
CH current version:
CH modifications:
CH added dpcom: 7/16/88 jdb
CH
C ROUTINES CALLED:
CC
CC
C COMMON MEMORY USED:
CM
CM DPCOM -- see dpcommon.f and dpcom.f
CM
C----------------------------------------------------------------------
C written for: The CASCADE Project
C Oak Ridge National Laboratory
C U.S. Department of Energy
C contract number DE-AC05-840R21400
C subcontract number 37B-07685C S13
C organization: The University of Tennessee
C----------------------------------------------------------------------
C THIS SOFTWARE IS IN THE PUBLIC DOMAIN
C NO RESTRICTIONS ON ITS USE ARE IMPLIED
C----------------------------------------------------------------------
C
C
INCLUDE 'Parameter.f'
C
DOUBLE PRECISION A(SIZE,SIZE)
DOUBLE PRECISION B(SIZE,SIZE)
DOUBLE PRECISION C(SIZE,SIZE)
DOUBLE PRECISION D(SIZE,SIZE)
C
INTEGER INTEG(SIZE)
INTEGER NSTATS
INTEGER NINPS
INTEGER NOUTS
INTEGER ERRCD
C
CHARACTER*(*) OLDFIL
CHARACTER*(*) NEEDSI
CHARACTER*(*) NEWFIL
CHARACTER*(*) ERRMSG
C
INCLUDE 'dpcom.f'
C
ERRCD=0
C
CALL INSYS(OLDFIL,NINPS,NOUTS,NSTATS,
2 SIZE,SIZE,SIZE,A,B,C,D,ERRCD)
IF (ERRCD .NE. 0) THEN
ERRMSG = 'AUGINP: Fatal error from INSYS '//
1 'accessing '//OLDFIL
RETURN
END IF
CLOSE (UNIT=UNIT1)
C
C--read in vector of integers representing number of integrators to
C--be added at each input
C
OPEN (UNIT=UNIT1,FILE=NEEDSI,STATUS='OLD',ERR=999)
C
READ (UNIT1,*,END=999,ERR=999) (INTEG(I),I=1,NINPS)
C
CLOSE (UNIT=UNIT1)
C
C--check to see if the new system matrices will fit into arrays of
C--parameter size.
C
N = 0
DO 10, I = 1,NINPS
N = N + INTEG(I)
10 CONTINUE
C
IF (N+NSTATS .GT. SIZE) THEN
ERRCD = 1000 + N
ERRMSG = 'AUGINP: Fatal error; Addition of '//
2 'integrators exceeds parameter size.'
RETURN
END IF
C
C--for each input, add the appripriate number of integrators
C
DO 20, I = 1, NINPS
C
IF (INTEG(I) .NE. 0) THEN
C
C--move the ith column of b into the right side of a
C
CALL CPYCOL (SIZE,SIZE,NSTATS,B(1,I),A(1,NSTATS+1))
C
C--if adding more than one integrator, add integ(i)-1 columns of
C--zeros to a
C
IF (INTEG(I) .GT. 1) THEN
DO 30, J = NSTATS+2, NSTATS+INTEG(I)
DO 30, K = 1, NSTATS+INTEG(I)
A(K,J) = 0.D0
30 CONTINUE
END IF
C
C--zero out the lower rows of a; note that the last step took
C--care of columns nstats+2 thru nstats+integ(i)
C
DO 40, J = 1, NSTATS+1
DO 40, K = NSTATS+1, NSTATS+INTEG(I)
A(K,J) = 0.D0
40 CONTINUE
C
C--add the ones, if necessary (integ(i) .gt. 1), to the
C--lower right block of a
C
DO 50, J = NSTATS+2, NSTATS+INTEG(I)
A(J-1,J) = 1.D0
50 CONTINUE
C
C--now form new b
C--start by zeroing lower partition of b
C
DO 60, J = 1, NINPS
DO 60, K = NSTATS+1, NSTATS+INTEG(I)
B(K,J) = 0.D0
60 CONTINUE
C
C--zero ith column of b
C
DO 70, J = 1, NSTATS
B(J,I) = 0.D0
70 CONTINUE
C
C--and add a 1.d0 to the last element of the ith column
C
B(NSTATS+INTEG(I),I) = 1.D0
C
C--now form c; all that has to be done here is zero out the right partition
C
DO 80, J = NSTATS+1, NSTATS+INTEG(I)
DO 80, K = 1, NOUTS
C(K,J) = 0.D0
80 CONTINUE
C
C--finally, update n_states
C
NSTATS = NSTATS+INTEG(I)
C
C--end of processing of ith column
C
END IF
20 CONTINUE
C
C--write the new system to file new_file_spec; make sure default
C--extension is .SYS
C
CALL OUTSYS(NEWFIL,NINPS,NOUTS,NSTATS,
2 SIZE,SIZE,SIZE,A,B,C,D,ERRCD)
IF (ERRCD .NE. 0) THEN
ERRMSG = 'AUGINP: Fatal error from OUTSYS '//
1 'accessing '//NEWFIL
RETURN
END IF
CLOSE (UNIT=UNIT1)
C
C--end of processing
C
RETURN
C
C--eof error processing
C
999 ERRCD=1
ERRMSG = 'AUGINP: Fatal Error trying to access '//
1 'integrators file '//NEEDSI
RETURN
END